De novo transcriptome assembly for rudimentary leaves in Litchi chinesis Sonn. and identification of differentially expressed genes in response to reactive oxygen species

Background Litchi is an evergreen woody tree widely cultivated in subtropical and tropical regions. Defective flowering is a major challenge for litchi production in time of climate change and global warming. Previous studies have shown that high temperature conditions encourage the growth of rudimentary leaves in panicles and suppress litchi flowering, while reactive oxygen species (ROS) generated by methyl viologen dichloride hydrate (MV) promote flowering and abortion of rudimentary leaves. To understand the molecular function of the ROS-induced abortion of rudimentary leaves in litchi, we sequenced and de novo assembled the litchi transcriptome. Results Our assembly encompassed 82,036 unigenes with a mean size of 710 bp, and over 58% (47,596) of unigenes showed significant similarities to known sequences in GenBank non-redundant (nr) protein database. 5,865 unigenes were found to be differentially expressed between ROS-treated and un-treated rudimentary leaves, and genes encoding signaling components of plant hormones such as ABA and ethylene were significantly enriched. Conclusion Our transcriptome data represents the comprehensive collection of expressed sequence tags (ESTs) of litchi leaves, which is a vital resource for future studies on the genomics of litchi and other closely related species. The identified differentially expressed genes also provided potential candidates for functional analysis of genes involved in litchi flowering underlying the control of rudimentary leaves in the panicles. Electronic supplementary material The online version of this article (doi:10.1186/1471-2164-15-805) contains supplementary material, which is available to authorized users.


Background
Litchi is one of the most important subtropical evergreen fruit trees in southern Asia. A major factor determining litchi crop production is the competition between vegetative and reproductive growth during floral differentiation. Floral initiation in litchi could be triggered by low temperatures and enhanced by droughts in autumn and winter [1,2]. In the following spring, the apical buds will break and elongate as the air temperature and soil moisture rise. Next, the axillary or apical panicle primordia will emerge and become visible as "whitish millets" [3].
At this millet stage, floral buds are considered to be mixed buds containing axillary or apical panicle primordia, leaf primordia and rudimentary leaves. Whether these mixed buds could develop into flowers are largely dependent on the environmental conditions. Under normal climate conditions, the growth of panicle primordia will prevail and the rudimentary leaves will abscise. However, if the buds are exposed to high temperature, the rudimentary leaves could develop into fully expanded leaves and the panicle primordia will cease to develop and shrink [4]. Suppressing the growth of the rudimentary leaves encourages panicle development. Warm winter and hot spring potentially resulted from global warming present a major threat to litchi flowering. In order to develop a counter measure, it is important to understand the genetics behind litchi bud development, and the abortion of the rudimentary leaves.
We have previously shown that ethylene could promote abortion of rudimentary leaves. It is associated with an increase in H 2 O 2 [5], a type of reactive oxygen species (ROS). It is well known that environmental stresses stimulate ROS production in plants [6]. They act as key signals in response to stresses [7]. In Arabidopsis, H 2 O 2 has been shown to be elevated in leaves at the time of floral transition [8]. It has been proposed that H 2 O 2 bursts generate a signal in leaves associated with either the induction of flowering or leaf senescence [9]. The ethylene-induced H 2 O 2 might act as a signal which induces abortion of rudimentary leaves in the panicle and promotes flowering in litchi. Methyl viologen dichloride hydrate (MV) is a ROS producer in plants. It can generate superoxide by accepting an electron from PSI to become a reduced free radical, which is immediately reoxidized by dioxygen, producing superoxide in chloroplast [10]. MV also induces the increase of superoxide production in mitochondria, where complexes I and III are the major electron donors [11]. Superoxide is then transformed to H 2 O 2 by superoxide dismutase [12,13]. When the litchi leafy panicles were treated with MV, it was found that ROS accumulated, the rudimentary leaves abscised and the numbers of flowers per panicle increased [14]. We have also shown that ROS increased the expression of LcLFY and LcAP1 in panicles [14,15]. Studies on Arabidopsis and other plants indicated that LFY (LEAFY) is a transcription factor which determines the floral meristem identity and is strongly expressed in the flower buds [16,17]. Constitutive expression of Arabidopsis AP1 (APETALA1) has also been shown to promote flowering in citrus [18]. AP1 is involved in the transition from floral induction to flower formation and constitutes a hub in the corresponding network of regulatory genes [19,20]. Beside a few ROS responsive EST clones derived from a suppression subtractive hybridization (SSH) library screen [21], little is known about the transcriptional network controlling litchi flowering.
Without a litchi reference genome, de novo transcriptome assembly using Illumina short RNA-Seq reads is the most cost effective approach for generating a large collection of ESTs suitable for subsequent transcriptome analysis. This method has been successfully applied to Chinese bayberry (Myrica rubra) and watermelon (Citrullus lanatus (Thunb.) Matsum. & Nakai var. lanatus) for fruit development and ripening studies [22,23], pear (Pyrus pyrifolia) for bud dormancy analysis [24,25], and litchi (Litchi chinesis) and melon (Cucumis melo) for fruit abscission study [26,27]. Though litchi fruit transcriptome sequencing data were published [27], those of leaves are unknown. In this study, we have constructed a litchi reference transcriptome for rudimentary leaves using Illumina RNA Sequencing (RNA-seq). We also used digital gene expression assay to profile the transcriptome dynamics of ROS treated rudimentary leaves, in order to elucidate genetic network of the ROS-induced abortion.

