Identification of novel microRNAs in Hevea brasiliensisand computational prediction of their targets

Classification of small RNAs

Annotation of the Hevea small RNAs was attempted by BLASTN on the Arabidopsis genome. Most of the sequences did not map onto the Arabidopsis genome sequence. Of the mapped sequences, 26,721 were not annotated. Annotated sequences predicted 6,618 products of mRNA degradation, 6,800 transposable elements, 399 rRNAs, 138 products of pseudogene degradation, 73 ncRNAs, 39 tRNAs, 5 sn/snoRNAs and 57 miRNAs corresponding to 14 families.

Further mapping was carried out with all the miRNA sequences available in the PMRD database. In all, 18,430 sequences matched (0 or 1 mismatch) to conserved miRNAs in other species. The accessions were classified into 48 miRNA families (Table 1). Since most sequences did not map against PMRD, the Hevea EST database was scanned using the LeARN pipeline. This strategy led to the identification of 10 additional miRNA families, based on the stem-loop structure, which were not described in other species (Table 1). The distribution of the 670,645 unique accessions showed that the most abundant accessions had in descending order 27, 19, 24, 26 and 17 nucleotides (Figure 1). Accessions with 17 and 19 nucleotides came from products of mRNA degradation and transposable elements according to the mapping results against the Arabidopsis genome sequence (data not shown).

Characteristics of conserved and putative novel miRNA families in Hevea brasiliensis

Reads of conserved and putative novel miRNAs were distributed according to their nucleotide lengths (Figure 2). The 48 conserved miRNA families mostly had 24 and 23 nucleotides (Figure 2A). For the 10 putative novel miRNA families, the most represented lengths were 24, 27 and 26 nucleotides (Figure 2B).

Of the 48 conserved miRNA families, 12 were predominant with more than 1,000 reads, including 3 families that displayed more than 20,000 reads. The latter corresponded to HbmiR159/319 and HbmiR408 (Figure 3A). By contrast, 36 families were less represented in the dataset with 18 families having fewer than 100 reads: HbmiR1310, HbmiR162, HbmiR168, HbmiR1863, HbmiR2118, HbmiR2910, HbmiR2914, HbmiR2915, HbmiR2916, HbmiR393, HbmiR394, HbmiR395, HbmiR399, HbmiR444, HbmiR476, HbmiR482, HbmiR828 and HbmiR845 (Figure 3B).

A detailed analysis of length distribution within the Hevea miRNA families revealed some discrepancies (Figure 4 and Figure 5). Firstly, the 20-22 nucleotide and 23-27 nucleotide size classes were detected for both conserved and novel miRNA families. Secondly, the most abundant size for miRNAs was 21-24 nucleotides, except for HbmiR169 (27 nucleotides), HbmiR2911 (26 nucleotides), HbmiR482 (25 nucleotides), HbmiR472 (25 nucleotides), HbmiRn8 (25 nucleotides) and HbmiRn9 (27 nucleotides). Thirdly, the putative novel miRNAs occurred at low levels with fewer than 500 reads, and at extremely low levels for HbmiRn1 and HbmiRn2, with only three sequences (Figure 5).

Detection of single-nucleotide mutations in miRNA sequences

Identification of precursor transcripts for conserved and putative novel miRNA families from the hevea brasiliensistranscriptome

Conserved and putative novel miRNA precursors were sought in the Hevea clone PB260 transcript sequence databases obtained from various organs using the LeARN pipeline [51]. Eight conserved miRNA families that mapped against RNA sequences displayed a stem-loop structure (Table 4 and Additional Table 1). A single miRNA precursor transcript was found for the HbmiR156, HbmiR159, HbmiR396, HbmiR476 and HbmiR2910 families. Two miRNA precursor transcripts were identified for the HbmiR166, HbmiR319 and HbmiR408 families. In addition, the LeARN pipeline predicted 10 putative novel miRNAs with their unique precursors (Table 5 and Additional Table 2). Mapping the sequences of the small RNA dataset on the precursor sequences led to the identification of miRNA*for all eight conserved families (HbmiR156, HbmiR159, HbmiR166, HbmiR319, HbmiR396, HbmiR408, HbmiR476 and HbmiR2910) and for the five putatively new families HbmiRn3, HbmiRn4, HbmiRn8, HbmiRn9 and HbmiRn10. Stem-loop reverse transcription polymerase chain reactions (RT-PCR) were successfully performed on seven precursors for the conserved families and on nine for the putatively new families (Figure 6).

