Skip to main content


Springer Nature is making SARS-CoV-2 and COVID-19 research free. View research | View latest news | Sign up for updates

Molecular mechanism of ethylene stimulation of latex yield in rubber tree (Hevea brasiliensis) revealed by de novo sequencing and transcriptome analysis



Rubber tree (Hevea brasiliensis) is an important industrial crop cultivated in tropical areas for natural rubber production. Treatment of the bark of rubber trees with ehephon (an ethylene releaser) has been a routine measure to increase latex yield, but the molecular mechanism behind the stimulation of rubber production by ethylene still remains a puzzle. Deciphering the enigma is of great importance for improvement of rubber tree for high yield.


De novo sequencing and assembly of the bark transciptomes of Hevea brasiliensis induced with ethephon for 8 h (E8) and 24 h (E24) were performed. 51,965,770, 52,303,714 and 53,177,976 high-quality clean reads from E8, E24 and C (control) samples were assembled into 81,335, 80,048 and 80,800 unigenes respectively, with a total of 84,425 unigenes and an average length of 1,101 bp generated. 10,216 and 9,374 differentially expressed genes (DEGs) in E8 and E24 compared with C were respectively detected. The expression of several enzymes in crucial points of regulation in glycolysis were up-regulated and DEGs were not significantly enriched in isopentenyl diphosphate (IPP) biosynthesis pathway. In addition, up-regulated genes of great regulatory importance in carbon fixation (Calvin cycle) were identified.


The rapid acceleration of glycolytic pathway supplying precursors for the biosynthesis of IPP and natural rubber, instead of rubber biosynthesis per se, may be responsible for ethylene stimulation of latex yield in rubber tree. The elevated rate of flux throughout the Calvin cycle may account for some durability of ethylene-induced stimulation. Our finding lays the foundations for molecular diagnostic and genetic engineering for high-yielding improvement of rubber tree.


Rubber tree (Hevea brasiliensis (Willd. ex A. Juss.) Müll. Arg.), cultivated in many countries of tropical area, is the primary commercial source of natural rubber. Natural rubber (cis-1,4-polyisoprene) produced from rubber tree has many outstanding performance properties and is an important industrial material which cannot be reproduced by synthetic alternatives [1, 2].

Natural rubber is synthesized through the mevalonic acid (MVA) pathway with isopentenyl pyrophosphate (IPP) as the substrate and building block [15]. But there were evidences that the possibility that 1-dexoxy-D-xylulose 5-phosphate/2-C-methyl-D-erythritol 4-phosphate (DEX/MEP) pathway is involved in rubber biosynthesis cannot be excluded [69]. Whatever the isoprenoid biosynthesis pathway, sucrose is the only precursor of natural rubber [8, 10].

Natural rubber biosynthesis takes place in latex and the latex is the cytoplasm of laticifer cells (the highly specialized cells in phloem). The latex harvesting is conducted by tapping or cutting the bark of rubber tree at regular intervals ranging from 2 to 5 days. Since ethephon (chloro-2-ethylphosphonic acid, an ethylene generator or releaser) was known to enhance rubber production, stimulation with ethephon has become a common practice in worldwide rubber tree plantations [11, 12]. Usually, a 1.5- to 2-fold increase of latex yield can be achieved by the treatment with ethephon [13, 14].

Much effort has been devoted to the elucidation of the mechanism of ethylene action. Stimulation with ethephon was shown to prolong the flow of latex after tapping [13, 15]. Two ethylene-responsive aquaporin genes (HbPIP2;1 and HbTIP1;1) had been characterized and their functions possibly favor the prolongation of latex flow through regulation of water exchanges between inner liber and latex cells [16].

Physiological and biochemical studies showed that ethylene activated the general metabolism for latex regeneration between tappings, with adenylic pool, polysomes and rRNA contents accumulated in laticifers, and the activities of glutamine synthetase (GS) and chitinase up-regulated [13, 14, 17, 18]. Specifically, Tupy [1921] showed that bark application of auxins and ethylene enhanced the invertase activity and sucrose utilization in latex by increase of latex pH. Cytosolic alkalinization might be explained by the early activation of H+-translocating ATPase by ethylene [22].

In view of the fact that the laticifer is a strong sink for sucrose and sucrose importation into laticfer and quebrachitol absorption may be important for sustained sucrose demands for latex regeneration, Dusotoit-Coucaud et al. [2325] and Tang et al. [26] characterized several sucrose transporters and a polyol transporter, and suggested that they might be involved in ethylene-induced stimulation of latex production.

Liu et al. [27] constructed and screened two ethephon-induced latex SSH cDNA libraries, and found that the cDNAs associated with sucrose metabolism, regulation of coagulation, stability of lutoids and signal transduction were up-regulated and might be related to the ethephon action.

Due to the complexity of ethylene stimulation and the limitations of the methodologies, many researches have been undertaken on understanding the mechanism of ethylene action but it has still not resolved [28]. In contrast with the conventional methods such as single gene cloning and DNA microarrays which yield a limited amount of genetic information, RNA sequencing is powerful tool for analyzing differential gene expression with high resolution on the whole genome level [29, 30]. Particularly, transcriptome analysis can be employed to reveal relationships between plant gene expression and phenotypes [3133]. Transcriptome sequencing technology has been applied to investigate the biology of rubber tree for generation of tissue-specific transcriptomal data [34] and genome sequence [35], development of molecular markers [3639], identification of novel microRNAs [40], specific genes and gene families [4143].

In the present study, the first RNA sequencing project for deciphering the molecular mechanism of ethylene stimulation of latex yield in rubber tree was performed. Our objective was to identify the relevant metabolic pathways or major ethylene-responsive genes. It should be noted that, we treated the rubber tree bark with ethephon for 8 and 24 h, considering that it should be an effect of relatively long duration and the response to the stimulation should at last be transmitted to certain metabolic pathways.

Results and discussion

Sequencing and de novo assembly

Three cDNA libraries from bark tissue, C (control), E8 (ethephon treatment for 8 h) and E24 (ethephon treatment for 24 h), were sequenced by Illumina deep-sequencing and a total of approximately 55, 54 and 54 million raw reads for C, E8 and E24 were generated, respectively. After removal of low-quality reads, adaptor sequences and ambiguous reads, about 53, 51 and 52 million high-quality clean reads for C, E8 and E24 were obtained, respectively, with the Q20 (percentage of sequences with sequencing error rate lower than 1 %) over 98 % for the three samples (Table 1). Assembly of all trimmed reads produced 138,182-140,407 bp contigs from the three libraries with the average length exceeding 350 bp. The contigs were joined into unigenes based on the paired-end information, generating 80,800 (C), 81,335 (E8) and 80,048 (E20) unigenes, with an average length of 1,101 bp and N50 of 1,875 bp (50 % of the assembled bases were incorporated into unigenes of 1875 nt or longer) (Table 1). All unigenes were longer than 200 bp, the length of 29,541 (34.99 %) unigenes ranged from 201 to 400 bp, and unigenes with length longer than 2000 bp accounted for 16.78 % (14,165) of total unigenes. The length distribution of the unigenes is shown in Fig. 1.

Table 1 Overview of the sequencing and assembly
Fig. 1

The length distribution of the unigenes. The length distribution of the unigenes of C sample (a), E8 sample (b), E24 sample (c) and the all-unigenes (d)

Functional annotation and classification

