Genome-wide identification of rubber tree (Hevea brasiliensis Muell. Arg.) aquaporin genes and their response to ethephon stimulation in the laticifer, a rubber-producing tissue

Background Natural rubber, an important industrial raw material, is specifically synthesized in laticifers located inside the rubber tree (Hevea brasiliensis Muell. Arg.) trunk. Due to the absence of plasmodesmata, the laticifer water balance is mediated by aquaporins (AQPs). However, to date, the characterization of H. brasiliensis AQPs (HbAQPs) is still in its infancy. Results In this study, 51 full-length AQP genes were identified from the rubber tree genome. The phylogenetic analysis assigned these AQPs to five subfamilies, including 15 plasma membrane intrinsic proteins (PIPs), 17 tonoplast intrinsic proteins (TIPs), 9 NOD26-like intrinsic proteins (NIPs), 4 small basic intrinsic proteins (SIPs) and 6 X intrinsic proteins (XIPs). Functional prediction based on the analysis of the aromatic/arginine (ar/R) selectivity filter, Froger’s positions and specificity-determining positions (SDPs) showed a remarkable difference in substrate specificity among subfamilies. Homology analysis supported the expression of 44 HbAQP genes in at least one of the examined tissues. Furthermore, deep sequencing of the laticifer transcriptome in the form of latex revealed a key role of several PIP subfamily members in the laticifer water balance, and qRT-PCR analysis showed diverse expression patterns of laticifer-expressed HbAQP genes upon ethephon treatment, a widely-used practice for the stimulation of latex yield. Conclusions This study provides an important genetic resource of HbAQP genes, which will be useful to improve the water use efficiency and latex yield of Hevea. Electronic supplementary material The online version of this article (doi:10.1186/s12864-015-2152-6) contains supplementary material, which is available to authorized users.