New miR name	Accession Name	Mature sequence	Mature length	miRNA* (Accession No.)	Precursor length	EST
HbmiRn1	acc_584359	CAGGAACUGGUAUCAACCCAGC	22	nd	130	CL1Contig1853_r1
HbmiRn2	acc_495109	UAUUGUAGAAAUUUUCAGGAUC	22	nd	116	CL1Contig490_r1
HbmiRn3	acc_185377	UAAUGGGCUCUGCAUAGAUGG	21	Yes (3)	108	CL27Contig1_r1
HbmiRn4	acc_108425	UUGCAUAUCUCAGGAGCUUCA	21	Yes (4)	214	CL19922Contig1_r1
HbmiRn5	acc_644211	AAACGGCUACCACAUCCAA	19	nd	87	FYGXE6I01B755W_r1
HbmiRn6	acc_501419	UAGGAUGUAGAAGAGCAUAA	20	nd	90	F5VNCTM02H9XL3_r1
HbmiRn7	acc_35715	UAGUUUGUUUGAUGGUAUC	19	nd	113	FYGXE6I01EKFYA_r1
HbmiRn8	acc_246127	GAUUGACAGACUGAGAGCUC	20	Yes (16)	93	FYGXE6I01B7GPL_r1
HbmiRn9	acc_103082	UUUAUGAAAGACGAACAACUG	21	Yes (6)	87	CL1Contig2798_r1

The presence of the miRNA* sequence in the small RNA dataset is indicated. Pictures of the stem-loop structures are presented in Additional Table 1.

MiRNA name	Accession Name	Mature sequence	Mature length	miRNA* (Accession No.)	Precursor length	EST
Hbmir156	acc_462135	UUGACAGAAGAUAGAGAGC	19	Yes (10)	119	FYGXE6I01DJH4X_r1
Hbmir159	acc_17403	UUUGGAUUGAAGGGAGCUCUA	21	Yes (5)	221	GETTCZF01AWUSL_r1
Hbmir166	acc_28021	UCGGACCAGGCUUCAUUCC	19	Yes (14)	127	CL2671Contig1_r1
	acc_150486	UCGGACCAGGCUUCAUUCCCCC	22	Yes (14)	120	CL1Contig15359_r1
Hbmir319	acc_19786	UUGGACUGAAGGGAGCUCCCU	21	Yes (7)	221	FRZQES201D34W3_r1
Hbmir396	acc_112787	UUCCACAGCUUUCUUGAACUG	21	Yes (15)	154	CL1Contig13811_r1
Hbmir408	acc_393516	AAGACUGGGAACAGGCAGAGCA	22	Yes (294)	128	FYGXE6I01CVDV9_r1
	acc_39581	ACUGGGAACAGGCAGAGCAUGA	22	Yes (294)	120	CL3908contig1_r1
Hbmir476	acc_508145	UAAUCCUUCUUUGCAAAGUC	20	Yes (1)	126	CL1Contig11471_r1
Hbmir2910	acc_251816	GAGCGAUUUGUCUGGUUAAUC	21	Yes (5)	126	FYGXE6I01EDXGX_r1

The presence of the miRNA* sequence in the small RNA dataset is indicated (nd: not detected). Pictures of the stem-loop structures are presented in Additional Table 2.

Comparison of the gene ontology classification of predicted targets for conserved and putative novel miRNA families

Computational prediction of miRNA targets led to the identification of 1,083 sequences for the 48 conserved families and 705 for the 10 putative novel miRNA families. Their involvement in biological processes and their molecular functions were attributed using Gene Ontology (Figure 7A, B, C and 7D). Most of the GO terms were represented in the same proportions for predicted targets from conserved and putative novel miRNA families. Target genes of the putative novel miRNA families all had GO terms of conserved miRNA families for biological processes and molecular functions. In biological process terms, the putative novel miRNA families had two additional terms (growth, cell wall organization or biogenesis). The number of miRNA target genes was smaller for putative novel than for conserved miRNAs for some GO terms: cellular process (25% as opposed to 27%), metabolic process (22% as opposed to 27%) and biological regulation (10% as opposed to 12%) (Figures 7A and 7C). This number was larger for putative novel miRNA targets compared to conserved targets for cellular component organization (5% as opposed to 3%), localization (7% as opposed to 4%) and response to stimulus (10% as opposed to 8%) (Figures 7A, B, C and 7D). As regards the GO terms for molecular function, although the distribution for antioxidant activity was the same for the 2 classes of miRNAs (1%), a decrease in the proportion of transcription regulator activity (4% as opposed to 7%) and an increase in structural molecule activity (4% as opposed to 2%) and transporter activity (7% as opposed to 4%) were observed for the putative novel miRNA families (Figures 7B and 7D).