Function annotation of the generated unigenes was performed by searching the reference sequences using BLASTX against NT (non-redundant NCBI nucleotide database), NCBI’s non-redundant protein databases (NR), SwissProt, Kyoto Encyclopedia of Genes and Genomes (KEGG) and COG (Cluster of Orthologous Groups). A total of 59,452 significant BLAST hits (70.42 % of all unigenes) were returned. Among them, 55,936 (66.26 %), 54,741 (64.84 %), 35,557 (42.12 %), 33,721 (39.94 %), 23,044 (27.30 %), 45,515 (53.91 %) unigenes were found in NR, NT, Swiss-Prot, KEGG, COG and GO database, respectively. Summary of functional annotations of unigenes by aligning to the Nr, Nt, Swiss-Prot, COG, GO and KEGG databases is shown in Additional file 1.

COG classification revealed that 44,030 out of all the assembled unigenes were clustered into 25 functional categories (Fig. 2). The largest category was “General function prediction only” (8185, 18.59 %), followed by “Transcription” (4143, 9.41 %), “Replication, recombination and repair” (3787, 8.60 %), “Signal transduction mechanisms” (3319, 7.54 %) and “Posttranslational modification, protein turnover, chaperones” (3102, 7.05 %). “Nuclear structure” (9, 0.02 %) was the smallest group and 2202 (5.00 %) unigenes were annoted as “Function unknown”.

Fig. 2

COG classification of the total assembled unigenes

GO assignments were used to classify the predicted functions of the unigenes and total of 367,851 unigenes were classified into three main categories: biological process (178,004, 48.39 %), cellular components (135,748, 36.90 %) and molecular function (54,099, 14.71 %) (Fig. 3). In the “biological process” category, the top six largest categories were “cellular process” (28,776, 7.82 %), “metabolic process” (27,921, 7.59 %), “single-organism process” (19,263, 5.24 %), “response to stimulus” (14,053, 3.82 %), “biological regulation” (11,538, 3.14 %) and “regulation of biological process” (10,666, 2.90 %). As for the molecular function, unigenes with binding (23,778, 6.46 %), catalytic activity (22, 256, 6.05 %) formed the largest groups.

Fig. 3

GO classification of the total assembled unigenes. The left y-axis indicates the percentage of genes for each functions; the right-axis indicates the correspondent number of genes

Analysis of KEGG pathways and differentially expressed genes

KEGG pathway-based analysis was conducted to obtain a better understanding of the biological functions of the unigenes. The results showed that 33,721 annotated transcripts were mapped to 128 KEGG pathways. Of the 33,721 unigenes, 7470 (22.15 %), 3517 (10.43 %), 2289 (6.79 %) and 1717 (5.09 %) were involved in the metabolic pathways, the biosynthesis of secondary metabolites, the plant-pathogen interaction and the plant hormone signal transduction, respectively.

FPKM (Fragments Per Kb per Million reads) method was used to calculate the expression levels of the unigenes to identify differentially expressed genes (DEGs). A total of 5,326 up-regulated unigenes and 4,890 down-regulated unigenes were identified in E8 compared with C, and 4,440 and 4,934 unigenes were up-regulated and down-regulated in E24 compared with C. DEGs between C and E-8, and C and E-24 are shown in Additional file 2. The top 20 most up-regulated and down regulated genes between C and E-8, and C and E-24 are shown in Additional file 3.

By performing the KEGG pathway enrichment analysis, 23 and 17 significantly enriched metabolic pathways or signal transduction pathways of DEGs in E8 and E24 were identified, respectively (Additional file 4). To examine the molecular basis of the ethylene stimulation of latex yield in rubber tree, we concentrated on glycolysis, IPP biosynthesis pathway, C3 carbon fixation (Calvin cycle), plant hormone signal transduction and TFs in relation to rubber biosynthesis.

DEGs involved in glycolysis and IPP pathway

Glycolysis pathway is shown in Fig. 4. In the glycolysis pathway, the genes of PFK (in step 3), FBA (in step 4), PGM (in step 8) and PK (in step 10) were significantly up-regulated in both E8 and E24 compared to C, and the gene expression of GAPDH (in step 6) was up-regulated in E24 compared to C (Fig. 4, Additional file 5). It is well known that the PFK reaction is the first irreversible step committed to glycolysis and PK is responsible for producing pyruvate as the end product of the pathway. Both PK and PFK reactions are crucial points in regulating the rate of glycolysis and in other words these steps are rate-determining steps [44, 45]. In the “bottom up” regulatory mode of plant glycolysis, primary and secondary regulation are exerted at the levels of PEP and fructose 6-phosphate, respectively [44, 45]. In addition, transgenic study demonstrated that PGM also plays an important role in the control of glycolysis [46, 47].

Fig. 4

Differential expression of unigenes involved in glycolysis in E8 and E24 compared to C samples of Hevea brasiliensis. Glycolysis comprises 10 step reactions. In step 1, the enzyme hexokinase (HK) phosphorylates glucose by transferring a phosphate group from ATP to glucose forming glucose 6-phosphate. In step 2, glucose 6-phosphate is converted to its isomer fructose 6-phosphate by phosphoglucoisomerase (PGI). In step 3, the enzyme phosphofructokinase (PFK) catalyzes the conversion of fructose 6-phosphate into fructose 1,6-bisphosphate using another ATP to transfer a phosphate group to fructose 6-phosphate. In step 4, Fructose 1,6-bisphosphate aldolase (FBA) splits fructose 1,6-bisphosphate into dihydroxyacetone phosphate (DHAP, also glycerone phosphate) and glyceraldehyde 3-phosphate (G3P). In step 5, triose-phosphate isomerase (TPI or TIM) catalyzes the reversible interconversion of DHAP and G3P. In step 6, G3P is converted to 1,3-bisphosphoglycerate (1,3-BPG) by glyceraldehyde 3-phosphate dehydrogenase (GAPDH). In step 7, phosphoglycerate kinase (PGK) catalyzes the reversible transfer of the phosphate group from 1,3-BPG to ADP generating 3-phosphoglycerate (3-PG) and ATP. In step 8, the conversion of 3-PG to 2-phosphoglycerate (2-PG) is catalyzed by phosphoglycerate mutase (PGM). In step 9, enolase (or phosphopyruvate hydratase) catalyzes the dehydration of 2-PG to form phosphoenolpyruvate (PEP). In step 10, pyruvate kinase (PK) catalyzes the transfer of the phosphate group from PEP to ADP, yielding pyruvate and ATP. In an acetyl-CoA metabolic pathway, acetyl-CoA synthetase (ACS) (or acetate : CoA ligase) catalyzes the interconversion between acetyl-CoA and acetate (Step 1′), aldehyde dehydrogenase (ALDH) reversibly catalyze the conversion of acetate into acetaldehyde (Step 2′) and alcohol dehydrogenase (ADH) facilitates the interconversion between aldehydes and alcohols (Step 3′). Cells with gray border lines in the upper rows represent differentially expressed unigenes in E8 compared to C and cells with green border lines in the lower rows represent differentially expressed unigenes in E24 compared to C. Relative levels of expression are showed by a color gradient from low (blue) to high (red)

Glycolysis pathway not only plays a crucial in energy generation for latex generation but also provides carbon building blocks for the biosynthesis of rubber and other organic constituents of latex [20, 44, 45]. Particularly, the breakdown of glucose by glycolysis produces the acetyl coenzyme A (acetyl-CoA) as the direct precursor in the MVP pathway, as well as the G3P and pyruvate which are the precursors in DEX/MEP pathway [48]. Up-regulation of the key enzyme genes described above possibly leads to an elevated rate of flux throughout the entire glycolytic pathway. Thus, the ethephon treatment enhances the latex yield through the rapid acceleration of the glycolysis for replenishing the precursors consumed in the biosynthesis of IPP and natural rubber molecules.