P. trichocarpa, one more subfamily named X intrinsic protein (XIP) subfamily that contains three subgroups is also found [7,12].
AQPs assemble in tetramers in the cell membrane, although each monomer can act as a water channel [13]. Even though the overall pairwise sequence similarity can be low, AQPs share some structural features such as harboring six transmembrane helices (TM1-TM6) connected by five loops (LA-LE). LB and LE from opposite sides dip into the membrane and form two half helices (HB and HE), at the N-termini of which, two highly conserved NPA (Asn-Pro-Ala) motifs form one selectivity region. And another region that determines the substrate specificity is known as the aromatic/arginine (ar/R) selectivity filter (H2 in TM2, H5 in TM5, LE1 and LE2 in LE) [13]. The NPA motifs create an electrostatic repulsion of protons and act as a size barrier, where the ar/R filter renders the pore constriction site diverse in both size and hydrophobicity [13]. In addition to water-conducting AQPs, certain AQPs called aquaglyceroporins (GLPs) transport glycerol instead. Statistical analysis indicated that GLPs feather five highly conserved amino acid residues (named Froger's positions: P1-5): an aromatic residue at P1, an acidic residue at P2, a basic residue at P3, a proline followed by a nonaromatic residue at P4 and P5, as Y 108 -D 207 -K 211 -P 236 -L 237 observed in the Escherichia coli glycerol facilitator GlpF in contrast to A 103 -S 190 -A 194 -F 208 -W 209 in the pure water channel AqpZ [14]. Very recently, nine specificity-determining positions (SDPs) for non-aqua substrates, i.e. urea, boric acid, silicic acid, ammonia (NH 3 ), carbon dioxide (CO 2 ) and hydrogen peroxide (H 2 O 2 ) were also proposed for each group via a comprehensive analysis of functionally characterized AQPs [15].
The para rubber tree (Hevea brasiliensis Muell. Arg.) is a perennial tropical plant species that belongs to the Euphorbiaceae family. Although it is native to the Amazon basin, the economic importance and increasing demand of natural rubber has prompted its wide-domestication to Southeast Asia. To date, the rubber tree is still the only commercial source of natural rubber for its high production and quality of rubber. Natural rubber (cis-1,4-polyisoprene) is specifically synthesized in the highly differentiated vessels termed laticifers, which are located in the secondary phloem of the tree trunk and are periodically differentiated from the vascular cambium [16]. The rubbercontaining latex which represents the cytoplasmic content of laticifers is harvested by tapping the bark every 2-5 days, and the latex yield is determined by the initial flow rate, duration of latex flow and latex regeneration between two tappings [17].
Ethylene, a gas phytohormone, plays crucial roles in numerous aspects of growth, development, and response to biotic and abiotic stresses in plants [18]. In Arabidopsis, the ethylene signaling pathway has been well established and its involvement in certain agronomically important processes such as seed dormancy, fruit ripening, abscission and senescence has made ethylene a target for manipulation by chemical and biotechnological methodologies [19]. Ethephon (also known as Ethrel), an ethylene releaser, was initiatively tested for rubber yield promotion as early as the 1970s and now is widely used in rubber production all over the world. Although the molecular mechanism on ethephon stimulation of latex yield is still unclear, researches showed that the treatment of rubber tree barks with ethephon could significantly decrease latex dry rubber content (DRC) or total solid content (TSC), extend the bark drainage area, and prolong latex flow duration [20][21][22][23]. These effects are mainly benefited from water influx toward laticifers and the resultant latex dilution.
Lately, the rubber tree genome was sequenced by Rahman et al. [30] and Rubber Research Institute, Chinese Academy of Tropical Agricultural Sciences. In addition, more than 50,000 ESTs (expressed sequence tags) and a high number of RNA sequencing reads derived from several tissues such as shoot apex, leaf, laticifer, bark, root and somatic embryogenesis are also available in NCBI SRA [30][31][32][33][34][35]. These datasets provide a good chance to analyze the rubber tree AQP gene family from a global view. In the present study, a genome-wide search was carried out to identify AQP genes encoded in the rubber tree genome. Further, functional prediction was performed based on the analysis of the ar/R filter, Froger's positions and SPDs, and deep transcriptome sequencing and qRT-PCR expression analysis were also adopted to identify the most important AQP members expressed in the laticifer.

Identification and classification of rubber tree aquaporin genes
Via a comprehensive homology analysis, 57 or 54 loci putatively encoding AQP-like genes were identified from the rubber tree genome of clone RY7-33-97 or RRIM600, respectively (data not shown). Since all AQP-encoding loci identified in the RRIM600 genome were found in the genome of RY7-33-97 and some genes from RRIM600 are incomplete and/or the sequences have a high number of "N"s, the AQP genes identified from the RY7-33-97 genome were selected for further analyses. After discarding loci encoding partial AQP-like sequences which are truncated and lacking any of the NPA motifs, 51 fulllength AQP genes were retained and the gene models are available in Additional file 1.

Analysis of exon-intron structure
The exon-intron structures of the 51 HbAQP genes were analyzed based on the gene models. Although the ORF (open reading frame) length of each gene is similar (684-927 bp), the gene size (from start to stop codons) is distinct (720-13833 bp, Table 1; Fig. 2). The introns of HbAQP genes harbor an average length of 404 bp, with the minimum of 71 bp in HbNIP2;1 and the maximum of 13000 bp in HbSIP2;1. Genes in different subfamilies harbor distinct exon-intron structures. All members of the HbPIP subfamily feature three introns (92-736 bp, 78-1650 bp and 80-186 bp, respectively). Except for HbTIP1;1, HbTIP1;2, HbTIP1;3 and HbTIP1;4 that contain only one intron, other HbTIP genes contain two introns instead. HbNIP genes usually have four introns except for HbNIP5;1 containing three introns. Most HbSIP genes contain two introns except for HbSIP1;1 without any intron. Subgroups of HbXIP genes vary in the number of introns: one intron for subgroup one, two or zero for subgroups two and three, respectively (Fig. 2).

Structural features of HbAQPs
Sequence analysis showed that the 51 deduced HbAQPs consist of 227-305 amino acids, with a theoretical molecular weight of 23.78-32.28 kDa and a pI value of 4.59-9.74. Homology analysis revealed a high sequence diversity existing within and between the five subfamilies. The sequence Topological analysis showed that all HbAQPs were predicted to harbor six transmembrane helical domains (Table 2), which is consistent with the results from multiple alignments with structure proven AQPs (see Additional file 5). The subcellular localization of each HbAQP was also predicted ( Table 2). HbPIPs with an average pI value of 8.47 are localized to plasma membranes. HbTIPs with an average pI value of 5.81 are mainly localized to vacuoles (known as lutoids in laticifers with a natural pH of about 6), though several members were predicted to be localized to endoplasmic reticulum (ER), chloroplast and cytosol. HbNIPs with an average pI value of 7.60 are mostly localized to plasma membranes, but HbNIP2;1 and HbNIP3;1 were predicted to be localized to the membrane of vacuole and chloroplast, respectively. Two members (HbSIP1;2 and HbSIP1;3) of the SIP subfamily (with an average pI value of 9.06) were predicted to be localized to plasma membranes, whereas HbSIP1;1 and HbSIP2;1 are localized to the membrane of vacuole and chloroplast, respectively. Although the XIP subfamily harbors only six members (with an average pI value of 7.95), the predicted localizations are diverse, including the vacuole, chloroplast, plasma membrane and cytosol. To learn more about the putative function of HbAQPs, the conserved residues typical of dual NPA motifs, the ar/R filter, five Froger's positions and nine SDPs were also identified (Tables 2 and 3).

HbPIP subfamily
All HbPIPs were identified to have similar sequence length, however, HbPIP2s (278-288 residues) can be distinguished from HbPIP1s (269-287 residues) by harboring relatively shorter N-terminal and longer Cterminal sequences (Additional file 5). The five HbPIP1s have sequence similarities of 79.2-99.3 %, whereas the similarity percents of ten HbPIP2s are 77.4-98.3 %. Between HbPIP1 and HbPIP2 members, sequence similarities of 59.1-65.9 % are observed (Additional file 4). The dual NPA motifs, ar/R filter (F-H-T-R), and four out of five Froger's positions are highly conserved in HbPIPs (Table 2). In contrast, the P1 position is more variable with the appearance of an E, Q or M residue ( Table 2). In addition, two phosphorylation sites corresponding to S115 and S274 in SoPIP2;1 [13] are invariable in HbPIP2s, and the former one is even highly conserved in all HbPIPs, HbTIPs and HbXIPs except for the S → T substitution in several members (Additional file 5), implying their regulation by phosphorylation.

HbXIP subfamily
HbXIPs vary from 256 to 305 residues in length (   (Table 2). Similar to most XIPs [7], two highly conserved C residues in the LGGC motif of LC and the NPARC motif of LE were also found in HbXIPs except for HbXIP1;2, in which the F residue is located at the corresponding position of LE (Additional file 5).

Transcriptional profiles of HbAQP genes in the laticifer and their response to ethephon stimulation
The rubber tree laticifer is a single-cell-type tissue specifically for natural rubber biosynthesis. To identify the AQP genes expressed in the laticifer and determine the most important members in the laticifer water balance, the latex representing the laticifer cytoplasm was collected and high-quality total RNAs (260/280 value between 1.95 and 2.00, 28S/18S value between 3.0 and 3.2 and RIN value between 8.9 and 9.1) were isolated from three biological replicates, respectively.

High abundance and diversity of HbAQPs
A total of 51 full-length AQP genes were identified from the rubber tree genome, which is comparable to 55 members reported in poplar (a tree species also belongs to Malpighiales) [7,36]; more than 23 in grapevine [6], 33 in rice [5], 35 in Arabidopsis [3], 36 in maize [4], 41 in potato [10] and 47 in tomato [9]; less than 66 in soybean [11] and 71 in cotton [8]. Since the AQP genes in Arabidopsis and poplar were well characterized, their deduced proteins were added in the phylogenetic analysis of HbAQPs, which assigned them to five subfamilies. With the exception of the XIP subfamily, the further classification of HbAQP subfamilies into subgroups is consistent with Arabidopsis, i.e. two PIP subgroups, five TIP subgroups, seven NIP subgroups and two SIP subgroups. Nevertheless, classing AtNIP2;1 and AtNIP3;1 into the NIP1 subgroup was proposed. As shown in Fig. 1, AtNIP2;1 and AtNIP3;1 were clustered with the NIP1 subgroup, sharing the highest similarity with Fig. 4 qRT-PCR analysis of 19 laticifer-expressed HbAQP genes upon ethephon stimulation. Data are mean ± SD (n = 9). Different letters mean significant difference over the time course AtNIP1;2 in Arabidopsis, HbNIP1;2 or HbNIP1;1 in rubber tree, PtNIP1;2 or PtNIP1;1 in poplar, respectively. Thereby, no NIP2s and NIP3s were retained in Arabidopsis as seen in rubber tree and poplar (Fig. 5). Since no XIP homologs were found in the Arabidopsis genome, the nomenclature for poplar proposed by Lopez et al. [36] was adopted to divide HbXIPs into three subgroups. Besides supported by high bootstrap values, XIP1s are characterized by the ar/R filter of V-M-V/P/ A-R, XIP2s by I-I-V-R and XIP3s by V-K-A-R.
Gene pairs were identified not only in rubber tree, but also in poplar and Arabidopsis (Fig. 1). For example, five AtPIP1s were clustered together apart from PIP1s of rubber tree and poplar; HbPIP1;1 and HbPIP1;2 were clustered with PtPIP1;1 and PtPIP1;2. These results suggest the occurrence of more than one gene duplication events. Previous studies indicated that poplar underwent one whole-genome triplication event (designated 'γ') and one doubling event, whereas Arabidopsis underwent the same γ event and two independent doubling events, though the Arabidopsis genome encodes relatively less AQP genes due to massive gene loss and chromosomal rearrangement after genome duplications [40][41][42]. The γ duplication occurred at approximate 117 million years ago, shortly before the origin of core eudicots [43]. As a core eudicot plant, the rubber tree appears to share the γ duplication. However, another one as the data suggested is likely to be a doubling event independent from both Arabidopsis and poplar, probably occurred after the divergence of Euphorbiaceae and Salicaceae. A genome-wide comparative analysis may provide more information.

Functional inference of HbAQPs
Although plant AQPs firstly raised considerable interest for their high water permeability, when heterologously expressed in Xenopus oocytes or yeast cells, increasing evidence has shown that some of them are also participated in the transport of other small molecules such as glycerol, urea, boric acid, silicic acid, NH 3 , CO 2 and H 2 O 2 [44]. Based on atomic resolution structures and molecular dynamics stimulations of GlpF, AqpZ, AQP1 and other MIPs, several structural features determining their transport selectivity were identified, e.g. the two opposite NPA motifs, the ar/R filter and the amino acid residues at Froger's positions for discriminating between AQPs and GLPs [14,45]. As shown in Table 2, most HbAQPs exhibit an AqpZ-like Froger's positions to favor the permeability of water. In contrast, HbSIP2;1 and NIP subfamily members possess mixed key residues of GlpF for P1 and P5, and AqpZ for P2-P4. The glycerol permeability of GmNOD26 and Arabidopsis NIPs was reported [46,47], however, the potential glycerol transport ability of HbSIP2;1-like SIPs have not be confirmed by experimental means yet.
In addition to high permeability to water, plant PIPs were reported to transport urea, boric acid, CO 2 and H 2 O 2 [48]. As shown in Table 2, all HbPIPs represent the F-H-T-R ar/R filter as observed in AqpZ which harbors an extremely narrow and hydrophilic pore (diameter 2.8 Å) [45], suggesting their high water permeability. However, when expressed in Xenopus, extremely low water permeability of HbPIP1 members such as HbPIP1;1 and HbPIP1;4 was observed [22,26] as seen in many other plant species [49]. Based on the SDP analysis proposed by Hove and Bhave [15], all HbPIPs represent urea-type SDPs  (Table 3), supporting their similar functionality.
Although highly variable in the ar/R filter, plant TIPs were shown to transport water as efficiently as PIPs [21]. Additionally, they also allow urea, NH 3 and H 2 O 2 through [50]. As shown in Table 3, all HbTIPs except for  HbTIP2;4 and HbTIP4;1 represent urea-type SDPs (H-P Besides glycerol and water, plant NIPs have been found to transport urea, boric acid, silicic acid, NH 3 and H 2 O 2 [50][51][52]. As shown in Table 3 and Additional file 6, HbNIP5;1 is promised to be a urea and boric acid transporter with nine SDPs of H-P-I-A-L-P-G-S-N or T-I-H-P-E-L-L-A-P. HbNIP2;1 represent typical urea SDPs (H-P-T-A-M-P-G-S-N), and SDPs of V-V-H-P-E-I-I-A-P with the substitution of V for I at SDP2 in comparison to typical boric acid SDPs (T/V-I-H-P-E-I/L-I/L-A/T-A/G/P). Compared with typical urea and boric acid SDPs, HbNIP6;1 seems to represent novel SDPs-types with the substitution of Q for A/P at SDP6 or Q for A/P/G at SDP9. Although characterized as an NIP III member, the silicic acid transport ability of HbNIP2;1 needs to be experimentally validated since it seems to represent novel SDPs (S-F-V-H-G-N-R-T-Q in contrast to typical C/S-F/Y-A/E/L-H/R/Y-G-K/N/T-R-E/S/T-A/K/P/T) similar to that of GmNIP2;1 and GmNIP2;2 (S-Y-E-R-G-N-R-T-P) [53]. Although GmNOD26 was reported to transport NH 3 [54], whether its close rubber tree homologs (i.e. HbNIP1;1, HbNIP1;2, HbNIP3;1, HbNIP4;1 and HbNIP4;2) represent novel SDPs-types still needs to be tested. HbNIP3;1, HbNIP4;2 and HbNIP5;1 represent H 2 O 2 -type SDPs (A/S-A-L-L/V-I/V-L-Y-V-P) slightly different from AtNIP1;2 (S-A-L-L-V-L-Y-V-P) [50].
A crucial role of HbPIPs in the water balance of laticifers As a unique site for rubber biosynthesis, the laticifers are present in a wide variety of rubber tree tissues, including shoots, roots, stems, leaves, flowers, fruits, cotyledons, inner seed coats, etc., and can be divided into primary and secondary laticifers according to their origin [16]. Compared with the procambium-derivation of primary laticifers, the secondary laticifers, mainly located in the soft inner bark of the rubber tree trunk, are periodically differentiated from the vascular cambium and serve as a sole source for the commercial latex [56]. During the differentiation and maturation process, laticifer mother cells articulate with each other and further anastomose together into a successive vertical network (called rings or mantles) arranged as concentric sheaths in the secondary phloem [56]. Unlike other cells such as neighboring parenchyma cells, the mature laticiferous cells are totally devoid of plasmodesmata [25], and thus its water exchanges with surrounding cells are mainly governed by AQPs. Upon bark tapping, the laticifer cytoplasm is expelled in the form of latex due to the high turgor pressure inside [57]. Generally, latex flow can continue for several hours until coagulation processes lead to the plugging of severed laticifers [58]. During the latex flow, a progressive decrease in DRC was observed [21][22][23], indicating rapid water influx and latex dilution inside laticifers caused by the activity of HbAQPs. Given that of HbPIPs and HbTIPs account for more than 62.7 % of the total HbAQPs and their AqpZ-like Froger's positions favoring the high water permeability, we initially prospect that these two subfamilies may play important roles in the laticifer water balance: the plasma membrane-targeted HbPIPs facilitate the water transport from the extracellular space to the laticifer cytoplasm, whereas the lutoid-targeted HbTIPs play an essential role in maintaining the cell osmotic balance as observed in most plant cells [59]. However, in contrast to the mature plant cells characterized by a large central vacuole which occupies 80 % or more of the intracellular space, the lutoids in laticifers are polydispersed microvacuoles occupying only 12 % of the total latex [60], arguing the central role of HbTIPs in the laticifer water balance, though their potential role in the lutoid stability and latex vessel plugging should be noted. To address this issue, the transcriptome of such a single-cell-type tissue was deeply sequenced. Results showed that PIP members were the main AQP genes expressed in the laticifer (similar results were also observed when the recently available laticifer transcriptome of clone RRIM928 was analyzed, see Additional file 7), suggesting their crucial role, especially the highly abundant HbPIP2;7, HbPIP1;4 and HbPIP2;5, in the laticifer water balance. When expressed in Xenopus oocytes, our previous study showed that HbPIP2;5 could transport water as efficiently as HbPIP2;1 [21,23]; in contrast, HbPIP1;4 and HbPIP2;7 were shown to be less efficient [22]. In addition, as a PIP1 member, the poor efficiency of HbPIP1;1 was also observed [26]. Therefore, the exact role of HbPIP2;7, HbPIP1;4 and HbPIP2;5 in the water balance of rubber tree laticifers needs further investigations.
To profile the AQP genes in response to ethylene stimulation in laticifers, the latex at different time points after ethephon treatment was collected from rubber tree clone PR107. Similar to PB217, PR107 clone is characterized as a relatively late mature variety which has a high TSC, short latex flow duration and low latex metabolism, however, ethephon stimulation could significantly prolong its latex flow duration and enhance latex yield [22,23]. Our qRT-PCR analysis showed that the expression levels of most laticifer-expressed genes significantly changed at least one tested time point after ethephon application (Fig. 4), indicating their involvements in the ethephon enhanced water influx into laticifers. Among these time points, the latex collected at 24 and 40 h (especially 24 h) after ethephon treatment was shown to harbor the most abundant transcripts, which include four of the five highly abundant HbPIP1;3, HbPIP2;3, HbPIP1;4 and HbPIP2;5, corresponding to the significantly decreased TSC, the longest latex flow duration and the highest latex yield as reported by Wang et al. who utilized the same materials [22]. Besides, similar effects of ethephon on latex yield and latex TSC of the PB217 clone were also observed by Tungngoen et al., although they used mature virgin trees as materials [21].

Conclusions
To our knowledge, this is the first genome-wide study of the rubber tree AQP gene family and using systematic nomenclature assigned 51 HbAQPs into five subfamilies based on the sequence similarity and phylogenetic relationship with their Arabidopsis and poplar counterparts. Furthermore, their structural and functional properties were investigated based on the analysis of the ar/R filter, Froger's positions and SPDs, which suggested the potentially key role of HbPIPs and HbTIPs in the laticifer water balance. Most importantly, the laticifer transcriptome was deeply sequenced to identify the most important AQPs in such a single-cell-type tissue, and qRT-PCR analysis was also performed to investigate the expression profiles of laticifer-expressed HbAQP genes upon ethephon stimulation. Our results revealed that HbPIPs were the mainly AQP genes expressed in the laticifer. Among 19 HbAQP genes detected in the laticifer, most of them were significantly regulated by ethylene. Consistent with the significantly decreased TSC and increased latex yield, most laticifer-expressed PIP genes were considerably induced at the time point of 24 h after ethephon application, supporting their crucial roles in the water balance of laticifers in the case of ethephon stimulation. This study provides an important genetic resource of HbAQP genes, which will be useful to improve the water use efficiency and latex yield of Hevea.

Identification of rubber tree aquaporin genes
The deduced amino acids of HbAQPs available in the NCBI GenBank were used as queries to search the available RRIM600 genome and our in-house RY7-33-97 genome for rubber tree homologs. Sequences with an Evalue of less than 1e −5 in the tBlastn search [61] were selected for further analyses. The gene structures were firstly predicted using GeneMark.hmm [62], and the gene models were further validated with ESTs and raw RNA sequencing reads available at GenBank. The exonintron structures of AQP genes detected in the laticifer transcriptome were also confirmed by aligning the cloned cDNAs to the corresponding gene sequences. Gene structures were displayed using GSDS [63]. Homology search for nucleotides or Sanger ESTs was performed using Blastn, and sequences with an identity of more than 98 % were taken into account. RNA sequencing reads were mapped using Bowtie 2 [64] with default parameters, and mapped read number of more than one was counted as expressed. Unless specific statements, the tools used in this study were performed with default parameters.

Sequence alignments and phylogenetic analysis
Multiple sequence alignment using deduced proteins was performed with ClustalX [65], and the unrooted phylogenetic tree was constructed by the maximum likelihood method using MEGA6 [66]. The reliability of branches in the resulting tree was supported with 1,000 bootstrap resamplings. Classification of AQPs into subfamilies and subgroups was done as described before [3,36].

Plant materials and field experiments
PR107, the male parent of rubber tree clone RY7-33-97, was planted at the experimental farm of Chinese Academy of Tropical Agricultural Sciences (Danzhou, China) in 2002. Six batches of three trees with similar growth performance and latex yield were selected for this study. The trees had been tapped for 3 years on the s/2 d 3 system (tapping every 3 days with half spiral) without ethephon stimulation. For ethephon stimulation, five batches of trees were treated with 1 g of 2.5 % (w/w) ethephon in carboxyl methyl cellulose (CMC, 1 %) for 6, 16, 24, and 40 h before the sampling. The sixth batch was treated with 1 % CMC as a control.

Latex collection and total RNA extraction
The latex was collected through tapping the bark at around 6:00 am, and samples representing three biological replicates were subjected for total RNA isolation as described by Tang et al. [71]. Briefly, the latex within the first 45 min was dropped into liquid nitrogen after discarding the first 5 drops. The frozen latex was suspended with extraction buffer (0.3 M LiCl, 0.01 M disodium salt EDTA, 10 % (W/V) SDS, 0.1 M Tris-HCl), and equal volume of water-saturated phenol/chloroform/isoamyl alcohol (PCI) (25:24:1) was added and vigorously shaken. Then, the mixture was centrifuged at 12,000 × g for 10 min at 4°C, and the aqueous phase was collected and subjected to one more PCI and one chloroform/isoamyl alcohol (24:1) extraction. The supernatant was precipitated with 8 M LiCl solution for twice. The pellet was dissolved with H 2 O, and 3 M NaAc (pH 5.2) and absolute alcohol were added to precipitate the RNA. After washed with 75 % ethanol, the RNA was dissolved with H 2 O. The concentration and integrity of total RNA was confirmed using a 2100 Bioanalyzer (Agilent, Palo Alto, CA, USA).
For the expression analysis, the first-strand cDNA was synthesized from 2 μg of total RNA to a final 20 μL reaction mixture using PrimeScript® RT reagent kit with gDNA Eraser (Takara, Dalian, China) according to the manufacture's instruction, and then stored at −20°C.
For Illumina sequencing, magnetic beads with biotin-Oligo (dT) were used to isolate poly(A) mRNA according to the manufacturer's protocol of Illumina TruSeqTM RNA sample preparation kit (Qiagen GmbH, Hilden, Germany).

Expression analysis based on Illumina sequencing
RNA sequencing was performed as described previously [72] using Illumina HiSeq™ 2000 (Illumina Inc., San Diego, CA, USA) at Beijing Genomics Institute (Shenzhen, China). The raw data were filtered by the Illumina pipeline to remove adaptor sequences, adaptor-only reads, reads with "N" rate larger than 10 % ("N" representing ambiguous bases) and low quality reads containing more than 50 % bases with Q-value ≤ 5. Assembly of clean reads was carried out using SOAP de novo [73] (Luo et al. 2012). The trimmed reads were mapped to Unigenes using Bowtie 2 [64], and the RPKM (reads per kilo bases per million reads) method [74] was used for the expression annotation.

qRT-PCR analysis
HbYLS8, the most stably expressed genes in response to ethephon stimulation [75], was selected as the reference gene in this study. The gene-specific primers are listed in Additional file 8, and the PCR reaction was performed using the SYBR-green Mix (Takara, Dalian, China) and the Real-time Thermal Cycler (Type 5100, Thermal Fisher Scientific Oy, Finland). All qRT-PCR assays were performed in triplicate for each biological sample. The amplification efficiency of each primer pair was estimated via melting curve analysis, and PCR products were confirmed by Sanger sequencing. The relative abundance of each transcript was estimated with the 2 ΔΔCt method after normalization against HbYLS8 using PikoReal2.0 software unless otherwise specified. Statistical analyses were executed using the Data Processing System software v11.0. The differences among means were tested following Duncan's one-way multiple-range post hoc ANOVA (P < 0.05).