Comparison of predicted target genes for conserved and putative novel miRNA families involved in "response to stimulus", and in "antioxidant" and "transcription regulation" activities

Twelve members of the AP2/ERF domain-containing transcription factor family were putatively targeted by HbmiR156, HbmiR159, HbmiR172, HbmiR393, HbmiR395, HbmiR396, and HbmiR408. Two HD-ZIP III domain-containing protein members (HbmiR165/166, HbmiR172), a transcription factor with a bZIP domain (HbmiRn4), three NAC domain-containing proteins (HbmiR164, HbmiR165 and HbmiRn7), three Squamosa promoter binding proteins (HbmiR156/157), four WRKY transcription factors (HbmiR156/157, HbmiR390, HbmiR444), four Auxin Response Factor (ARF); HbmiR160 and HbmiRn5) and five GRAS family transcription factors (HbmiR170, HbmiR171, HbmiR396 and HbmiRn5) were also potentially targeted.

Sequences encoding protein involved in protein degradation, such as cysteine protease (HbmiR396) and ubiquitin protein ligase (HbmiRn6 and HbmiRn7), were identified for genes classified in the "response to stimulus" category.

Three sequences related to ABA signalling were also identified: an ABA-insensitive protein (HbmiR166, and HbmiR2910), and an HVA22-like protein (HbmiR398). In addition, one sequence encoding CTR1, a negative regulator of the ethylene signalling pathway, was targeted by HbmiRn5. MiRNA biogenesis was also highlighted by the presence of a cap binding protein predicted to be targeted by HbmiR167. Interestingly, two enzymes involved in natural rubber biosynthesis, HMG-CoA reductase and 4-hydroxyl-3-methylbut-2-enyl diphosphate reductase, were predicted to be targeted by HbmiR159 and HbmiR444, respectively. For the "antioxidant activity" category, several ROS-scavenging enzymes were identified. The chloroplastic and peroxisomal isoforms of CuZnSOD were putatively targeted by HbmiR398 and HbmiR159, respectively. Another enzyme involved in glutathione biosynthesis, gamma glutamyl cysteine ligase (GCL), was predicted to be targeted by HbmiRn5. Rboh (or NADPH oxidase) is an enzyme involved in ROS production, which is likely to be targeted by two miRNAs (HbmiR2914 and HbmiR476). In addition, targets that are subjected to or are involved in the redox status of proteins, such as genes encoding proteins involved in lignin synthesis, cynnamoyl CoA reductase, lignin forming anionic peroxidase precursor and 4-coumarate:CoA ligase, were also found to be targeted by HbmiR159, HbmiR408 and HbmiR398, respectively.

Size distribution of small RNAs in Heveamay differ from others species

Next-generation high-throughput sequencing techniques are powerful tools for miRNA identification. The sequencing and analysis of the Hevea small RNA library generated more than 2 million sequences. In their vast majority (86%), the sequenced small RNAs were only represented by one or two reads. Similarly, 65% of unique miRNA sequences were found in Arabidopsis by 454 sequencing techniques [52]. Interestingly, we found that the size distribution of Hevea small RNAs differed from other species. The most abundant small RNAs had 17, 19 and 27 nucleotides in Hevea whereas they have 21 and 24 nucleotides in other species. Indeed, in Arabidopsis thaliana [52], Arachis hypogea [45] and Oryza sativa [53], the fraction of 24 nucleotides has proved to be predominant, accounting for 60%, 45% and 36% of unique small RNAs, respectively. The predominant size has been found to be 21 nucleotides in Pinus contorta [53] and Taxus chiniensis [54]. In Hevea, the high percentage of 17- and 19-nucleotide small RNAs resulted from the increased amount of RNA degradation products and transposable elements in response to abiotic treatments according to the mapping results against the Arabidopsis genome (data not shown), as previously observed in Brachypodium [55]. Indeed, a change in small RNA size distribution has been reported in Brachypodium in response to cold treatment. Under normal growth conditions, the 21- and 24-nucleotide classes were the most abundant small RNA families, whereas it was the 19-nucleotide class after cold treatment [55]. Our dataset was generated from pooled small RNAs of various tissues subjected to several types of abiotic stress. The size distribution may have been differentially affected by each type of abiotic stress and may explain the difference in size distribution in Hevea.