Moreover, we found that, in the metabolic pathway of acetyl-CoA into alcohol, the ACS (in step 1′) and ALDH (in step 2′) were down-regulated in both E8 and E24 compared to C (Fig. 4, Additional file 5). The down-regulation of the ACS and ALDH leads to the inhibition of the acetyl-CoA metabolic pathway under ethephon treatment, which may favor for the acetyl-CoA flowing toward the rubber synthesis.

In contrast, the IPP pathway was not found to be significantly enriched in this study and specifically, 3-hydroxy-3-metylglutaryl coenzyme A reductase (HMGR), which plays only limiting role in plant isoprenoid biosynthesis of MVA pathway [49], was slightly down-regulated (Additional files 6 and 7). This implied that carbohydrate utilization by glycolysis and availability of the precursors for IPP biosynthesis may critically determine the rate of rubber production, and the IPP pathway accounts little for the effect of ethephon stimulation. Although biosynthesis of natural rubber begins with the IPP pathway but the result that ethylene has little direct effect on accelerating rubber biosynthesis has been reported [28]. Interestingly, in one report about the mechanisms underlying super productivity of rubber tree, many rubber-biosynthesis-pathway genes showed no differential expressions between the super-high-yielding clone and the control [50].

DEGs involved in carbon fixation

The C3 carbon fixation (Calvin cycle) is demonstrated in Fig. 5. The expressions of RubisCo(in step 1), GADPH(in step 3), FBA(in step 5), aldolase(in step 8) and SBP(in step 9) were shown to be up-regulated in response to ethylene in both E8 and E24 samples (Fig. 5, Additional file 8). In the Calvin cycle, extensive studies demonstrated that RubisCO [51, 52], SBP [5355] and aldolase [56, 57] dominate control of photosynthetic carbon fixation and they represent the rate-limiting enzymes in the Calvin cycle (reviewed by Raines [58] and Skitt et al. [59]). The up-regaulation of these enzymes in governing substrate flux through the Calvin cycle in response to the application of ethephon speeds up carbon fixation and further enhances the sustainable rubber productivity.

Fig. 5

Differential expression of unigenes involved in carbon fixation in E8 and E24 compared to C samples of Hevea brasiliensis. The C3 carbon fixation proceeds through 13 steps in three phases. In carboxylation phase, ribulose-1,5-bisphosphate carboxylase/oxygenase (RubisCo) adds carbon dioxide to ribulose-1,5-bisphosphate (RuBP) to generate 3-PG (Step 1). Reduction phase include two reactions: PGK catalyses the phosphorylation of 3-PG by ATP to form 1,3-BPG (or 3-phospho-D-glyceroyl phosphate) and ADP (Step 2); GADPH catalyses the reduction of 1,3-BPG by NADPH to produce G3P (Step 3). The regeneration phase of the cycle comprises 9 steps and is initiated by the enzyme TPI which catalyses the interconversion of G3P and DHAP (Step 4). Then fructose-1,6-bisphosphate aldolase (FBA) reversibly catalyses the aldol condensation of G3P and DHAP into fructose-1,6-bisphosphate (Step 5). The resulting fructose-1,6-bisphosphatase is converted into fructose 6-phosphate by fructose-1,6-bisphosphatase (FBP) (Step 6). Transketolase (TK) catalyzes the transfer of a 2-carbon fragment from fructose 6-phosphate to G3P, affording erythrose-4-phosphate and xylulose-5-phosphate (Step 7). Erythrose-4-phosphate and a DHAP are converted into sedoheptulose-1,7-bisphosphate by aldolase (Step 8). The cleavage of sedoheptulose-1,7-bisphosphate into sedoheptulose-7-phosphate and an inorganic phosphate ion is catalyzed by sedoheptulose-1,7-bisphosphatase (SBP) (Step 9). The reversible interconversion of sedoheptulose-7-phosphate and G3P into ribose 5-phosphate and xylulose 5-phosphate is catalyzed by TK (Step 10). Ribose-5-phosphate isomerase (RPI) interconverts ribose-5-phosphate and ribulose-5-phosphate (Step 11). Finally, phosphoribulokinase (PRK) phosphorylates ribulose-5-phosphate into RuBP (Step 12). Cells with gray border lines in the upper rows represent differentially expressed unigenes in E8 compared to C and cells with green border lines in the lower rows represent differentially expressed unigenes in E24 compared to C. Relative levels of expression are showed by a color gradient from low (blue) to high (red)

In plants, carbon-containing compounds including sucrose, IPP and natural rubber, ultimately comes from photosynthesis or from stored photosynthetic products. Wititsuwannakul [60] found that the HMGR activity in latex of rubber tree showed a diurnal variation pattern, and suggested that the regulation of the light-dark phase was probably due to physiological processes associated with photosynthesis and one or more of the three components of the rubber biosynthesis (acetyl-CoA, NADPH and ATP) become a limiting factor causing the decline in HMGR activity during the dark period. Therefore, to a certain extent, the photosynthetic carbon fixation may play an important role in ethylene-induced stimulation of latex production by constantly supplying carbohydrate for glycolytic transformation.

DEGs of hormone signaling components and transcription factors

The responses of plant metabolism, growth, development and defense to exogenous ethylene application can hardly do without the ethylene perception and signal transduction [61]. Generally, hormone responses are the output of multiple pathway interaction and crosstalk [6164]. Ethylene signaling may exert its functions via interaction with other major plant hormones such as brassinosteroid (BR), gibberellin (GA),auxin, cytokinin, salicylic acid (SA), jasmonate (JA), abscisic acid (ABA). The transcripts of signaling components such as ETHYLENE RESPONSE (ETR), ETHYLENE RESPONSE FACTOR 1 (ERF1), ETHYLENE INSENSITIVE 3 (EIN3) and EIN3 binding F-Box protein 1/2 (EBF1/2) for ethylene, BRASSINOSTEROID-INSENSITIVE 2 (BIN2), BRASSINAZOLE RESISTANT 1/2 (BZR1/2) and TCH4 for BR, DELLA for GA, INDOLE ACETIC ACID (IAA) for auxin, type-A response regulator (A-ARR) for cytokinin, TGA for SA, and MYC2 for JA were found to be remarkably accumulated after ethephon treatment (Additional files 9 and 10). Although all these hormones have been reported for being linked to growth regulation in certain manner [64], how the ethylene interacts and coordinates with other hormones in relation to the stimulation of latex production in rubber tree still remains to be elucidated.

Transcription factors (TFs) play important roles in the control of many of the biological processes in a cell or organism by the regulation of gene expression [65]. In addition, TFs also play crucial roles in the cross-talk between hormone signalling pathways [66]. A total of 1752 DEGs of putative TFs were identified in this study (Additional file 11) and the top five up-regulated and down-regulated TFs in E8 and E24 compared to C were listed in Additional file 12. Among them, the largest gene family was the The v-myb avian myeloblastosis viral oncogene homolog family (MYB) (88, 9.10 %), followed by MYB related (67, 7.44 %), the basic helix-loop-helix family (bHLH) (58, 6.44 %), C2H2 family (52, 5.77 %), and the ethylene-responsive element binding factor family (ERF) (44, 4.88 %). Although their functions relevant to ethylene responses for the production enhancement of rubber remained to be unknown but the data will be a valuable resource for the discovery of candidate genes related to the complex regulatory networks involved in the response.

