De novo Transcriptome Assembly of Floral Buds of Pineapple and Identification of Differentially Expressed Genes in Response to Ethephon Induction

A remarkable characteristic of pineapple is its ability to undergo floral induction in response to external ethylene stimulation. However, little information is available regarding the molecular mechanism underlying this process. In this study, the differentially expressed genes (DEGs) in plants exposed to 1.80 mL·L−1 (T1) or 2.40 mL·L−1 ethephon (T2) compared with Ct plants (control, cleaning water) were identified using RNA-seq and gene expression profiling. Illumina sequencing generated 65,825,224 high-quality reads that were assembled into 129,594 unigenes with an average sequence length of 1173 bp. Of these unigenes, 24,775 were assigned to specific KEGG pathways, of which metabolic pathways and biosynthesis of secondary metabolites were the most highly represented. Gene Ontology (GO) analysis of the annotated unigenes revealed that the majority were involved in metabolic and cellular processes, cell and cell part, catalytic activity and binding. Gene expression profiling analysis revealed 3788, 3062, and 758 DEGs in the comparisons of T1 with Ct, T2 with Ct, and T2 with T1, respectively. GO analysis indicated that these DEGs were predominantly annotated to metabolic and cellular processes, cell and cell part, catalytic activity, and binding. KEGG pathway analysis revealed the enrichment of several important pathways among the DEGs, including metabolic pathways, biosynthesis of secondary metabolites and plant hormone signal transduction. Thirteen DEGs were identified as candidate genes associated with the process of floral induction by ethephon, including three ERF-like genes, one ETR-like gene, one LTI-like gene, one FT-like gene, one VRN1-like gene, three FRI-like genes, one AP1-like gene, one CAL-like gene, and one AG-like gene. qPCR analysis indicated that the changes in the expression of these 13 candidate genes were consistent with the alterations in the corresponding RPKM values, confirming the accuracy and credibility of the RNA-seq and gene expression profiling results. Ethephon-mediated induction likely mimics the process of vernalization in the floral transition in pineapple by increasing LTI, FT, and VRN1 expression and promoting the up-regulation of floral meristem identity genes involved in flower development. The candidate genes screened can be used in investigations of the molecular mechanisms of the flowering pathway and of various other biological mechanisms in pineapple.