MV induced abortion of rudimentary leaves
In our previous study, we found that an early sign of abortion of rudimentary leaves was downward growth of the leaves [4]. To confirm the effect of MV treatment, shoot cuttings were treated with water or MV in a growth chamber. The proximal angle (α) and the distal angle (β) of the third and the fourth MV-treated rudimentary leaves were measured ( Figure 1A). We used this experimental system for the easy control of light intensity and temperature. Furthermore, in our preliminary study, it was confirmed that detached shoots could survive in water without wilting for at least 2 d. To avoid the immediate effect of the cutting from trees, treatments were carried out after the shoots were placed in water for 2 h. The results showed that proximal angle α had no significant change during the treatments, while the distal angle β significantly increased after the treatment ( Table 1). The ROS-treated rudimentary leaves showed epinasty as characterized by downward curvature of leaves ( Figure 1B, C, D), presenting an early sign of abortion [4].

Sequence, assembly and annotation of a litchi reference transcriptome
To obtain a reference litchi transcriptome for the rudimentary leaves, a RNA-Seq library has been constructed using RNA from all leaf samples. As shown in Table 2, we have generated 6.13 G of total nucleotides with a Q20 percentage of 97.8%. The Trinity package assembled 82,036 unigenes with a mean size of 710 bp ( Table 3). The size distributions of these unigenes are shown in Additional file 1. 58% (47,596/82,036) of unigenes could be annotated by BLASTx (E-value < 1e −5 ) using the NCBI nr database, while 34,368 were annotated using the Swiss-Prot protein database. In addition, 13,728 and 16,700 unigenes could be annotated according to the Kyoto Encyclopedia of Genes and Genomes (KEGG) and Cluster of Orthologous Groups of protein (COG) database, respectively. About 10% (8,191/82,036) unigenes could be assigned to a homolog in all four databases ( Figure 2A). Based on the NCBI nr database, 23.0% of the unigenes showed homology (1e −20 < E-value < 1e -5 ), 50.0% of those showed strong homology (1e −100 < Evalue < 1e −20 ) and the remaining 27.0% were very strong homology (E-value < 1e −100 ) to available plant sequences ( Figure 2B). As shown in Figure 2C, 28,556 unigenes were annotated to 4 top-hit species, including Glycine max, Arabidopsis thaliana, Medicago truncatula and Populus trichocarpa. 23,278 unigenes could be classified into 3 gene ontology (GO) categories: cellular component, biological process, and molecular function (Additional file 2). 13,728 unigenes of the rudimentary leaves in litchi were mapped into 274 KEGG pathways (Additional file 3). The maps with the highest unigene representation were ribosome pathway (ko03010) with 629 ungenes counted, followed by protein processing in endoplasmic reticulum (ko04141), starch and sucrose metabolism (ko00500), RNA transport (ko03013), purine metabolism (ko00230), Spliceosome (ko03040) and plant hormone signal transduction pathway (ko04075).