Long conserved and novel miRNAs of 23-27 nucleotides consisted of short miRNA sequences with an extension. This lmiRNA structure suggested that they were generated from the same MIR genes. Our observation suggests dual functions for 14 canonical MIR genes producing 20-22 nt and 23-27 nt lmiRNA in Hevea, namely HbMIR169, HbMIR2911, HbMIR482, HbMIR472, HbMIR408 and HbMIRn3-10. Across the plant kingdom, lmiRNA were previously classed as heterochromatic small interfering RNAs (hc-siRNA) [28]. In Arabidopsis, the hc-siRNA derived from the miRNA-generating site is dependent on DCL3, RDR2 and PolIV associated with AGO4 and direct DNA methylation for some target sites [28]. As regards the higher number of reads for the 24 nucleotide compared with the 20-22 nucleotide miRNAs, the 14 miRNA families producing 20-22 nt and 23-27 nt miRNAs might mostly regulate their targets through direct DNA methylation in response to abiotic stress. This may confirm the results of a recent study showing the existence of site-specific abiotic stress-induced DNA methylation for five genes in Hevea [56].

Ancient and recently differentiated Hevea brasiliensismiRNAs provide additional information for evolutionary studies

The 48 Hevea brasiliensis conserved miRNA families were more or less similar to those of other species. In Hevea, five miRNA families (HbmiR319, HbmiR156/157 and HbmiR165/166) have been found in more than forty plant species [57], and nine correspond to ancient miRNA already present in the common ancestor of Embryophytes (HbmiR156, HbmiR159/319, HbmiR160, HbmiR165/166, HbmiR171, HbmiR408 and HbmiR395), two in the common ancestor of Tracheophytes (HbmiR397 and HbmiR398) and nine in the common ancestor of Spermaphytes (HbmiR162, HbmiR164, HbmiR168, HbmiR169, HbmiR172, HbmiR393, HbmiR394, HbmiR399 and HbmiR827) [58]. Interestingly, the so-called monocotyledon-specific HbmiR444 family was identified in this study [58]. This suggests that HbmiR444 appeared before the speciation between monocotyledons and dicotyledons. Three woody species-specific miRNA families were also found in Hevea brasiliensis: HbmiR473 (Citrus sinensis and Populus trichocarpa), HbmiR476 (Populus trichocarpa) and HbmiR479 (Populus trichocarpa, Vitis vinifera, Gossypium hirsutum and Citrus sinensis) and were probably recruited independently by woody species during their evolution. The 10 novel miRNA families identified in this work are likely to be more recent, after Hevea diversification, as they have not been identified up to now in other euphorbiaceous plants. Furthermore, this information complements previous studies on euphorbiaceous species [49] in which only 25 conserved miRNAs were identified.

Our data could be used to study the evolution of MIR genes. We showed that 70% of the putatively new MIR genes were able to produce lmiRNA of 23-27 nucleotides, compared to 29% for conserved MIR genes. This corresponds to what has been observed in Arabidopsis, where it was previously suggested that young MIR genes consistently produced siRNAs, while ancient MIR genes produced predominantly canonical miRNAs, showing that the evolution of MIR genes is associated with changes in the use of DCL, resulting in specific miRNA size classes [59]. In addition, several so-called plant species-specific miRNAs were also found in Hevea, such as HbmiR1310 initially found only in Pinus taeda [53], HbmiR2910, HbmiR2914, HbmiR2915 and HbmiR2916 found in Populus euphratica [60] and HbmiR476 in Populus trichocarpa [61]. Our data suggest that those miRNAs should no longer be considered as species-specific miRNA.