INTRODUCTION
Pineapple [Ananas comosus (L.) Merr.], which is indigenous to Brazil, Argentina and Paraguay, has been introduced worldwide and is the leading edible member of the Bromeliaceae family (Smith and Downs, 1974). It is an economically important tropical herbaceous fruit with attractive sensorial and nutritional characteristics (Coppens d'Eeckenbrugge and Leal, 2003). It is generally consumed fresh as table fruit or in desserts and is used for juice preparation due to its delicate flavor, overall acceptability and nutritional richness in terms of vitamins and minerals (Chauhan et al., 2009).
Flowering is a crucial developmental stage in the plant life cycle. The transition from vegetative growth to flowering is of importance as the first step of sexual reproduction and fruit setting (Bernier et al., 1993). A number of different factors, from environmental to chemical, can trigger flowering. Environmental conditions that promote natural flowering increase the sensitivity of plants to external flowering induction agents. Flowering induction can also be achieved via the use of a range of commercial agents. In pineapple and other bromeliads, it has been proposed that flowering is triggered by a small burst of ethylene production in the meristem in response to environmental cues (Trusov and Botella, 2006). This feature is exploited worldwide by commercial pineapple growers to synchronize flower and fruit development, resulting in a single harvest pass. This synchronization is of crucial importance in pineapple cultivation due to the non-climacteric nature of the fruit ( Van de Poel et al., 2009). Ethephon (2chloroethylphosphonic acid), an ethylene releaser, is universally accepted as the most effective stimulant to induce the flowering of pineapple plants and is used worldwide. Application of ethephon to the central cup induces more homogeneous flowering compared to application of ethylene gas to the entire plant (Van de Poel et al., 2009). With the use of ethephon for flower induction, pineapple can be fruited and harvested multiple times per year .
To date, few flowering genes have been cloned from pineapple. A homolog of the PISTILLATA (PI) gene, AcPI, which plays a crucial role in the regulation of flowering in angiosperms, was isolated from the pineapple cultivar "Comte de Paris". Its expression level was the highest at 40 days after flower induction, which is when multiple fruits and floral organs form (Lv et al., 2012a). A homolog of the FLOWERING LOCUS T (FT) gene, AcFT, was also isolated from pineapple. qPCR analysis revealed that its expression was elevated in the flesh and absent in the leaves and that it was expressed at the highest level at 40 days after flower induction (Lv et al., 2012b). In other studies of the mechanism of flowering stimulation by ethephon, three aminocyclopropane carboxylic acid (ACC) synthase genes (AcACSs), which play vital roles in the biosynthesis of ethylene, were isolated from pineapple (Botella et al., 2000;Trusov and Botella, 2006;Yuri and José, 2006;Choudhury et al., 2008). Transgenic pineapple plants carrying a silenced AcACS2 gene showed a marked delay in flowering compared to nonsilenced transgenic plants and control non-transformed plants. The AcACS2 gene may be one of the key contributors to the triggering of natural flowering in mature pineapple under commercial field conditions (Trusov and Botella, 2006;Yuri and José, 2006). Compared with other monocots, such as rice, wheat and cereal, little information is available regarding the role of flowering genes in the regulation of the vegetative to flowering transition and flower initiation in pineapple. In addition, few reports have investigated the genes related to the stimulation of pineapple flowering by ethephon and their functions.
Recently, Illumina sequencing techniques have enabled fascinating discoveries in the life sciences and have dramatically improved the efficiency of gene discovery (Zhang et al., 2013). RNA-seq technology allows for the determination of genomewide expression levels as well as the identification of novel transcripts and isoforms, and it has been used successfully in evaluations of numerous plant species, including Arabidopsis (Zhu et al., 2013), rice (Huang et al., 2014), soybean (Stamm et al., 2014), maize (Dukowic-Schulze et al., 2014;Thakare et al., 2014), and non-model species, such as wild strawberry (Sánchez-Sevilla et al., 2014). With regard to pineapple, Ong et al. (2012) performed transcriptome sequencing of ripe yellow pineapple fruit flesh using Illumina technology to determine the mechanisms and processes underlying fruit ripening. However, there is still a dearth of data on the molecular mechanisms of stimulation of pineapple flowering by ethephon. Further transcriptome sequencing studies of pineapple flowering and related genes induced by ethephon during the floral transition are needed.
Accordingly, in this study, RNA-seq and gene expression profiling were performed using Illumina technology to investigate the flowering pathway. In addition, various screened candidate genes related to flowering transition were analyzed by qPCR. The findings of this study contribute to the understanding of the molecular mechanisms of ethephon-mediated stimulation of the floral transition in pineapple.

Site Description
To investigate the molecular mechanisms of the ethephonmediated induction of the flowering transition in pineapple plants, a potted culture experiment was conducted at an experimental orchard at the Institute of Fruit Tree Research of the Guangdong Academy of Agricultural Sciences located in the Tianhe District (113.35 • N, 23.12 • E), Guangzhou, Guangdong Province, P. R. China. The local climate is a southern subtropical monsoon climate, with a mean annual precipitation of more than 1800 mm and an average annual air temperature of 21.5∼22.2 • C. The monthly average temperature is the highest between June and August, at approximately 28∼35 • C, and it is the lowest between December and February, at approximately 10∼15 • C.

Plant Materials and Treatments
In this study, the pineapple cultivar "Smooth Cayenne, " categorized as Cayenne, was used. In mid-October 2013, thirty pots containing pineapple plants with similar heights and leaf numbers were transported to a plastic shed. Five days later, flowering was induced in the pineapple plants with 100 mL of ethylene solution containing 1.80 mL·L −1 (T1) or 2.40 mL·L −1 (T2) ethephon (v/v 40%), which was applied to the central cups. The same volume of cleaning water was used as a control (Ct). Ten replicates of the two treatments and control were used, resulting in a total of 30 pots. Fifty days later, when the central cups of the flowering-induced plants opened, three shoot apical meristems (a total of 9) were randomly sampled from the T1, T2, and Ct groups, respectively. All of the samples were placed in liquid nitrogen promptly for further investigation.

RNA Extraction
Total RNA was isolated from each sample using a Trizol Kit (Promega, USA), following the manufacturer's instructions. Subsequently, the total RNA was treated with RNase-free DNase I (Takara Bio, Japan) for 30 min at 37 • C to remove residual DNA. RNA quality was verified using a 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA) and was also evaluated by RNasefree agarose gel electrophoresis. The concentration of total RNA was measured using a 2100 Bioanalyzer at 260 and 280 nm. In the process of RNA extraction, two additional extractions were performed for each treatment or control sample. Among the three RNA extractions from each treatment or control, those with a 260/280 nm ratio between 1.8 and 2.0 were considered to have the best quality and were used for subsequent analyses. Equal amounts of RNA from each sample were mixed for the subsequent steps.

RNA-seq Library Construction and Sequencing
Total RNA was isolated from each of the samples using a Trizol Kit (Promega, USA), following the manufacturer's instructions. Then, the total RNA was treated with RNase-free DNase I (Takara Bio, Japan) to remove residual DNA. Poly (A) mRNA was isolated using oligo-dT beads (Qiagen). Fragmentation buffer was added to the mRNA to generate short fragments (200 nt). First-strand cDNA was synthesized by random hexamer-primed reverse transcription, and then second-strand cDNA was generated using RNase H and DNA polymerase I. The cDNA fragments were purified using a QIAquick PCR Extraction Kit. These purified fragments were then washed with EB buffer for end reparation poly (A) addition and ligated to sequencing adapters. Following agarose gel electrophoresis and extraction of cDNA from the gels, the cDNA fragments (200 ± 25 bp) were purified and enriched by PCR to construct the final cDNA library. The cDNA library was sequenced with an Illumina sequencing platform (Illumina HiSeq ™ 2000) using paired-end technology. The processing of the original image to determine the sequences, in addition to the base calling and quality value calculations, were performed using Illumina GA Pipeline (version 1.6), from which 90 bp paired-end reads were obtained (Li et al., 2013).

Illumina Read Processing and Assembly
A Perl program was written to select clean reads by removing low-quality sequences (over 50% of bases with a quality score of lower than 20 in one sequence), reads with more than 5% N bases (unknown bases) and reads containing adaptor sequences. Then, the clean reads were assembled using Trinity to construct unique consensus sequences.

Determination of Gene Expression Levels
The sequencing reads were mapped to reference sequences using SOAPaligner/soap2 (Ernst and Bar-Joseph, 2006). Reads that could be uniquely mapped to a gene were used to calculate expression levels. Gene expression levels were measured according to the number of uniquely mapped reads per kilobase of exon region per million mapped reads (RPKM) using the following equation: RPKM = 10 6 C NL/10 3 in which C is the number of reads uniquely mapped to a given gene, N is the number of reads uniquely mapped to all genes, and L is the total length of the exons in the given gene. For genes with more than one alternative transcript, the longest transcript was selected to calculate the RPKM. The RPKM method eliminates the influences of different gene lengths and sequencing discrepancies on gene expression calculations. Therefore, RPKM values can be used directly for comparing differences in gene expression among samples.

Variations in Gene Expression
To identify the differentially expressed genes (DEGs) among T1, T2, and Ct, the number of clean tags for each gene was calculated, and the genes that were differentially expressed were identified based on the method described by Audic and Claverie (1997) for determining the threshold p-value, FDR and fold change (log Ratio 2 ). DEGs were defined as those with an FDR of ≤ 0.001 and an absolute value of log Ratio 2 ≥ 1.

Go and Pathway Enrichment Analyses of DEGs
To determine the main biological functions of the DEGs, they were first annotated using the Gene Ontology (GO) database (http://www.geneontology.org/) and Blast2GO (Conesa et al., 2005) according to their numerical order in the nr database. Blast2GO is an all-in-one tool for the functional annotation of novel sequences and analysis of annotation data. It has been cited by other articles over 150 times and is a widely recognized GO annotation software program. After obtaining GO annotations for each of the DEGs, WEGO software (Ye et al., 2006) was used to obtain GO functional classifications. GO enrichment analysis was performed using hypergeometric testing to identify significantly enriched GO terms among the DEGs compared to the genome background.

Validation of Subset of DEGs by qPCR
Total RNA was extracted as described above. Each RNA sample was treated with RNase-free DNase (Takara) following the manufacturer's protocol to remove any residual genomic DNA (gDNA). DNase-treated RNA (2 mg) was subjected to reverse transcription using arbitrary primers and PrimeScript ™ Reverse Transcriptase (Takara) according to the manufacturer's protocol. The gene expression levels were determined using Ct values with the formula 2 − Ct . The actin gene of A. comosus (HQ148720.1) was used as a reference gene. Each qPCR analysis was performed in triplicate. The sequences of the specific primer sets are listed in Additional File 1.

DGE Library Sequencing and Evaluation
Shoot apices were sampled from the Ct, T1, and T2 plants and were sequenced to construct three digital gene expression (DGE) libraries. The results of the sequencing quality assessment and alignment analysis are presented in Table 1. Low-quality reads, including adaptor reads and reads containing Ns, represented 1.27, 1.24, and 1.24% of the reads in the three libraries, respectively.

Sequence Assembly
To broadly investigate the genes associated with floral induction in pineapple, an overview of the pineapple transcriptome was obtained, which is summarized in Table 2. A cDNA library was prepared using a mixture of equal amounts of RNA isolated from the stem apices of Ct, T1, and T2 plants, and pairedend sequences were generated using an Illumina Hiseq ™ 2000 platform. After the cleaning and removal of low-quality reads (quality scores < 20), 65,825,224 clean reads with an average length of 100 bp were identified in the three libraries (66,656,448 raw reads), with a Q20 percentage (base quality score of over 20) of 98.03%. The GC percentage was 50.15%.
Based on the high-quality reads, 148,728 contigs were assembled, with 6,582,522,400 bp and an N50 of 2023 bp. The sizes and length distributions of the contigs and unigenes are presented in Figures 1, 2. The lengths of the contigs ranged from 201 to 11,724 bp, with an average length of 1225 bp. A total of 1225 contigs were >500 bp, and 2038 were between 400 and 500 bp in length. Most of the contigs were between 75 and 400 bp. Among them, 93,379 were between 300 and 400 bp in length, which represented 62.76% of the total contigs. However, 10,140 contigs (6.82%) were between 200 and 300 bp, 16, 014 (10.77%) were between 100 and 200 bp, and 29,195 (19.63%) were between 75 and 100 bp.
Using the paired-end reads, these contigs were further assembled into 129,594 unigenes by Trinity. The lengths of the unigenes ranged from 201 to 12,421 bp, with an average of 1173 bp.

Functional Annotation
The statistical results for annotation of the pineapple unigenes are summarized in Table 3. A search of the reference sequences and a BLASTX search against nr, SwissProt, COG, KEGG, GO and Pfam databases produced significant results for 129,594 unigenes.
The results of COG functional annotation of the pineapple unigenes are shown in Figure 3. A total of 35,523 unigenes were classified into 25 COG functional categories, among which the category "General function prediction" was predominant. A high percentage of genes were also assigned to "Transcription, " "Replication, recombination and repair, " and "Translation, ribosomal structure and biogenesis." Only a few genes were assigned to the categories "Extracellular structures" and "Nuclear structure" (Additional File 2).
COG analysis revealed that the identified genes are involved in various biological processes. For example, 5612 unigenes were annotated to "Signal transduction mechanisms" and are thus involved in signal transduction pathways, particularly   the transduction of plant hormones; in addition, 5194 were annotated to "Carbohydrate transport and metabolism" and are therefore associated with carbohydrate transport and metabolism pathways, and 5220 were annotated to "Cell cycle control, cell division, and chromosome partitioning" and are involved in cell division and chromosome partitioning. Pathways associated with the transduction of plant hormones, carbohydrate transport and metabolism, cell division, and chromosome partitioning are involved in floral induction. Thus, identification of the genes in these  pathways is crucial and will aid in the elucidation of the molecular mechanisms of floral induction by ethephon in pineapple.

GO Classifications of Unigenes
In this study, GO analysis was performed to categorize the functions of the predicted genes in pineapple, and the predominant GO categories among the assembled unigenes are presented in Figure 4. A total of 35,485 annotated genes were analyzed and divided into 44 functional GO categories within the three main ontologies of biological processes, cellular component, and molecular function. Cellular and metabolic processes were the most highly represented GO categories in the biological process ontology, with 18,640 and 18,627 unigenes, respectively. In addition, 6115 unigenes were annotated to the category "Response to stimulus, " and 2109 and 2066 unigenes were annotated to "Reproduction" and "Response to reproduction, " respectively. The GO annotations of these unigenes may aid in the identification of the molecular mechanisms of floral induction by ethephon in pineapple. The GO categories "Cell, " "Cell part, " and "Organelle" were predominant in the cellular component ontology, with 25,950, 25,950, and 21,578 unigenes, respectively. These unigenes may be associated with cell cycle control, cell division, and chromosome partitioning. The identification of these unigenes may be useful for understanding the molecular mechanisms of floral induction in pineapple by ethephon.

KEGG Pathway Assignments
The predominant KEGG pathways of the assembled unigenes are presented in Figure 5. In this study, 65,535 annotated sequences were mapped to reference pathways described in KEGG. In total, 24,775 unique sequences were assigned to 122 KEGG pathways (Additional File 3). Among the unigenes, the highest percentage (6081/24,775, 24.54%) was assigned to "Metabolic pathways, " followed by "Biosynthesis of secondary metabolites" (2857/24,775, 11.53%). These pathway assignments provide a valuable resource for investigating specific processes, functions, and pathways in pineapple.

Identifying Pineapple Flowering-Associated Genes
Based on the annotations, a few pineapple genes associated with flowering time were identified, including photoperiod pathway genes, such as EARLY FLOWERING 3 (ELF3), CIRCADIAN TIMEKEEPER (CTK), and COSTANS (CO); vernalization pathway genes, such as VERNALIZATION (VRN1), FRIGIDA (FRI), and MADS-box genes; floral integrator pathway genes, such as CAULIFLOWER (CAL), AP2, FT, and TERMINAL FLOWER 1 (TFL1); and the floral meristem identity genes LFY, AGMOUS (AG), and AP1. These results provide important insights into the mechanism of ethephon-mediated floral induction, floral development, and flower organ formation in pineapple.

Variations in Gene Expression
Variations in gene expression were examined by comparing the T1 and Ct, T2 and Ct, and T1 and T2 groups. The numbers of DEGs in each comparison are presented in Figure 6. The comparison of T1 and Ct, as well as that of T2 and Ct, revealed significant differences in gene expression. A total of 3788 DEGs were identified in T1 compared to Ct, including 1402 upregulated and 2386 down-regulated genes. In addition, a total of 164 DEGs were specifically expressed in T1, and 211 were specifically expressed in Ct. Further, 3062 DEGs were identified in T2 compared to Ct, including 1321 up-regulated and 1741 down-regulated DEGs. A total of 77 DEGs were specifically expressed in T2, and 116 were specifically expressed in Ct. Finally, 758 DEGs were detected in T2 compared to T1. Among them, 438 were up-regulated, and 320 were down-regulated. A total of 24 DEGs were specifically expressed in T2, and 105 were specifically expressed in T1. The DEGs that were specifically expressed in either the T1, T2 or Ct group are listed in Additional Files 4-6.

GO Classifications of DEGs
The GO classifications of the top DEGs are presented in Figure 7. In the biological process ontology, most of the DEGs were annotated to the metabolic process and cellular process GO categories (Figure 8). In total, 607 and 575 DEGs were annotated to metabolic process and cellular process, respectively, in T1 compared to Ct, in addition to 474 and 460, respectively, in T2 compared to Ct, and 133 and 121, respectively, in T2 compared to T1.
In the cellular component ontology, most of the DEGs were annotated to the categories cell, cell part and organelle (Figure 7). In total, 844, 844, and 687 DEGs were annotated to cell, cell part and organelle, respectively, in T1 compared to Ct, in addition to 620, 620, and 488, respectively, in T2 compared to Ct, and 182, 182, and 152, respectively, in T2 compared to T1.
In the molecular function ontology, most of the DEGs were annotated to the catalytic activity and binding categories (Figure 7). In total, 586 and 503 DEGs were annotated to catalytic activity and binding in T1 compared to Ct, in addition to 489 and 391, respectively, in T2 compared to Ct, and 121 and 113, respectively, in T2 compared to T1.

Candidate DEGs Related to Pineapple Floral Induction by Ethephon and qPCR Validation
Thirteen DEGs were identified as candidate DEGs involved in the process of floral induction by ethephon in pineapple. Figure 8 and Table 4 summarize the expression of these thirteen DEGs.

ERF-like and ETR-like
Three differentially expressed ERF-like (ethylene-responsive factor) genes, namely Unigene0021445, Unigene0115271, and Unigene0070140, were isolated and identified as candidate genes (Figure 8 and Table 4). Among them, Unigene0021445 and Unigene0115271 were identified in the comparison of T1 with Ct and were found to be up-regulated by 1.36-and 1.13-fold, respectively. Unigene0070140 was identified in the comparison of T2 with Ct and was found to be up-regulated by 2.54-fold.
A differentially expressed ETR-like (ethylene receptor) gene, Unigene0063713, was identified as a candidate gene in the T2 group and was found to be up-regulated by 1.75-fold (Figure 8 and Table 4). This gene was also up-regulated by 1.51-fold in the T1 group but was not considered to be differentially expressed because of the high FDR value (0.001977). These results indicate that floral stimulation by ethephon induces expression of the ETR gene.

LTI-like
In this study, a differentially expressed LTI-like (low temperature-induced) gene, Unigene0089653, was isolated and identified as a candidate gene. The expression of this gene was induced by the T1 and T2 treatments, which resulted in 12.74-and 12.25-fold up-regulation, respectively (Figure 8 and Table 4). These results showed that the treatment of pineapple with ethephon induced the expression of an LTI-like gene.

FT-like
Three FT-like genes were isolated and identified (Unigene0109587, Unigene0032386, and Unigene0032387), and their expression was found to be up-regulated in the T1 and T2 groups compared with the Ct group. In particular, the FT-like gene Unigene0032387 was identified as a DEG in the comparison of T2 with Ct, exhibiting 12.08-fold increased expression, and it was screened as a candidate gene (Figure 8 and Table 4).

VRN1-like
A VRN1-like gene, Unigene0077614, was isolated and identified as a DEG. Ethephon treatment induced the expression of this gene, which exhibited 15.69-and 16.55-fold up-regulation in the T1 and T2 groups, respectively (Figure 8 and Table 4); thus, it was screened as a candidate gene.

AP1-Like
An AP1-like gene, Unigene0058964, was isolated and identified as a DEG in the comparison of T2 with Ct. Its expression was up-regulated by 14.73-fold after ethephon stimulation in the T2 group (Figure 8 and Table 4).

CAL-Like
A CAL-like gene, Unigene0132144, was isolated and identified as a DEG. Its expression was up-regulated in the T1 and T2 groups by 11.41-and 11.30-fold, respectively (Figure 8 and Table 4).

AG-Like
One AG-like gene, Unigene0077607, was isolated and identified as a DEG. Its expression was up-regulated in the T1 and T2 groups by 3.52-and 4.89-fold, respectively (Figure 8 and Table 4).
To validate the accuracy and reproducibility of the transcriptome analysis results, qPCR was performed to assess the thirteen DEGs described above. Figure 9 shows the expression levels of the candidate DEGs determined using qPCR and RNA-seq. The expression profiles of the thirteen genes revealed

Illumina Sequencing of Pineapple and Sequence Annotation
The pineapple is an economically important tropical fruit with a delicate taste, and it has widespread consumer acceptance in both fresh and processed forms (Pino and Queris, 2010). Despite its economic importance, little research has been conducted on the genetics of this crop fruit because of the low number of pineapple genome sequences that are available in public databases. In a previous study, Ong et al. (2012) generated 4.7 million Illumina paired-end reads that were then assembled into 28,728 unique transcripts. This group performed NCBI searches and GO and KEGG pathway analyses and revealed the complete transcriptomic profile of pineapple.
In this study, RNA-seq and differential gene expression profiling analyses were performed using Illumina sequencing, which generated 65,825,224 clean reads that were assembled into 129,594 unigenes with an average sequence length of 1173 bp. The numbers of high-quality reads and assembled unigenes surpassed those of the previous study (Ong et al., 2012). The results of this study are strong, as the pineapple unigenes were assembled using only shoot apical meristem tissue, and the unigenes were assembled without the use of a reference genome.
The unigenes assembled in this study were subjected to a BLASTx similarity search and annotation against the nr, SwissProt, COG, KEGG, GO, and Pfam databases. A total of 19.12-58.33% of the unigenes were annotated, suggesting that the flower transition of pineapple involves many unique processes and pathways. Similar results have been found in studies of the transcriptome of an Orchidaceae species associated with floral development and the response of cassava to cold stress Zhang et al., 2013). The COG functional annotation of the unigenes revealed that the category "General function prediction" was predominant among the 25 functional categories identified. A large number of unigenes were also annotated to the category "Signal transduction mechanisms." Similar results have been obtained in studies of genes associated with floral development in Orchidaceae and Bambusoideae (Zhang et al., , 2013. In addition, a large number of unigenes were annotated to the categories "Carbohydrate transport and metabolism" and "Cell cycle control, cell division, and chromosome partitioning." Pathways involving carbohydrate transport and metabolism and cell division and chromosome partitioning are involved in plant floral induction (Rady and El-Yazal, 2013;Vaddepalli et al., 2015). These results provide a foundation for further studies of gene expression and the genomics of pineapple in relation to flower induction and reproductive growth.

Differential Gene Expression Profiling Related to Flower Induction
In this study, DEGs between each of the treatments and control were identified. More DEGs were identified in the comparisons of T1 and T2 with Ct, and only a few were detected in the comparison of T1 with T2. Most of the enriched KEGG pathways identified in the comparison of T1 with Ct were similar to those detected in the comparison of T2 with Ct. This result suggests that the two ethephon treatments (T1 and T2) have common functions in the induction of pineapple flowering. GO analysis of the transcriptomes revealed that the DEGs were predominantly annotated to subcategories such as metabolic process, cellular process, cell, cell part, organelle, and catalytic activity and binding. Similar results have been obtained in studies of differential gene expression in an Orchidaceae species in response to floral development, as well as in litchi in response to reactive oxygen species (Zhang et al., 2013;Lu et al., 2014). KEGG pathway analysis revealed that the pathways "Metabolic pathways" (ko01100), "Biosynthesis of secondary metabolites" (ko01110), and "Plant hormone signal transduction" (ko04075) were significantly enriched among the DEGs. Similar results have been obtained in studies of differential gene expression associated with floral development in Bambusoideae species, Paphiopedilum orchid, longan, and Tapiscia sinensis Lai and Lin, 2013;Li D. M. et al., 2015;Zhou et al., 2015). The DEGs assigned to the pathway "Plant hormone signal transduction" are involved in the signal transduction of several plant growth regulators, auxin, cytokinin, gibberellin, abscisic acid, ethylene, brassinosteroid, jasmonic acid (JA), and salicylic acid (SA) and thus regulate plant floral induction, reproductive growth and stimulus responses (Aloni et al., 2006;Li L. F. et al., 2015;Vaddepalli et al., 2015). These results indicate that ethephonmediated stimulation of floral induction in pineapple activates "Plant hormone signal transduction" pathways. This information provides a foundation for further identification of pineapple floral induction-associated genes.

Response of Pineapple Floral Induction-Associated Genes to Stimulation by Ethephon
In general, pineapple flowering induction is triggered by a small burst of ethylene production (ethephon solution) applied to the central cup (Trusov and Botella, 2006; Van de Poel et al., 2009). ACS (ACC synthase), which is the rate limiting enzyme in the ethylene biosynthesis pathway, is induced by exogenous application of ethylene in pineapple, and the AcACS2 gene is thought to be one of the key contributors to the triggering of natural flowering in mature pineapple (Botella et al., 2000;Trusov and Botella, 2006;Yuri and José, 2006;Choudhury et al., 2008). In this study, an ACS-like gene of pineapple was isolated, and RPKM analysis indicated that this gene was induced by ethephon stimulation. However, it was not considered to be differentially expressed because of the high P-value and FDR value. Further, three ERF-like (ethylene-responsive factor) genes and one ETR-like (ethylene receptor) gene were found to be differentially expressed, exhibiting up-regulation in response to ethephon stimulation. These results are in agreement with those of previous studies of the expression of ETR genes in cut rose during flower opening (Ma et al., 2006) and of the expression of ERF genes in rice during flowering (Hu et al., 2008) in response to exogenous ethylene. It is known that the expression of certain ERF genes is induced by ethylene, jasmonic acid (JA), abscisic acid (ABA), and environmental stress, suggesting that ERF proteins could play important roles in the processes regulated by these signaling pathways (Fujimoto et al., 2000;Gutterson and Reuber, 2004;Jin and Liu, 2008;Shinshi, 2008). Auxin has also been reported to promote ethylene synthesis, and both are related to the switch of pineapple from the vegetative to the flowering stage (Burg and Burg, 1966;Min and Bartholomew, 1997;Van de Poel et al., 2009). In this study, some auxin-related genes, such as genes encoding auxin-induced protein, auxin-responsive protein, and auxin transporter-like protein, were isolated and were found to be differentially expressed (Additional Files 4, 5).
These findings support those of the abovementioned previous studies.
The floral transition is a vital process for completion of the life cycles of higher plants. In Arabidopsis, previous molecular genetic studies have elucidated five genetically defined pathways that control flowering: the vernalization, photoperiod, autonomous, gibberellin, and age pathways (Srikanth and Schmid, 2011;Klepikova et al., 2015). These pathways converge to regulate the expression of several floral pathway integrators, including FT, SOC1/AGL20, LFY, and FLC (Boss et al., 2004;Parcy, 2005;Lei et al., 2013). They activate the downstream floral meristem identity genes AP1, CAL, and FUL, which in turn initiate the transition from the vegetative to the reproductive stage (Tani et al., 2009;Yu et al., 2013).
Vernalization, which can be defined as acceleration of flowering in plants, is promoted by a prolonged period of low temperatures (as typically experienced during the winter season), and it is considered to be one of the most specific primary factors for inducing and controlling reproductive development in some plant species (Tarkowská et al., 2012). In pineapple, ethylene is thought to be responsible for natural flowering at low temperatures (Botella et al., 2000). Nevertheless, there are few reports of genes related to low-temperature floral induction. In this study, an LTI-like (low temperature-induced) gene was found to be differentially expressed, exhibiting upregulation in response to the two ethephon treatments. Recent molecular genetic approaches to understanding the regulatory pathway that promotes flowering in Arabidopsis have led to the conclusion that the main essential component of florigen is FLOWERING LOCUS T (FT) protein. FT-like genes are crucial regulators of flowering in angiosperms. In different plant species, FT homologs are involved in the earliest stages of flower development (Zeevaart, 2006). In this study, three FT-like genes were isolated from the plants exposed to the T2 treatment (the higher concentration of ethephon). In addition, one FT-like gene was found to be differentially expressed in the comparison of T2 with Ct. Similar results have been observed in a previous study (Lv et al., 2012b).
In wheat, barley and other monocotyledons species, the vernalization response is controlled by a MADS-box gene, VERNALIZATION1 (VRN1; Kim et al., 2009;Alonso-Peral et al., 2011;Hemming and Trevaskis, 2011;Oliver et al., 2013;Yu et al., 2014). The VRN1 gene of monocotyledons such as cereals plays a crucial and dual role in flowering; first, it induces expression of the cereal FT homolog through the vernalization pathway, and second, it acts as a floral meristem identity gene (Li and Dubcovsky, 2008). In this study, a VRN1like gene was isolated from pineapple (a monocotyledon) and was found to be differentially expressed, exhibiting marked up-regulation in response to treatments T1 and T2 compared to Ct. These results indicate that the genes related to vernalization in pineapple are more closely related to vernalization genes in other monocotyledons than to those in dicotyledonous species and that ethephon induces expression of the VRN1 gene.
In Arabidopsis and beet, the vernalization process promotes flowering by repressing the gene FLOWERING LOCUS C (FLC). It is known that the expression of this gene is increased by the action of the FRIGIDA (FRI) gene. The FRI and FLC genes, in contrast, act as constitutive floral inhibitors (Schmitz and Amasino, 2007;Crevillen and Dean, 2011;Pin et al., 2012;Yruela, 2015). In this study, three FRI-like genes were isolated from pineapple and were found to be differentially expressed, exhibiting down-regulation in response to the T1 and T2 ethephon treatments. These results indicate that stimulation with ethephon results in down-regulation of the expression of FRI genes.
The results of this study, including the observed up-regulation of the LTI-like and VRN1-like genes and down-regulation of the FRI-like genes, suggest that stimulation by ethephon may mimic the vernalization process induced by low temperatures in pineapple plants. Notably, further study is needed to elucidate the specific mechanism of ethephon-mediated stimulation rather than vernalization (low temperatures) in induction of flowering in pineapple.
After flowering initiation, floral meristems require the action of several FMI genes, such as LFY, AP1, AG, and CAL, the upregulation of which controls floral organ identity and initiates the generation of floral meristems Kaufmann et al., 2010;Han et al., 2014;Jung et al., 2014;Li D. M. et al., 2015). The AP1, CAL, and FUL genes act redundantly to control inflorescence architecture by influencing LFY and TERMINAL FLOWER 1 (TFL1) expression as well as their relative activity levels (Ferrándiz et al., 2000). Stamens and carpels are converted to petaloid organs, which is related to AGAMOUS (AG) gene expression (Akita et al., 2011;Sharifi et al., 2015). AP1-like genes specify floral meristem and organ identities (Shulga et al., 2011). CAL partially functions in the regulation of floral meristem identity and has a unique role in regulating cell differentiation during fruit development (Ferrándiz et al., 2000). These floral meristem identity genes cooperate to promote the transition from vegetative growth to reproductive growth (Zhang et al., 2013). In this study, three floral meristem identity genes, including AP1-like, AG-like, and CAL-like, were found to be differentially expressed in the comparison of T2 with Ct. Two of them, AGlike and CAL-like, were found to be differentially expressed in the comparison of T1 with Ct. The expression of these genes was upregulated in response to stimulation by ethephon. These results suggested that stimulation by ethephon initiated the process of floral meristem identity and converted stamens and carpels to petaloid organs in pineapple by inducing specific genes. It also revealed the vital roles of the AP1-like, AG-like, and CAL-like genes in floral organ development and indicated that stimulation by ethephon induced the flowering of pineapple, perhaps by causing the increased expression of these flowering-related genes when it reaches the shoot apex (Turnbull et al., 1999). This information regarding flowering integration genes and FMI genes will facilitate more detailed studies of the mechanism of floral differentiation in pineapple.
In terms of the efficiencies of the two different concentrations of ethephon used for floral induction in this study, T1 and T2 induced floral transition in similar manners (Additional File 10) by up-regulating related genes. However, the T2 treatment, which contained the higher concentration of ethephon, solely up-regulated the FT-like and AP1-like genes, confirming that the efficiency of ethephon in floral induction in pineapple is determined by the ethephon concentration, the precise plant maturity level and climatic factors (Van de Poel et al., 2009).
The findings of this study provide a valuable resource for future pineapple genomic studies, and they can be applied in parallel studies of other closely related species. Notably, further studies are needed to verify the functions of the candidate genes identified in this study in floral induction in pineapple.