Identify differentially expressed genes using digital gene expression tag (DGE)
To identify differentially expressed genes in response to ROS, 6 DGE libraries have been generated from ROStreated rudimentary leaves at 0 h, 5 h and 10 h post treatment, with two biological replicates. The libraries produced over 1.41 G of 49 nt single-end read data with a Q20 percentage of about 98%. The percentage of unassigned base "N" is 0.01% and the average GC content is around 45% ( Table 2). To assay the normality of the RNA-Seq data in the 6 DGE libraries, we calculated the distribution of unique reads in each DGE libraries. This value is the ratio of the number of bases in a gene covered by unique mapping reads to the total bases in that from our transcriptome reference database. The distribution over different reads abundance categories showed similar patterns among all six libraries. Above 36% of the sequences have coverage more than 80% (Additional file 4). Next we calculated the unigene expression using the uniquely mapped DGE reads and normalized the results to RPKM. Results from the two biological replicates are highly similar, suggesting good reproducibility of the method (Additional file 5). We performed a pairwise comparison using 0 h as the control, and 5 h, or 10 h as the treatments. We also identified differentially expressed unigenes with FDR (false discovery rate) ≤ 0.001 and absolute value of log 2 Ratio ≥ 1. As a result, 5,865 unigenes were referred to as differentially expressed genes (DEGs) and used for the subsequent analysis (Additional file 6).
DEGs in the signal transduction pathways in response to ROS Table 5 shows the number of the DEGs involved in the plant hormone signal transduction pathway during 0 to 10 h of ROS-treatment. Numbers of those DEGs in the 4 significantly different expression patterns were calculated. A total of 57 DEGs were annotated in the plant hormone signal transduction pathways, including auxin, cytokinine, gibberellin, abscisic acid, ethylene, brassinosteroid and jasmonic acid.
In the auxin signal transduction pathway, 4 unigenes encoding auxin influx transport protein (AUX1) were found to be differentially expressed, among which two were annotated to profile 6, and one to profile 7 showing different pattern of up-regulation. It was found that 7 DEGs encode gretchenhagen-3 (GH3) or GH3 family protein. Four of them were clustered to profile 1, and 1 to profile 0 showing down-regulated trends. Only 1 DEG belonging to profile 1 was annotated to auxin-induced protein AUX/IAA. In the auxin signal transduction pathway, 7 out of the 15 DEGs showed down-regulated trends, and 3 of those showed up-regulated trends, indicating that DEGs of the down-regulated trends in this pathway were more than those of up-regulated trends. Similar results were found in the cytokinine, gibberellin, brassinosteroid and jasmonic acid signal transduction pathway. In the cytokinine signaling pathway, 2 unigenes encoding cytokinin receptor 1B or 1 (CRE1) and 1 encoding type-a response regulator (A-ARR) showed down-regulated trends, while the Unigene0026793 encoding histidine phosphotransfer protein (AHP) showed an up-regulated trend. In the gibberellin signal transduction pathway, 4 unigenes encoding gibberellic acid receptor or DELLA protein were found to be differentially expressed, where 1 was identified as a down-regulated profile, and the other as up-regulated  profiles. In the brassinosteroid signaling pathway, 6 out of the 8 DEGs showed down-regulated trends, while 1 showed as an up-regulated trend (Table 5, Figure 4). In the abscisic acid signal transduction pathway, 11 out of the 12 DEGs were clustered to profile 6 or 7 showing up-regulated trends. They encode abscisic acid receptors PYR (PYRabactin resistance)/PYL (PYR1-like), type 2C protein phosphatase (PP2C), SNF1 related protein kinase 2 (SnRK2), and ABA responsive element binding factor (ABF). Only Unigene0076273 encoding abscisic acid receptor showed a down-regulated trend. In the ethylene signal pathway, 11 out of 13 DEGs were also clustered to up-regulated profiles, including 3 unigenes encoding ethylene receptor (ETR), 2 encoding ein3binding F-box protein (EBF1/2), 3 encoding ethyleneinsensitive 3b (EIN3) and 2 encoding ERF transcription factor ERF1/2. Only 2 unigenes (Unigene0022597 and Unigene0022598) encoding mitogen activated protein kinase (MPK6) showed down-regulated trends (Table 5, Figure 4). These results showed that DEGs of the up-regulated trend in the abscisic acid and ethylene signaling pathway were much more than those of down-regulated trend, suggesting that most genes involved in these hormone signal transduction pathway were induced by ROS.