Our study also revealed some conserved and divergent mechanisms at post-transcriptional level. Indeed, four families of miRNA were subjected to RNA editing events at the sites previously observed in Arabidopsis. However, for the thirteen families left, RNA editing was not observed under our experimental conditions. Under our experimental conditions, the -1 + UU modification was not observed in Hevea for the nine families selected from Ebhardt et al. [50]. Once the Hevea genome sequence is available, it will be possible to conduct an exhaustive analysis of the discrepancies between members of a multigenic family and RNA editing events.

Predicted target/miRNA pairs conserved between plant species and putative new ones in Hevea brasiliensis

In comparison with other plants, several regulation mechanisms seem to be conserved in Hevea brasiliensis, such as the regulation of HD-ZIP III protein by HbmiR166 [62], Squamosa promoter binding protein by HbmiR156 [63], NAC domain protein by HbmiR164 [64], APETALA2-like protein by HbmiR172 [65–67], CCAAT-binding transcription factor by HbmiR169 and cysteine proteinase by HbmiR396 [44]. These ancient miRNAs regulate ancestral transcription factors that coordinate highly conserved functions, such as organ polarity and separation, cell division, or hormonal control [68, 69]. We predicted potentially new miRNA/target pairs in Hevea. New target/miRNA pairs were found for transcription factors, such as, for example, the NAC domain protein by HbmiR165 and HbmiRn7, the GRAS domain protein by HbmiR170, HbmiR396 and HbmiRn5, WRKY by HbmiR396 and HbmiR444, and the APETALA2-like transcription factor by HbmiR156, HbmiR159, HbmiR393, HbmiR396 and HbmiR408. (Duan et al., submitted). In addition, we also detected new pairs involved in protein degradation (ubiquitin protein ligase by HbmiRn6 and HbmiRn7), in response to hormones (ABA-insensitive by HbmiR166, ABA responsive element binding protein by HbmiR2910, and HVA22-like protein by HbmiR398, CTR by HbmiRn5) and in the regulation of miRNA (cap binding protein by HbmiR167). Further analyses are needed to validate the inhibition of these targets at transcriptional or translational level and to unravel specific regulation in Hevea brasiliensis regarding the response to stress in latex cells, and latex production.

Conserved and divergent regulation of redox homeostasis in Hevea brasiliensis

Several enzymes involved in ROS production or detoxification were found in our analysis. Two ROS-producing enzymes, Rboh, are likely to be targeted by HbmiR2914 and HbmiR476. They are key regulators in the ROS signalling network as they integrate different signal transduction pathways, such as calcium, protein phosphorylation and lipid signalling with ROS production [6]. Careful attention was also paid to targets classified in the GO term "antioxidant activity". For conserved miRNAs, two CuZnSOD chloroplastic and peroxisomal isoforms were predicted to be targeted by HbmiR398 and HbmiR159 respectively. Surprisingly, the CuZnSOD cytosolic isoform was not predicted by the LeARN pipeline and, even manually, no correct alignment was possible with HbmiR398. So far, validation by 5' RACE PCR has failed under our experimental conditions. By contrast, previous studies in Arabidopsis showed that HbmiR398 targeted both chloroplastic and cytosolic isoforms [34].

Some miRNA families have been reported to be H₂O₂-responsive in rice, such as OsmiR169, OsmiR397, OsmiR528, OsmiR827, OsmiR319 and OsmiR408 [70]. In rice, OsmiR397 targets laccase genes involved in lignin biosynthesis and those targets were also found for HbmiR397 in Hevea.

Production of in vitro plantlets

Embryogenic friable callus from line CI05538 of clone PB 260 was maintained on MM medium in the culture room at 27°C, with 60% humidity in the dark [71]. In vitro plants were regenerated according to the procedure described in [72]. Briefly, somatic embryogenesis was initiated by sub-culturing 1 gram of callus in 250 mL flasks containing 50 mL of a semi-solid Induction Medium. Pro-embryos were then developed in a temporary immersion system (RITA^®, CIRAD, France) with 1 min of immersion per day in the liquid development media DEV1 and DEV2. For plant regeneration, only well-shaped mature embryos were collected after 8 weeks of culturing [73]. They were transferred to glass tubes on a semi-solid DEV3 medium for one month under a light intensity of 60 μmol m^-2 s^-1 and a 12 h light/dark photoperiod for embryo conversion into plants. Plantlets were then acclimatized in the greenhouse at 28°C with 60% relative humidity.