The results obtained by qRT-PCR analysis matched the expression found by transcriptome analysis (Fig. 6). The consistent expression patterns further increased the confidence of our findings and the transcriptome dataset reported here which will be a valuable supplement to the rubber tree genomic and transcriptome information.

Fig. 6

qRT-PCR validation of three DEGs involved in glycolysis and carbon fixation of Hevea brasiliensis bark. The axis indicates treatments; the y-axis indicates relative expression level. a PFK; b Rubisco; c GADPH

Differential gene expression, especially the up-regulated expression of the key genes in the rubber-biosynthesis-related metabolic pathways in response to ethylene allowed us to find the rubber-yield limiting factors. Based on our finding, it is possibly to develop a molecular diagnostic for rubber yield and to improve the rubber tree for high yield through the genetic modification of those key genes [67].


In the present study, the molecular mechanism behind the ethylene stimulation of rubber yield was revealed. Up-regulation of the key enzymes in the glycolysis pathway and the C3 carbon fixation may be responsible for enhancing the rubber production by ethephon treatment. Identification of the rubber-yield limiting factors possibly leads to development of a molecular approach to assay and predict the productivity of rubber tree, and obtainment of new high yielding cultivars through genetic engineering. Moreover, our transcriptome data provides an useful resources for gene mining for high production of rubber.


Plant material and RNA extraction

Rubber trees of clone PR107 were planted at the Experimental Farm of the Chinese Academy of Tropical Agriculture Science in 2000, and opened for tapping in 1995, on the s/2 d/4 system (half spiral tapped every 4 days) and with 1.5 %-ethephon treatment 1 day before tapping at intervals of 15 days.

Trees of equal size in the same plot were selected and divided into three groups, each with six trees. Ethephon (2-chloroethane phosphonic acid; 1.5 % v/v) was applied using a brush on the tapping panel. Bark samples were collected from two groups by tapping at 8 and 24 h after ethephon application, respectively. Concurrently, bark was collected from a corresponding group of trees at 8 h after treatment in the same way with water only as control. The samples were immediately frozen in liquid nitrogen and shipped on dry ice to BGI Life Tech Co., Ltd (Shenzhen, China) for Illumina sequencing.

Bark RNA was isolated using the TRIzol® Reagent (Invitrogen) following the protocol in the manufacturer’s instructions. RNA integrity was confirmed by a 2100 Bioanalyzer (Agilent Technologies).

cDNA library construction and sequencing

Three cDNA libraries, C (control), E8 (Ethephon treatment for 8 h) and E24 (Ethephon treatment for 24 h), were generated using the mRNA-Seq 8 sample prep Kit (Illumina) according to the manufacturer’s instructions. Magnetic beads containing poly-T molecules was used to isolate the poly(A) mRNA from 20 μg of total RNA. Following purification, the samples were fragmented into small pieces using divalent cations at 94 °C for 5 min and converted into the first and second-strand cDNA with the SuperScript double-stranded cDNA synthesis kit (Invitrogen, CA). Then, the synthesized cDNA was subject to end repair and adenylation of 3′ ends and purified using the QIAquick PCR Purification Kit (QIAGEN). Afterward, Illumina paired end adapters were ligated to the resulting cDNA fragments. Each cDNA library was finally constructed with an insert size of 200 bp. After quality validating by an Agilent Technologies 2100 Bioanalyzer, deep sequencing was performed with Illumina HiSeqTM 2000 (Illumina Inc., San Diego, CA, USA).

De novo assembly and gene annotation