Confirm unigenes expression using real-time quantitative reverse transcription PCR
To confirm the accuracy and reproducibility of the transcriptome analysis results, 7 unigenes were selected for real-time quantitative reverse transcription PCR (qRT-PCR) validation. RNA samples from the ROStreated rudimentary leaves were used as templates. Primers of the candidate unigenes are shown in Additional file 8. The expression profiles of the candidate unigenes revealed by qRT-PCR data were consistent with those derived from sequencing ( Figure 5). Linear regression analysis of the fold change of the gene expression ratios between RNA-seq and qRT-PCR showed significantly positive correlation (Figure 6), confirming our transcriptome analysis.

Expression levels of the candidate genes in response to ROS
Based on the KEGG pathway enrichment analysis, it was found that most of the unigenes encoding the ABA and ethylene signal transduction components were induced by ROS. We then selected 10 candidate unigenes encoding these components and analyzed their expression levels during 0 to 10 h of ROS treatment. Primers of the candidate unigenes are shown in Additional file 8. The ABA and ethylene signal transduction components are shown in Additional file 9. Analysis of the gene expression revealed by qRT-PCR shows that genes encoding the ABA signal transduction components, including the abscisic acid receptor PYR1, protein phosphatase 2C, ABA responsive element-binding protein 2 and serine/ threonine-protein kinase SRK2E-like isoform 1 increased while that of the PYL5-like decreased during 0 to 10 h of treatment. Genes encoding the ethylene signal transduction components, including the ethylene response 3, Ein3-binding f-box protein 3, ethylene-insensive 3b and ERF transcription factor 4 increased while that of the mitogen-activated protein kinase decreased by the treatment (Figure 7).