Plant material cultivation and treatments

Budded plants and two-year-old in vitro plants of clone PB260 were grown in the greenhouse at 27°C, with 50% humidity under natural light conditions. Budded plants and two-year-old in vitro plants of clone PB260 were exposed to cold at 4°C for 8 h, light (PAR = 1000-1500 μmol/m²/s for 4 h), drought (4 weeks for young budded plants and 10 days for two-year-old in vitro plants without watering), flooding (covered with water for 4 days), salinity (300 mM NaCl for 4 days) and wounding for 4 h (Table 2). Hormonal treatments were applied in hermetically sealed boxes, either with 5 ppm of ethylene for 4 h or 0.3 μM of methyl jasmonate for 4 h as described in [74, 75].

Isolation of small RNAs

For all treatments, leaves, bark and roots were sampled in liquid nitrogen and then frozen at -80°C pending small RNA isolation. Small RNAs were purified with the mirPremier microRNA Isolation Kit (Sigma-Aldrich, St. Louis, MO, USA), following the manufacturer's instructions for plant tissues. Briefly, 100 mg of the plant tissue powder was mixed with a lysis mix that released small RNA and, at the same time, inactivated ribonucleases and interfering secondary metabolites that might exist in plant tissues. Large RNA and genomic DNA remained insoluble and were removed from the lysate along with other cellular debris by brief centrifugation. Small RNA was then captured on a silica binding column in the presence of alcohol. Residual impurities were removed by wash solutions, and purified small RNA was eluted in RNase-free water.

Deep sequencing of small RNAs

Small RNA samples from callus and from control and treated plants were pooled. A Hevea small RNA library was constructed by the SKULTECH Company (France). Briefly, small RNAs were separated on a denaturing 15% polyacrylamide gel. The gel slices containing RNA with a size of about 15-50 nucleotides were excised and the RNA was eluted. Small RNAs were ligated to 5' and 3' adapters using T4 RNA ligase. Products from the second ligation were gel purified then reverse transcribed and amplified using PCR to produce cDNA. The PCR was performed with two primers that annealed to the ends of the adapters. PCR products were separated on 15% polyacrylamide gel with ethidium bromide staining. The gel slices containing DNA with a size of about 92 nucleotides were excised and the DNA was eluted. The resulting cDNA was sequenced using Solexa/Illumina technology (Illumina Inc., San Diego, CA, USA). Raw reads were trimmed using Trim_Adaptor application (non Open source application) with the minimum clip length fixed at 6 nucleotides.

Bioinformatic analyses

This work was supported by the high performance cluster of the SouthGreen Bioinformatics platform (southgreen.cirad.fr - UMR AGAP - CIRAD) comprising 26 compute nodes (208 Nehalem cores) and a 16 to storage capacity connected through a low latency Infiniband network. Sun Grid Engine was used as the workload manager to schedule jobs on the different compute nodes.

Identification of the conserved heveamiRNAs

Adapters and low-quality sequences were removed from the raw reads obtained from Illumina sequencing using cross-match software http://www.phrap.org/phredphrpconsed.html; the non-redundant dataset was produced using WU nrdb software. The unique small RNA sequences were then compared with BlastN [76] against the Arabidopsis genome (Phytozome V6) and the PMRD database [38]http://bioinformatics.cau.edu.cn/PMRD in order to classify the small RNAs and identify mature miRNAs in Hevea already known in several species.

Identification of HeveamiRNA precursors