Raw reads were filtered by the Illumina pipeline before the assembly. The reads with more than 20 % of bases with a quality value ≤10, unknown nucleotides higher than 5 % and adaptor contamination were removed. Transcriptome de novo assembly was then conducted with short reads assembling program – Trinity (release 20130225; running parameter: seqType fq --min_contig_length 100) ( [68, 69]. The resulting sequences of trinity, termed Unigenes, from each sample’s assembly were taken into further process of sequence splicing and redundancy removing to acquire non-redundant Unigenes as long as possible. Then, the Unigenes were split into two classes: clusters (CL, with some Unigenes which similarity between them was higher than 70 %) and singletons (Unigene). At last, Blast X (v2.2.26 + x64-linux) alignment (an E-value <0.00001; running parameter: -F F -e 1e-5 -p blastx) between the Unigenes and protein databases like NR (release 20130408), SwissProt (release 201303), KEGG (release 63.0) and COG (release 20090331) was carried out. Sequence direction of and functional annotations to the Unigenes were decided and assigned with the best aligning results. The sequence direction of the Unigene unaligned to non of the above databases was determined by a software named ESTScan (v3.0.2) [69]. With the NR annotation, the Blast2GO program (v2.5.0; release 20120801) [70] was used to get the GO annotation for the Unigenes. The GO functional classification for the Unigenes was produced by the WEGO software [71] and the pathway assignments were performed with the help of KEGG database [72].

Differential gene expression analysis

The normalized expression levels of the Unigenes were calculated using the FPKM method [73]. Then, the identification of DEGs was conducted between C and E8, and C and E24 using a computational method [74] included in SOAPaligner/soap (; v2.21; running parameter: -m 0 -x 500 -s 40 -l 35 -v 5 -r 1), a tool of the Short Oligonucleotide Analysis Package (SOAP) for the RNA-Seq data analysis. The significance of differential transcript abundance was judged with the false discovery rate (FDR) value [75]. Only those DEGs with FDR ≤0.001 with the absolute fold change ≥2 were reserved. The enrichment analysis was performed to find GO terms in which DEGs are significantly enriched comparing to the whole transcriptome, using the hypergeometric test. For the pathway enrichment analysis, all DEGs were mapped to terms in the KEGG database to identify significantly over-represented metabolic pathways or signal transduction pathways. Then the hypergeometric test was performed for the statistical analysis, while the Bonferroni correction was adopted for the multiple testing correction, with a q-value cutoff of ≤0.05. From the results, we mainly focused on the differentially regulated pathways closely relevant to the biosynthesis of natural rubber. In order to evaluate expression levels of individual genes in these pathways, we re-mapped the Unigenes to the sequences of each gene and all FPKM values were added together. Then the final FPKM values were used as the expression levels of the individual gene [76]. For these genes, Log2 fold changes were calculated between C and E8, and between C and E24, respectively, while the results were presented and illustrated in the pathways.

Plant Transcription Factor Database (PlantTFDB, v3.0) is a comprehensive database of transcription factors (TFs) in plants, with detailed annotation and classification information [77]. The sequences of all existing TFs were retrieved from PlantTFDB. Then, Blast X search with PlantTFDB was performed and the TFs were identified with an E-value cutoff of ≤1E-5.

Validation of gene expression by qRT-PCR

To validate our transciptome results, expression of three randomly chosen key genes (PFK, Rubisco and GADPH) in glycolysis and carbon fixation with significant expression changes in the transcriptome data was verified by qRT-PCR. Total RNA was isolated from the equal amount of bark tissues of three rubber trees of each treatment by following the protocol described by Venkatachalam et al. [78]. The first strand cDNA was synthesized from 2.5 μg of total RNA through a RevertAidTM Premium first strand cDNA synthesis kit (Fermentas). The standard curve for each target gene was obtained by qRT-PCR with series cDNA dilutions of cDNA. The reaction mixture (20 μl) for qRT-PCR comprised of 10 μl SYBR Premix Ex TaqII, 6.8 μl EASY Dilution, 0.4 μl 10 μM Forward primer and 0.4 μl 10 μM Reverse primer. The PCR reactions were performed on an CFX96TM Real-Time PCR Detection System (Bio-Rad) as follows: 95 °C for 30 s, followed by 40 cycles of 95 °C for 5 s, and then annealing at 60–95 °C for 30 s. Expression analysis of each gene was confirmed in three independent reactions with Actin gene as an internal control for normalization of the expression levels of the chosen transcripts. The relative expression of the genes was calculated using the 2Ct method. The primer sequences used for qRT-PCR are listed in Table 2.

Table 2 The forward and reverse primers used in validation experiment of gene expression by qRT-PCR

Availability of data and materials

The datasets supporting the conclusions of this article are included within the article and its additional files. The clean reads for C, E8 and E24 have been deposited in the National Center for Biotechnology Information (NCBI) Sequence Read Archive (SRA) under accession numbers SRX1134637, SRX1134639 and SRX1134640, respectively. The RNA-seq data in this article has been deposited in the NCBI’s Gene Expression Omnibus (GEO) under accession number GSE78145.



type-A response regulator


abscisic acid


acetyl-coenzyme A


acetyl-coenzyme A synthetase or acetate, CoA ligase


alcohol dehydrogenase


aldehyde dehydrogenase


the basic helix-loop-helix family


brassinosteroid-insensitive 2






brassinazole resistant 1/2


Cluster of Orthologous Groups


differentially expressed genes


1-dexoxy-D-xylulose 5-phosphate/2-C-methyl-D-erythritol 4-phosphate


dihydroxyacetone phosphate or glycerone phosphate


EIN3 binding F-Box protein 1/2


ethylene insensitive 3


ethylene response factor 1


the ethylene-responsive element binding factor family


ethylene response


fructose 1,6-bisphosphate aldolase




false discovery rate


fragments per kb per million reads



G3P or GA3P:

glyceraldehyde 3-phosphate


glyceraldehyde 3-phosphate dehydrogenase


gene ontology


glutamine synthetase




3-hydroxy-3-metylglutaryl coenzyme A reductase


indole acetic acid


isopentenyl diphosphate




SwissProt, Kyoto Encyclopedia of Genes and Genomes


mevalonic acid


v-myb avian myeloblastosis viral oncogene homolog


National Center for Biotechnology Information


NCBI’s non-redundant protein databases


non-redundant NCBI nucleotide database






phosphoglycerate kinase


phosphoglycerate mutase








Plant Transcription Factor Database


pyruvate kinase




ribose-5-phosphate isomerase


ribulose-1,5-bisphosphate carboxylase/oxygenase




salicylic acid




short oligonucleotide analysis package


sequence read archive


transcription factors




triose-phosphate isomerase


  1. 1.

    Cornish K. Biochemistry of natural rubber, a vital raw material, emphasizing biosynthetic rate, molecular weight and compartmentalization, in evolutionarily divergent plant species. Nat Prod Rep. 2001;18:182–9.

  2. 2.

    Whalen M, McMahan C, Shintani D. Development of crops to produce industrially useful natural rubber. In: Bach TJ, Rohmer M, editors. Isoprenoid synthesis in plants and microorganisms: new concepts and experimental approaches. New York: Springer Science+Business Media; 2013. p. 329–45.

  3. 3.

    Hepper CM, Audley BG. The biosynthesis of rubber from hydroxyl-methyl-glutaryl coenzyme A in Hevea brasiliensis latex. Biochem J. 1969;114:379–86.

  4. 4.

    Skilleter DN, Kekwick RGO. The enzymes forming isopentenyl pyrophosphate from 5-phosphomevalonate (mevalonate-5-phosphate) in the latex of Hevea brasiliensis. Biochem J. 1971;124:407–15.

  5. 5.

    Kekwick RGO. The formation of isoprenoids in Hevea latex. In: d’Auzac J, Jacob JL, Chrestin L, editors. Physiology of rubber tree latex. Boca Raton: CRC Press; 1989. p. 145–64.

  6. 6.

    Ko JH, Chow KS, Han KH. Transcriptome analysis reveals novel features of the molecular events occurring in the laticifers of Hevea brasiliensis (para rubber tree). Plant Mol Biol. 2003;53(4):479–92.

  7. 7.

    Seetang-Nun Y, Sharkey TD, Suvachittanont W. Molecular cloning and characterization of two cDNAs encoding 1-deoxy-d-xylulose 5-phosphate reductoisomerase from Hevea brasiliensis. J Plant Physiol. 2008;165(9):991–1002.

  8. 8.

    Chow KS, Wan KL, Isa MNM, Bahari A, Tan SH, Harikrishna K, Yeang HY. Insights into rubber biosynthesis from transcriptome analysis of Hevea brasiliensis latex. J Exp Bot. 2007;58:2429–40.

  9. 9.

    Chow KS, Mat-Isa MN, Bahari A, Ghazali AK, Alias H, Mohd-Zainuddin Z, Hoh CC, Wan KL. Metabolic routes affecting rubber biosynthesis in Hevea brasiliensis latex. J Exp Bot. 2012;63:1863–71.

  10. 10.

    d’Auzac J. Mise en evidence de la glycolyse et de ses relations avec la biosynthêse du caoutchouc au sein du latex d’ Hevea brasiliensis. Rev Gén Caoutch Plast. 1964;41:1831–4.

  11. 11.

    Abraham PD, Wycherley PR, Pakianathan SW. Stimulation of latex flow in Hevea brasiliensis of 4-amino-3,5,6-trichloropicolinic acid and 2-chloroethane-phosphonic acid. J Rubber Res. 1968;20:291–305.

  12. 12.

    d’Auzac J, Ribailler D. L’éthylène, nouvel agent de la production de latex chez l’Heuea brasiliensis. CR Acad Sci Paris. 1969;268:3046–9.

  13. 13.

    Coupé M, Chrestin H. Physico-chemical and biochemical mechanisms of the hormonal (ethylene) stimulation: early biochemical events induced, in Hevea latex, by hormonal bark stimulation. In: d’Auzac J, Jacob JL, Chrestin H, editors. Physiology of rubber tree latex. Boca Raton: CRC Press Inc; 1989. p. 295–319.

  14. 14.

    Pujade-Renaud V, Clement A, Perrot-Rechenmann C, Prevôt J-C, Chrestin H, Jacob J-L, Guern J . Ethylene-induced increase in glutamine synthetase activity and mRNA levels in Hevea brasiliensis latex cells. Plant Physiol. 1994;105:127–32.

  15. 15.

    Audley BG, Archer BL, Carruthers IB. Metabolism of ethephon (2-chloroethylphosphonic acid) and related compounds in Hevea brasiliensis. Arch Environ Contam Toxicol. 1976;4:183–200.

  16. 16.

    Tungngoen K, Kongsawadworakul P, Viboonjun U, Katsuhara M, Brunel N, Sakr S, Narangajavana J, Chrestin H . Involvement of HbPIP2;1 and HbTIP1;1 aquaporins in ethylene stimulation of latex yield, through regulation of water exchanges between inner liber and latex cells in Hevea brasiliensis. Plant Physiol. 2009;151:843–56.

  17. 17.

    Amalou Z, Bangratz J, Chrestin H. Ethrel (ethylene releaser) induced increases in the adenylate pool and transtonoplast pH within latex cells. Plant Physiol. 1992;98:1270–6.

  18. 18.

    Tupy J. Ribosomal and polyadenylated RNA content of rubber tree latex, association with sucrose level and latex pH. Plant Sci. 1988;55(2):137–44.

  19. 19.

    Tupy J. The activity of the latex invertase and latex production of Heven brasiliensis. Physiol Veg. 1973;11:633–41.

  20. 20.

    Tupy J. Stimulatory effects of 2,4-dichlorophenoxyacetic acid and of 1-naphthylacetic acid on sucrose level, invertase activity and sucrose utilization in the latex of Hevea brasiliensis. Planta. 1969;88(2):144–53.

  21. 21.

    Tupy J. Control of carbohydrate metabolism by ethylene in latex vessels of Hevea brasiliensis Muel. Arg. In relation to rubber production. Biol Plant. 1976;18(5):373–84.

  22. 22.

    Gidrol X, Chrestin H, Mounoury G, D’Auzac J. Early activation by ethylene of the tonoplast H-pumping ATPase in the latex from Hevea brasiliensis. Plant Physiol. 1988;86:899–903.

  23. 23.

    Dusotoit-Coucaud A, Brunel N, Kongsawadworakul P, Viboonjun U, Lacointe A, Julien J-L, Chrestin H, Sakr S . Sucrose importation into laticifers of Hevea brasiliensis, in relation to ethylene stimulation of latex production. Ann Bot. 2009;104:635–47.

  24. 24.

    Dusotoit-Coucaud A, Kongsawadworakul P, Maurousset L, Viboonjun U, Brunel N, Pujade-Renaud V, Chrestin H, Sakr S. Ethylene stimulation of latex yield depends on the expression of a sucrose transporter (HbSUT1B) in rubber tree (Hevea brasiliensis). Tree Physiol. 2010;30(12):1586–98.

  25. 25.

    Dusotoit-Coucaud A, Porcheron B, Brunel N, Kongsawadworakul P, Franchel J, Viboonjun U, Chrestin H, Lemoine R, Sakr S. Cloning and characterization of a new polyol transporter (HbPLT2) in Hevea brasiliensis. Plant Cell Physiol. 2010;51(11):1878–88.

  26. 26.

    Tang C, Huang D, Yang J, Liu S, Sakr S, Li H, Zhou Y, Qin Y. The sucrose transporter HbSUT3 plays an active role in sucrose loading to laticifer and rubber productivity in exploited trees of Hevea brasiliensis (para rubber tree). Plant Cell Environ. 2010;33(10):1708–20.

  27. 27.

    Liu KC, Yang Y, Zhao LH, Zhang ZL. Primary analysis and construction of ethephon-induced latex SSH cDNA library from Hevea brasiliensis. Chin J Trop Crops. 2007;28:1–4.

  28. 28.

    Zhu J, Zhang Z. Ethylene stimulation of latex production in Hevea brasiliensis. Plant Signal Behav. 2009;4(11):1072–4.

  29. 29.

    Wang Z, Gerstein M, Snyder M. RNA-Seq: a revolutionary tool for transcriptomics. Nat Rev Genet. 2009;10:57–63.

  30. 30.

    Wilhelm BT, Landry JR. RNA-Seq—quantitative measurement of expression through massively parallel RNA-sequencing. Methods. 2009;48:249–57.

  31. 31.

    Metzker ML. Sequencing technologies-the next generation. Nat Rev Genet. 2010;11:31–46.

  32. 32.

    Ozsolak F, Milos PM. RNA sequencing: advances, challenges and opportunities. Nat Rev Genet. 2011;12:87–98.

  33. 33.

    Van Verk MC, Hickman R, Pieterse CMJ, Van Wees SCM. RNA-Seq: revelation of the messengers. Trends Plant Sci. 2013;18:175–9.

  34. 34.

    Xia Z, Xu H, Zhai J, Li D, Luo H, He C, Huang X.. RNA-Seq analysis and de novo transcriptome assembly of Hevea brasiliensis. Plant Mol Biol. 2011;77:299–308.

  35. 35.

    Rahman AY, Usharraj AO, Misra BB, Thottathil GP, Jayasekaran K, Feng Y, Hou S, Ong SY, Ng FL, Lee LS, Tan HS, Muhd Sakaff MKL, Teh BS, Khoo BF, Badai SS, Ab Aziz N,Yuryev A, Knudsen B, Dionne-Laporte A, Mchunu NP, Yu Q, Langston BJ, Freitas TAK, Young AG, Chen R, Wang L, Najimudin N, Saito JA, Alam M. Draft genome sequence of the rubber tree Hevea brasiliensis. BMC Genomics. 2013;14:75.

  36. 36.

    Triwitayakorn K, Chatkulkawin P, Kanjanawattanawong S, Sraphet S, Yoocha T, Sangsrakru D, Chanprasert J, Ngamphiw C, Jomchai N, Therawattanasuk K, Tangphatsornruang S. Transcriptome sequencing of Hevea brasiliensis for development of microsatellite markers and construction of a genetic linkage map. DNA Res. 2011;18:471–82.

  37. 37.

    Li D, Deng Z, Qin B, Liu X, Men Z. De novo assembly and characterization of bark transcriptome using Illumina sequencing and development of EST-SSR markers in rubber tree (Hevea brasiliensis Muell. Arg.). BMC Genomics. 2012;13:192.

  38. 38.

    Mantello CC, Cardoso-Silva CB, da Silva CC, de Souza LM, Scaloppi Junior EJ, de Souza Gonçalves P, Vicentini R; de Souza AP. De novo assembly and transcriptome analysis of the rubber tree (Hevea brasiliensis) and SNP markers development for rubber biosynthesis pathways. PLoS One. 2014;9(7):e102665.

  39. 39.

    Rippel Salgado L, Martins Koop D, Guariz Pinheiro D, Rivallan R, Le Guen V, Nicolás MF, De Almeida LGP, Ribeiro Rocha V, Magalhaes M, Gerber AL, Figueira A, De Mattos Cascardo JC, Ribeiro de Vasconcelos AT, Araújo Silva W, Coutinho LL, Garcia D. De novo transcriptome analysis of Hevea brasiliensis tissues by RNA-seq and screening for molecular markers. BMC Genomics. 2014;15:236.

  40. 40.

    Gébelin V, Argout X, Engchuan W, Pitollat B, Duan C, Montoro P, Leclercq J. Identification of novel microRNAs in Hevea brasiliensis and computational prediction of their targets. BMC Plant Biol. 2012;12:18.

  41. 41.

    Duan C, Argout X, Gébelin V, Summo M, Dufayard JF, Leclercq J, Hadi K, Piyatrakul P, Pirrello J, Rio M, Champion A, Montoro P. . Identification of the Hevea brasiliensis AP2/ERF superfamily by RNA sequencing. BMC Genomics. 2013;14:30.

  42. 42.

    Piyatrakul P, Yang M, Putranto R-A, Pirrello J, Dessailly F, Hu S, Summo M, Theeravatanasuk K, Leclercq J, Kuswanhadi, Montoro P.. Sequence and expression analyses of ethylene response factors highly expressed in latex cells from Hevea brasiliensis. PLoS ONE. 2014;9(6):e99367.

  43. 43.

    Liu J-P, Xia Z-Q, Tian X-Y, Li Y-J. Transcriptome sequencing and analysis of rubber tree (Hevea brasiliensis Muell.) to discover putative genes associated with tapping panel dryness (TPD). BMC Genomics. 2015;16:398.

  44. 44.

    Givan CV. Evolving concepts in plant glycolysis: two centuries of progress. Biol Rev. 1999;74(3):277–309.

  45. 45.

    Plaxton W. The organization and regulation of plant glycolysis. Annu Rev Plant Physiol Plant Mol Biol. 1996;47:185–214.

  46. 46.

    Westram A, Lloyd JR, Roessner U, Riesmeier JW, Kossmann J. Increases of 3-phosphoglyceric acid in potato plants through antisense reduction of cytoplasmic phosphoglycerate mutase impairs photosynthesis and growth, but does not increase starch contents. Plant Cell Environ. 2002;25:1133–43.

  47. 47.

    Fernie AR, Carrari F, Sweetlove LJ. Respiratory metabolism: glycolysis, the TCA cycle and mitochondrial electron transport. Curr Opin Plant Biol. 2004;7(3):254–61.

  48. 48.

    Lichtenthaler HK, Schwender J, Disch A, Rohmer M. Biosynthesis of isoprenoids in higher plant chloroplasts proceeds via a mevalonate-independent pathway. FEBS Lett. 1997;400(3):271–4.

  49. 49.

    Rodríguez-Concepción M, Campos N, Ferrer A, Boronat A. Biosynthesis of isoprenoid precursors in Arabidopsis. In: Bach TJ, Rohmer M, editors. Isoprenoid synthesis in plants and microorganisms: new concepts and experimental approaches. New York: Springer Science+Business Media; 2013. p. 439–56.

  50. 50.

    Tang C, Xiao X, Li H, Fan Y, Yang J, Qi J, Li H. Comparative analysis of latex transcriptome reveals putative molecular mechanisms underlying super productivity of Hevea brasiliensis. PLoS One. 2013;8(9):e75307.

  51. 51.

    Stitt M, Schulze D. Does Rubisco control the rate of photosynthesis and plant-growth – an exercise in molecular ecophysiology. Plant Cell Environ. 1994;17:465–87.

  52. 52.

    Mate CJ, Hudson GS, Von Caemmerer S, Evans JR, Andrews TJ. Reduction of ribulose-bisphophate carboxylase activase levels in tobacco (Nicotiana tabacum) by antisense RNA reduces ribulose-bisphosphate carboxylase carbamylation and impairs photosynthesis. Plant Physiol. 1993;102:1119–28.

  53. 53.

    Harrison EP, Olcer H, Lloyd JC, Long SP, Raines CA. Small decreases in SBPase cause a linear decline in the apparent RuBP regeneration rate, but do not affect Rubsico carboxylation capacity. J Exp Bot. 2001;52:1779–84.

  54. 54.

    Olcer H, Lloyd JC, Raines CA. Photosynthetic capacity is differentially affected by reductions in sedoheptulose-1,7-bisphosphatase activity during leaf development in transgenic tobacco plants. Plant Physiol. 2001;125:982–9.

  55. 55.

    Lefebvre S, Lawson T, Fryer M, Zakhleniuk OV, Lloyd JC, Raines CA. Increased sedoheptulose-1,7-bisphosphatase activity in transgenic tobacco plants stimulates photosynthesis and growth from an early stage in development. Plant Physiol. 2005;138:451–60.

  56. 56.

    Haake V, Zrenner R, Sonnewald U, Stitt M. A moderate decrease of plastid aldolase activity inhibits photosynthesis, alters the levels of sugars and starch, and inhibits growth of potato plants. Plant J. 1998;14:147–57.

  57. 57.

    Haake V, Geiger M, Walch-Liu P, Engels C, Zrenner R, Stitt M. Changes in aldolase activity in wild-type potato plants are important for acclimation to growth irradiance and carbon dioxide concentration, because plastid aldolase exerts control over the ambient rate of photosynthesis across a range of growth conditions. Plant J. 1999;17:479–89.

  58. 58.

    Stitt M, Lunn J, Usadel B. Arabidopsis and primary photosynthetic metabolism – more than the icing on the cake. Plant J. 2010;61(6):1067–91.

  59. 59.

    Raines CA. The Calvin cycle revisited. Photosynth Res. 2003;75:1–10.

  60. 60.

    Wititsuwannakul R. Diurnal variation of 3-hydroxy-3-methylglutaryl coenzyme A reductase activity in latex of Hevea brasiliensis and its relation to rubber content. Experientia. 1986;42:44–5.

  61. 61.

    Santner A, Calderon-Villalobos LI, Estelle M. Plant hormones are versatile chemical regulators of plant growth. Nat Chem Biol. 2009;5:301–7.

  62. 62.

    Santner A, Estelle M. Recent advances and emerging trends in plant hormone signalling. Nature. 2009;459:1071–8.

  63. 63.

    Wolters H, Jurgens G. Survival of the flexible: hormonal growth control and adaptation in plant development. Nat Rev Genet. 2009;10:305–17.

  64. 64.

    Depuydt S, Hardtke CS. Hormone signalling crosstalk in plant growth regulation. Curr Biol. 2011;21(9):R365–73.

  65. 65.

    Riechmann JL, Heard J, Martin G, Reuber L, Jiang C, Keddie J, Adam L, Pineda O, Ratcliffe OJ, Samaha RR, Creelman R, Pilgrim M, Broun P, Zhang JZ, Ghandehari D, Sherman BK, Yu G. Arabidopsis transcription factors: genome-wide comparative analysis among eukaryotes. Science. 2000;290:2105–10.

  66. 66.

    Jaillais Y, Chory J. Unraveling the paradoxes of plant hormone signaling integration. Nat Struct Mol Biol. 2010;17:642–64.

  67. 67.

    Chrestin H, Gidrol X, Kush A. Towards a latex molecular diagnostic of yield potential and the genetic engineering of the rubber tree. Euphyt. 1997;96:77–82.

  68. 68.

    Grabherr MG, Haas BJ, Levin JZ, Thompson DA, Amit I, Adiconis X, Fan L, Raychowdhury R, Zeng Q, Chen Z, Mauceli E, Hacohen N, Gnirke A, Rhind N, di Palma F, Birren BW, Nusbaum C, Lindblad-Toh K, Friedman N, Regev A. Full-length transcriptome assembly from RNA-Seq data without a reference genome. Nat Biotech. 2011;29(7):644–52.

  69. 69.

    Iseli C, Jongeneel CV, Bucher P. ESTScan: a program for detecting, evaluating, and reconstructing potential coding regions in EST sequences. Proc Int Conf Intell Syst Mol Biol. 1999;138-148.

  70. 70.

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

  71. 71.

    Ye J, Fang L, Zheng H, Zhang Y, Chen J, Zhang Z, Wang J, Li S, Li R, Bolund L, Wang J. WEGO: a web tool for plotting GO annotations. Nucleic Acids Res. 2006;34(Web Server issue):W293–297.

  72. 72.

    Kanehisa M, Araki M, Goto S, Hattori M, Hirakawa M, Itoh M, Katayama T, Kawashima S, Okuda S, Tokimatsu T, Yamanishi Y. KEGG for linking genomes to life and the environment. Nucleic Acids Res. 2008;36(Database issue):D480–484.

  73. 73.

    Mortazavi A, Williams BA, McCue K, Schaeffer L, Wold B. Mapping and quantifying mammalian transcriptomes by RNA-Seq. Nat Methods. 2008;5:621–8.

  74. 74.

    Audic S, Claverie JM. The significance of digital gene expression profiles. Genome Res. 1997;7:986–95.

  75. 75.

    Reiner A, Yekutieli D, Benjamini Y. Identifying differentially expressed genes using false discovery rate controlling procedures. Bioinformatics. 2003;19:368–75.

  76. 76.

    Yednock BK, Sullivan TJ, Neigel JE. De novo assembly of a transcriptome from juvenile blue crabs (Callinectes sapidus) following exposure to surrogate Macondo crude oil. BMC Genomics. 2015;16:521.

  77. 77.

    Jin J, Zhang H, Kong L, Gao G, Luo J. PlantTFDB 3.0: a portal for the functional and evolutionary study of plant transcription factors. Nucleic Acids Res. 2014;42:D1182–1187.

  78. 78.

    Venkatachalam P, Thulaseedharan A, Raghothama K. Molecular identification and characterization of a gene associated with the onset of tapping paneldryness (TPD) syndrome in rubber tree (Hevea brasiliensis Muell.) by mRNA differential display. Mol Biotechnol. 2009;41(1):42–52.

Download references


This work was funded by the Hainan Province Major Science and Technology Project (ZDZX2013023), the Western Plan and Subject Key Areas Construction Project of Hainan University (ZXBJH-XK001).

Author information

Correspondence to Jin-Ping Liu.

Additional information

Competing interests

The authors declare that they have no competing interests.

Authors’ contributions

LJP supervised the experiments, conducted the bioinformatics studies and wrote the manuscript; ZYF carried out qRT-PCR validation; GXL and LYJ participated in collection of plant materials. All the authors read and approved the final manuscript.

Additional files

Additional file 1:

Summary of functional annotations of unigenes by aligning to the Nr, Nt, Swiss-Prot, COG, GO and KEGG databases. (XLS 24813 kb)

Additional file 2:

List of DEGs between C and E-8, and C and E-24. (XLS 23067 kb)

Additional file 3:

The top 20 most up-regulated and down regulated genes between C and E-8, and C and E-24. (XLS 43 kb)

Additional file 4:

Pathways significantly enriched of DEGs (C vs. E-8, and C vs. E-24). (XLS 196 kb)

Additional file 5:

DEGs involved in the glycolysis in E8 and E24 compared to C samples of Hevea brasiliensis. (XLS 59 kb)

Additional file 6:

Differential expression of unigenes involved in IPP biosynthesis pathway in E8 and E24 compared to C samples of Hevea brasiliensis. AACT: acetoacetyl-CoA thiolase; HMGS: 3-hydroxy-3-metylglutaryl coenzyme A (HMG-CoA) synthase; HMGR: HMG-CoA reductase; MVK: mevalonate (MVA) kinase; PMVK: 5-phosphomevalonate (MVP) kinase; PMD: phosphatemevalonate decarboxylase; DPMD: 5-diphosphomevalonate (MVPP) decarboxylase; IPK: isopentenyl phosphate (IP) kinase; GPS: geranyl diphosphate synthase; IDI: isopentenyl diphosphate (IPP) isomerase; DXS: 1-deoxy-d-xylulose 5-phosphate (DXP) synthase; DXR: DXP reductoisomerase; MCT: 4-(cytidine 5′ -diphospho)-2- C- methyl-D-erythritol (CDP-ME) synthase; CMK: CDP-ME kinase; MDS: 2C-methyl-D-erythritol 2,4-cyclodiphosphate (MEcPP) synthase; HDS: 4-hydroxy-3-methylbut-2-enyl diphosphate (HMBPP) synthase; HDR: HMBPP reductase. Cells with gray border lines in the upper rows represent differentially expressed unigenes in E8 compared to C and cells with green border lines in the lower rows represent differentially expressed unigenes in E24 compared to C. Relative levels of expression are showed by a color gradient from low (blue) to high (red). (JPG 157 kb)

Additional file 7:

DEGs involved in the IPP biosynthesis in E8 and E24 compared to C samples of Hevea brasiliensis. (XLS 31 kb)

Additional file 8:

DEGs involved in the Calvin cycle in E8 and E24 compared to C samples of Hevea brasiliensis. (XLS 43 kb)

Additional file 9:

Differential expression of unigenes involved in hormone signaling in E8 and E24 compared to C samples of Hevea brasiliensis. Ethylene signalling pathway: ETR1: ETHYLENE RESPONSE 1; CTR1: CONSTITUTIVE TRIPLE RESPONSE 1; EIN2: ETHYLENE INSENSITIVE 2; EIN3: ETHYLENE INSENSITIVE 3; ERF1/2: ETHYLENE RESPONSE FACTOR 1/2; EBF1/2: EIN3 binding F-Box protein 1/2; BR signaling pathway: BRI1: Brassinosteroid-Insensitive 1; BAK1: BRI1-associated kinase 1; BKI1: BRI1 KINASE INHIBITOR 1; BSK: BR SIGNALING KINASE; BSU1: bri1 SUPPRESSOR 1; BIN2: BRASSINOSTEROID-INSENSITIVE 2; BZR1/2: BRASSINAZOLE RESISTANT 1/2; TCH: TOUCH genes; CYCD3: CYCLIN D3; GA signaling pathway: GID1: GIBBERELLIN INSENSITIVE DWARF 1; GID2: GIBBERELLIN INSENSITIVE DWARF 2; DELLAs: DELLA growth inhibitors; TF: transcriptional factor; Auxin signaling pathway: AUX1: AUXIN1; TIR1: TRANSPORT INHIBITOR RESPONSE 1; IAA: INDOLE ACETIC ACID; ARF: AUXIN RESPONSE FACTOR; SAUR: Small Auxin-Up RNA; G10H: geraniol 10-hydroxylase gene; Cytokinin signaling pathway: CRE1: CYTOKININ RESPONSE 1; AHP: histidine phosphotransfer protein; B-ARR: type-B response regulator (ARR); A-ARR: type-A response regulator (ARR); SA signalling pathway: NPR1: Non-expressor of pathogenesis-related genes 1; TGA: the bZIP transcription factors; PR1: pathogenesis related protein 1; JA signaling pathway: JAR1: JASMONATES RESISTANT 1; JA-Ile: jasmonoyl isoleucine; JAZ: Jasmonate ZIM-domain-containing protein; MYC2: a basic helix-loop-helix (bHLH) transcription factor; ORCA3: Octadecanoid-derivative Responsive Catharanthus AP2-domain gene; ABA signalling pathway: PYR1/PYLs: Pyrabactin Resistance Protein1/PYR-Like proteins; PP2Cs: protein phosphatases which fall under the category of type 2C; SnRK2: SNF1 (Sucrose-Nonfermenting Kinase1)-related protein kinase 2: ABF: ABA responsive element (ABRE) binding factors. Cells with gray border lines in the upper rows represent differentially expressed unigenes in E8 compared to C and cells with green border lines in the lower rows represent differentially expressed unigenes in E24 compared to C. Relative levels of expression are showed by a color gradient from low (blue) to high (red). (JPG 249 kb)

Additional file 10:

DEGs of hormone signaling components in E8 and E24 compared to C samples of Hevea brasiliensis. (XLS 126 kb)

Additional file 11:

DEGs of transcription factors in E8 and E24 compared to C samples of Hevea brasiliensis. (XLS 256 kb)

Additional file 12:

The top five up-regulated and down-regulated TFs in E8 and E24 compared to C. (XLS 23 kb)

Rights and permissions

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

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Liu, J., Zhuang, Y., Guo, X. et al. Molecular mechanism of ethylene stimulation of latex yield in rubber tree (Hevea brasiliensis) revealed by de novo sequencing and transcriptome analysis. BMC Genomics 17, 257 (2016).

Download citation


  • Rubber tree
  • Hevea brasiliensis
  • Ethephon treatment
  • Ethylene stimulation
  • Rubber production