Discussion
Stresses induce flowering in evergreen woody fruit trees [28][29][30]. ROS act as stress signals and are recognized as important signaling components in a wide range of processes. They are generated in chloroplast and peroxisomes in the light, in mitochondria in the dark and in non-green tissues [31,32]. These signals are perceived specifically by diverse mechanisms, such as the direct redox modification of transcription factors and other proteins, resulting in the regulation of transcriptome [33]. Our previous study showed that ROS, induced by MV, promoted continuative development of panicle primordia and abortion of rudimentary leaves in leafy panicle in litchi [14]. To understand the molecular mechanisms underlying flowering in litchi, identifying ROS responsive genes in rudimentary leaves is needed as well as those in panicle primordia. We had established an SSH library and had identified 93 ROS responsive genes in the panicle primordia [21]. In the present study, we sought to identify the stress-responsive genes in the litchi rudimentary leaves. To overcome the lack of a reference litchi genome, we used RNA-Seq data to de novo assemble a reference transcriptome. In total, our assembly contains 82,036 unigenes with a mean size of 710 bp. 47,596 unigenes were annotated to public protein databases. Using the transcriptome as a reference, we performed DESeq and identified 5,865 differentially expressed genes between un-treated (0 h) and ROS-treated (5 h or 10 h) rudimentary leaves. 2,052 unigenes showed up-regulated trends and 3,035 showed down-regulated trends from 0 to 10 h of treatments. Compared to the 93 ROS responsive genes identified by previous SSH experiment [21], RNA-Seq has identified significantly more DEGs in the rudimentary leaf libraries.
Plant hormones are signal molecules produced within the plant, and occur in extremely low concentrations, but regulate a wide range of processes, including determining the formation of flowers, stems, leaves, the shedding of leaves, the development and ripening of fruit, and in response to biotic and abiotic stresses. The plant hormone signals are perceived and transmitted to the nuclear by series signal transduction components to induce gene expression, resulting in a series of physiological processes. Our KEGG pathway enrichment analysis of the DEGs indicated that unigenes encoding the hormone signaling components were significantly enriched in the differentially expressed groups after MV treatment. These hormones included auxin, cytokinine, gibberellin, abscisic acid, ethylene, brassinosteroid, and jasmonic acid, suggesting that their signal components are responsive to ROS. It is believed that ROS signaling and redox balance is integrated with salicylic acid (SA) signaling [33]. SA-signaling pathway has been proved to have a role in controlling gene expression during senescence [34]. Though the SA-signaling pathway was not found to be significantly enriched in the ROStreated rudimentary leaves, some unigenes, such as NPR1 and PAD4 encoding SA-signaling components were found to be differentially expressed (Additional file 6, Unigene0034828 and Unigene0014177), suggesting that SA might also be involved in the ROS-induced rudimentary leaf abortion. Abscisic acid (ABA) is an essential hormone to control plant growth, development and adaptation to environmental stresses [35]. We found that 11 out of the 12 DEGs encoding abscisic acid signal components were up-regulated. The unigenes encoding PYR /PYL, PP2C, SnRK2, and ABF were induced by MV-driven ROS. Our gene expression levels of the components determined by qRT-PCR were consistent with those by RNA-seq, further confirming that the ABA signal transduction components were ROS responsive. ABA is essential for abscission and senescence of aged organs. It is involved in shading-induced abscission of apple fruits [36], and ethylene-associated abscission activation in citrus fruitlets [37]. In the present study, MV treatment induced downward growing of the rudimentary leaves, an early sign of abortion of the rudimentary leaves [4]. The crosstalk between ABA and ROS signaling is well known [38,39], and ROS is also involved in the ABA-enhanced LcAP1 expression in litchi [40]. Therefore, we hypothesize that the increase in the gene expression of ABA signaling components might play a role in the ROS-induced abortion of litchi rudimentary leaves.
We have also identified 13 differentially expressed genes encoding ethylene signaling components, and 11 of them were up-regulated after MV-treatment. Gene expression levels determined by qRT-PCR indicated that the unigenes encoding ethylene response 3, Ein3-binding f-box protein 3, ethylene-insensive 3b and ERF transcription factor 4 increased while that of the mitogen-activated protein kinase decreased by the ROS treatment, consisting with those revealed by RNA-seq. Our present study suggested that the ethylene signaling transduction components were ROS responsive. Ethylene as a gaseous hormone is involved in a variety of plant developmental adaptations including seed germination, organ senescence, fruit ripening, abscission and stress responses [41]. In agricultural practice, growers often use ethephon as ethylene producer to control rudimentary leaf growth in panicles and promote continual panicle development. We found that ethylene could increase H 2 O 2 levels in the rudimentary leaves [5]. We also showed that the petioles of rudimentary leaf displayed downward growing after MV-ROS treatment (Figure 1), similarly to those of the ethylene-treated leaves, which is a phenomenon of epinasty, suggesting that the MV-ROS could function through ethylene to induce abscission of litchi rudimentary leaves.
In Arabidopsis thaliana, MV was used to study the influence of chloroplastic ROS generation at the transcriptional level [42]. We compared our DEGs in response to MV with those identified by Scarpeci et al. [42]. To our surprise, 84% (267/316) of their differentially expressed   homology genes were found in our GEGs data set, including genes encoding some plant hormone signaling transduction components. Our study provided more information on the influence of MV generated-ROS on plant transcriptome.