Conserved miRNAs and putative novel miRNA precursors were identified in the PB260 transcriptome http://bassigny/cgi-bin/esttik_dev/quick_search.cgi using the LeARN pipeline http://symbiose.toulouse.inra.fr/LeARN/release_1;5.0[50], a platform for detecting, clustering and annotating non-coding RNAs. Novel miRNAs were separated into 4 classes according to their position on the stem-loop structure. Class 1 precursors only produced small RNAs corresponding to the predicted miR and miR* (i.e., small RNAs in a 25-bp region around the miR:miR* predicted region to be considered with two nucleotide overhangs linked to DCL action); class 2 comprised hairpins containing both predicted miR and miR* and additional small RNAs, of lower abundance, outside this region; class 3 involved precursors with reads on the defining small RNA region only (3a, at least two reads; 3b, only one read); class 4 consisted of precursors containing the defining small RNA plus reads of lower abundance outside the predicted miR:miR* duplex. Finally, class 5 contained hairpins that generated small RNAs of greater abundance than the defining RNA outside the predicted miR:miR* region [77]. Each miRNA precursor was checked manually following the rules described in [45]. Specifically, we considered dominant, mature sequences residing in the stem region of the stem-loop structure and ranging between 20-22 nt with a maximum free-folding energy of -25 kcal mol^-1. A maximum of six unpaired nucleotides between the miRNA and miRNA* was allowed. The distance between the miRNA and miRNA* was fixed between 5 and 240 nucleotides.

Identification of putative HeveamiRNA targets

Computational identification of miRNA targets was performed with the Miranda toolbox included in the pipeline with the default parameters (gap_value = 2, mm_value = 1, gu_value = 0.5, score_threshold = 3, min_length_alignment = 20 and no_mismatch_positions = 10;11) for conserved miRNA and secondly checked with psRNAtarget http://plantgrn.noble.org/psRNATarget/(length_alignment = 18), [78]. For the putative novel miRNAs, psRNAtarget server (length_alignment = 18) and Miranda (gap_value = 2, mm_value = 1, gu_value = 0.5, score_threshold = 3, min_length_alignment = 18 and no_mismatch_positions = 10; 11) were used to scan the Hevea EST databases (Duan et al., to be submitted, http://bassigny/cgi-bin/esttik_dev/quick_search.cgi). A first assembly set was generated from reads of leaves, bark, latex, embryogenic tissues, and roots separately to create tissue-specific transcript databases with sequence names CLxxcontigxx. Then, reads from all tissues were used to generate one transcript sequence database for the Hevea clone PB260 with the sequence name hevea_454_xxx.

Gene Ontology analyses were performed with the Blast2GO program [79, 80]. In detail, sequence analyses, such as alignments and BLAST, were performed with Geneious software [81].

Total RNA extraction

Total RNAs from leaves, bark and roots subjected to abiotic stress treatments were extracted using the caesium chloride cushion method adapted from Sambrook [82] by Duan et al. [75]. One gram of fresh matter was ground and transferred to a tube containing 30 ml of extraction buffer consisting of 4 M guanidium isothiocyanate, 1% sarcosine, 1% polyvinylpyrrolidone and 1% ß-mercapto-ethanol. After homogenization, the tubes were kept on ice and then centrifuged at 10,000 g at 4°C for 30 minutes. The supernatant was transferred to a new tube containing 8 ml of 5.7 M CsCl. Ultracentrifugation in a swinging bucket was carried out at 89,705 g at 20°C for 20 hours. The supernatant and caesium cushion were discarded, whilst the RNA pellet was washed with 70% ethanol. After 30 minutes of air drying, the pellet was dissolved in 200 μl of sterile water. RNAs were conserved at -80°C.

Experimental validation of miRNA targets

The GeneRacer kit (Invitrogen, Carlsbad, CA, USA) was used to ligate a 5'RNA adapter to pooled RNAs prior to reverse transcription. The resulting cDNAs were used as a template for PCR amplification using one Gene-Specific Primer (GSP) designed downstream of the predicted miRNA:target binding site and one GeneRacer 5' forward primer. The primers used were HbCuZnSODchl R: GGGTAACCAGCAAATGCAAGCAGC, HbARF R: CTTGTGTTGAGTTGCAGCGCG and HbSquamosa R: TGGAGCCTAATTGGCTTCCTTCTGC. The first PCR product was then used for nest amplification performed with a nested GSP primer and a nested GeneRacer 5' forward primer. The nested primers used were HbCuZnSODchl nested R: GCAGGGAACAATGGCTGCC, HbARF nested R: ATGTTCAGCTCAGTGGACACGG, and HbSquamosa nested R: TTGTTTGGCACCACGCTTGTGGAAG. The PCR products were gel purified, cloned into PCR2.1 TOPO TA vector (Invitrogen, Carlsbad, CA, USA) and sequenced.

Complementary DNA (cDNA) synthesis and stem-loop RT-PCR