CONCLUSION
In this study, we performed RNA-seq and gene expression profiling analyses of pineapple exposed to ethephon for stimulation of flowering to determine the genome-wide gene expression profile. An abundance of specific pineapple genes related to different biological functions were identified. Many DEGs were isolated and identified in the comparisons of the T1 and T2 ethephon treatments with Ct. The GO annotations of the DEGs revealed that the metabolic process and cellular process categories of the biological process ontology were predominant, as well as cell, cell part and organelle in the cellular component ontology, and catalytic activity and binding in the molecular function ontology. KEGG pathway analysis indicated that metabolic pathways, biosynthesis of secondary metabolites, and plant hormone signal transduction were enriched among the DEGs.
Thirteen DEGs were screened as candidate genes involved in ethephon-mediated stimulation and floral induction as well as floral meristem identity, including the ERF-like, ETR-like, LTIlike, FT-like, VRN1-like, FRI-like, AP1-like, CAL-like, and AGlike genes. The qPCR data obtained for these thirteen DEGs were in agreement with the corresponding RPKM values derived from RNA-seq, confirming the accuracy and credibility of the RNA-seq and gene expression profiling results. With regard to the comparison of T2 with Ct, two up-regulated DEGs, FTlike gene and AP1-like, were solely identified. Considering the screened candidate genes with increased expression, it appears that ethephon-mediated stimulation mimics the vernalization process of flowering induction in pineapple and that it promotes the expression of floral meristem identity genes involved in flower development.
The DEGs identified in this study are good candidates for functional analyses of pineapple flowering induced by stimulation with ethephon. The results of this study are valuable for future pineapple genomic studies, and they can be applied in parallel studies of other closely related species. However, further studies are needed to verify the functions of these candidate genes.

AUTHOR CONTRIBUTIONS
CL carried out the experiments, performed the bioinformatics analyses, and wrote the manuscript. CF participated in the design and sampled the experimental materials.