Conclusions
In summary, we reported a comprehensive litchi leaf EST dataset generated by de novo assembly of next generation sequencing data. It is a valuable resource for future litchi genomic studies and will also benefit researches in other closely related species with significant agricultural importance. The differentially expressed genes dataset will also provide useful candidate genes for functional analysis of litchi flowering underlying the control of rudimentary leaf.

Plant material and experiment procedures
Commercially cultivated thirteen-year-old litchi trees cv. Nuomici grafted onto 'Huaizhi' were selected. The trees were grown at the experimental orchard of South China Agricultural University. About 6 cm length of branches with new flushes were cut off from the trees and immediately placed in water. Two hours later, the cuttings were treated with water or solutions containing 40 μM MV (Sigma) according to the method of Zhou et al. [14]. All the cuttings were placed in a growth chamber at 160 μmol/m 2 s 1 photosynthetic photon flux density at 20°C. The third and the fourth of the 0 h, 5 h and 10 h MV-treated rudimentary leaves as shown in Figure 1A were sampled, frozen in liquid nitrogen and stored at −80°C for RNA-seq and qRT-PCR. The proximal angle α, and the distal angle β of the third rudimentary leaves as shown in Figure 1A were measured. The angular dimension was measured by degree.

De novo assembly and annotation
For the assembly library, raw reads were filtered to remove those containing adapter and reads with more than 5% unknown nucleotides. Low quality reads were also remove, in which the percentage of low Q-value (≤10) base was more than 20%. Clean reads were de novo assembled by the Trinity Program [43]. To annotate the Figure 6 Coefficient analysis of fold change data between qRT-PCR and RNA-seq. Seven unigenes were selected for qRT-PCR. Data indicating relative transcript level from qRT-PCR are means of three replicates, and RPKMs from RNA-seq are means of two replicates. Scatterplots were generated by the log 2 expression ratios from RNA-seq (x-axis) and qRT-PCR (y-axis  dataset is available from the NCBI Short Read Archive (SRA) with an accession number SRA158542 (http:// www.ncbi.nlm.nih.gov/sra).

Expression annotation
Raw sequence data of the libraries for DGE profiling analyses were filtered to remove those containing adapter and reads with more than 10% unknown nucleotides, and reads with more than 50% of low quality base (value ≤5). Clean reads were mapped into the transcriptome reference database using SOAPaligner/ soap2 software. Not more than 2 mismatch bases were permitted, and unique mapped reads were obtained. The number of unique-match reads was calculated and normalized to RPKM (reads per kb per million reads) for gene expression analysis.
Comparison of unigene expression between treatments was according to DESeq as described by Abders and Huber [44]. P-value corresponds to differential gene expression test. FDR is a method to determine the threshold of P-value. In this experiment, the DEGs between 0 h and 5 h, or between 0 h and 10 h of ROS treatment were restricted with FDR ≤ 0.001 and the absolute value of log 2 Ratio ≥ 1.

Real-time PCR confirmation of the RNA-Seq data
First-strand cDNA was generated from 1 μg total RNA isolated from the rudimentary leaves using the superscript first-strand synthesis system (Invitrogen, USA). Primers for quantitative reverse transcription PCR (qRT-PCR) were designed using Primer Premier 5.0 software (Premier, Canada) and synthesized by Sangon Biotech (Shanghai) Co., Ltd. The litchi homologue Actin (GenBank accession number:HQ588865.1) was selected as reference. All the primers are shown in Additional file 8. qRT-PCR was performed on a Bio-Rad iQ5 Optical System Real Time PCR System (Bio-Rad, USA) using a SYBR Green based PCR assay. Each reaction mixture was 20 μL containing 6 μL of diluted first-strand cDNAs and 250 nM of each primer, SYBR Green PCR Master Mix (TaKaRa, Japan) 10 μL. The qPCRs were run as follows: 50°C for 2 min, 95°C for 10 min, followed by 40 cycles of 95°C for 30 s, 56°C for 30 s, and 72°C for 30 s in 96-well optical reaction plates (Bio-rad, USA). Each qRT-PCR analysis was performed in triplicate. Expression levels of the tested reference genes were determined by CT values and calculated by 2 -△△Ct .