Post-transcriptional regulation of several biological processes involved in latex production in Hevea brasiliensis

Background Small RNAs modulate plant gene expression at both the transcriptional and post-transcriptional level, mostly through the induction of either targeted DNA methylation or transcript cleavage, respectively. Small RNA networks are involved in specific plant developmental processes, in signaling pathways triggered by various abiotic stresses and in interactions between the plant and viral and non-viral pathogens. They are also involved in silencing maintenance of transposable elements and endogenous viral elements. Alteration in small RNA production in response to various environmental stresses can affect all the above-mentioned processes. In rubber trees, changes observed in small RNA populations in response to trees affected by tapping panel dryness, in comparison to healthy ones, suggest a shift from a transcriptional to a post-transcriptional regulatory pathway. This is the first attempt to characterise small RNAs involved in post-transcriptional silencing and their target transcripts in Hevea. Methods Genes producing microRNAs (MIR genes) and loci producing trans-activated small interfering RNA (ta-siRNA) were identified in the clone PB 260 re-sequenced genome. Degradome libraries were constructed with a pool of total RNA from six different Hevea tissues in stressed and non-stressed plants. The analysis of cleaved RNA data, associated with genomics and transcriptomics data, led to the identification of transcripts that are affected by 20–22 nt small RNA-mediated post-transcriptional regulation. A detailed analysis was carried out on gene families related to latex production and in response to growth regulators. Results Compared to other tissues, latex cells had a higher proportion of transcript cleavage activity mediated by miRNAs and ta-siRNAs. Post-transcriptional regulation was also observed at each step of the natural rubber biosynthesis pathway. Among the genes involved in the miRNA biogenesis pathway, our analyses showed that all of them are expressed in latex. Using phylogenetic analyses, we show that both the Argonaute and Dicer-like gene families recently underwent expansion. Overall, our study underlines the fact that important biological pathways, including hormonal signalling and rubber biosynthesis, are subject to post-transcriptional silencing in laticifers.

The proteins involved in biogenesis are highly diversified among plant species, which may explain the emergence of new functions and regulatory pathways (D'Ario, Griffiths-Jones & Kim, 2017). The impact of sRNA-mediated regulations on plant phenotypes has been illustrated through several examples, including association between increased apple fruit size and MIR172 gene deregulation (Yao et al., 2016), and alleviation of siRNA-dependent repression of FAD2 expression resulting in increased oil content in olive fruits (Unver et al., 2017). In annual plant species, size classes and respective abundances of sRNA transcriptomes may be affected by a variety of environmental stresses (Ferdous, Hussain & Shi, 2015;Khraiwesh, Zhu & Zhu, 2012;Liu et al., 2017;Yang et al., 2017;Zhao et al., 2016). Similar observations have been made in cultivated perennial species (Beule et al., 2015;Li et al., 2017).
Previous work by our team on Hevea brasiliensis revealed alterations in the sRNA transcriptome of latex cells in response to environmental cues (Gébelin et al., 2013b). This tropical perennial crop is used for the production of latex containing natural rubber (cis-1,4 polyisoprene) (Compagnon et al., 1986). Severe environmental constraints and stresses associated with latex harvesting practices, such as wounding through tapping or ethephon stimulation, cause physiological disorders leading to tapping panel dryness (TPD). TPD is characterized by a halt in latex flow and in situ latex coagulation (for a review Zhang, Leclercq & Montoro, 2017). In a previous study comparing healthy and TPD-affected trees, we observed a shift in the respective proportions of different size classes of sRNA in latex (Gébelin et al., 2012;Pramoolkit et al., 2014), the 24 nt size class predominating in healthy trees, and the 21 nt size class in TPD-affected ones. This observation led to the need to annotate small RNAs to decipher their biological activity, in particular that of 20-22 nt (miRNA and ta-siRNA). Although not complete, previous studies characterised endogenous miRNAs in Hevea, thereby enabling the identification of 68 and 16 species-specific families, thanks to the availability of a reference transcriptome from clone PB 260 (Gébelin et al., 2012;Gébelin et al., 2013b;Kanjanawattanawong et al., 2014;Lertpanyasampatha et al., 2012). The aim of the present study was to complete the annotation of Hevea microRNAs and ta-siRNAs thanks to the availability of genomic sequences from clone PB 260 (Zhang et al., 2019). Previous computational prediction of miRNA targets led to the identification of 1,788 sequences (Gébelin et al., 2012;Pramoolkit et al., 2014). However, among these target miRNAs, few have been experimentally validated so far (Gébelin et al., 2012;Pramoolkit et al., 2014). Therefore, the real impact of rubber tree miRNAs on transcript regulation remains to be fully assessed.
Unlike animal miRNAs, the main mode of action of plant miRNAs is through transcript cleavage by endonucleolytic slicing (Arribas-Hernandez, Kielpinski & Brodersen, 2016), thereby generating degradation products that are collectively known as ''degradome''. To improve our knowledge of post-transcriptional regulatory mechanisms in the clone PB 260 rubber tree, we generated and analysed six degradome libraries originating from six different tissues to detect miRNA-and ta-siRNA-mediated cleaved targets. Their 3 ends were then sequenced so that all sRNA-directed cleavage sites could be identified simultaneously (German et al., 2009). The targets were then functionally annotated.
Here we discuss our results with a focus on the laticiferous tissue and the natural rubber biosynthetic pathway. Our results suggest that laticiferous tissues have the highest level of post-transcriptional regulation by mRNA cleavage, and that the natural rubber (NR) biosynthetic pathway is under strong post-transcriptional silencing. These results could be used to achieve functional validation in rubber trees, under the physiological stresses associated with latex harvesting.

Plant material and RNA samples
The plant material used for this study comes from multiple origins but all from Hevea clone PB 260. Samples were collected during in-vitro culture, from in-vitro plants, and from juvenile budded plants grown in the greenhouse at CIRAD (Montpellier). Other samples were collected from mature trees grown in the field at the Sembawa Rubber Research Center in Indonesia and the Thailand Chachoengsao Rubber Research Centre.
Total leaf, bark and root RNAs were generated from juvenile plants subject to cold, high light, drought, flooding, salinity and wounding stress by Gébelin et al. (2012). Total leaf and bark RNAs from juvenile plant treated with hormones (ethylene, methyl jasmonate) were generated by Duan et al. (2010) and Kuswanhadi et al. (2010). Total RNAs from latex and bark of healthy and TPD affected mature trees were prepared by Gébelin et al. (2013b), and total RNAs from male and female flowers by Piyatrakul et al. (2012). Total RNAs from somatic embryos were generated by Piyatrakul et al. (2012). All these total RNAs from in-vitro culture, juvenile and mature plants were used to build specific leaf, bark, root, embryo, flower and latex degradome libraries.
Differential expression analysis of RNAseq data in latex from mature trees were performed by Montoro et al. (2018). Pairwise comparisons were performed between untreated and ethephon-treated healthy trees (WH vs EH) and between healthy and TPD ethephon-treated trees (EH vs ET).

Nuclear genome resequencing of clone PB 260
Leaf nuclear DNAs from a budded clone PB 260 were sequenced by GATC (https: //www.eurofinsgenomics.eu/) using Illumina pair-end sequencing (2 × 150 bp). About 280 million reads were generated, which correspond to a total of 84 Gb among which 62 Gb were assembled through alignment against the clone Reyan7-33-97 reference genome (Tang et al., 2016). The Bowtie tool (Langmead et al., 2009) was used to build a Burrows-Wheeler transform index to map small RNA libraries with specific parameters [bowti e − f − v 1 -p 10 − a − m 50 -best -strata (Langmead et al., 2009)]. The statistics used to resequence the genome of clone PB 260 are detailed in (Zhang et al., 2019). The resequencing data are available under project number PRJCA001333 in the GSA (Wang et al., 2017a) and BIG Data Center (BIG Data Center Members, 2018).

Detection and GO classification of cleaved target transcripts
Raw sequencing reads originating from degradome libraries were trimmed with Cutadapt (Martin, 2011). The CleaveLand pipeline (Addo-Quaye, Miller & Axtell, 2009) was used in Mode 1 (the mode that allows alignment of degradome data, creation of degradome density file, and new small RNA query/transcriptome alignment with GSTAr) in order to identify sRNA-generated transcript cleavage products (Fig. 2). The abundance of reads starting at the cleavage site is defined by the degradome category. Category 4 has just one read at this position. Other categories have more than one read. In category 3, the number of reads is below or equal to the average depth of coverage (all positions that have at least one read) on the transcript. In category 2, the number of reads is equal to the average depth of coverage on the transcript. In category 1, the number of reads is equal to the maximum of the average depth of coverage on the transcript when there is more than one position at the maximum value. In category 0, the number of reads is equal to the maximum of the average depth of coverage on the transcript when there is just one position at the maximum value. The outputs are presented in Tables S6, S9-S13 for miRNAs and S15, S17-21 for ta-siRNA).

Post-transcriptional regulation analysis
Functional annotation of the identified targets was performed using InterProScan (v5.18-57.0) in STANDALONE mode. GO classification was then extracted from the resulting tabular output (Ashburner et al., 2000).

Prediction of trans-acting siRNA production sites
Small RNAs were aligned to the PB 260 genome assembly, and small RNAs that did not match were discarded. As input sRNA files, we used trimmed reads from latex cells from healthy and TPD-affected trees and juvenile plants, with a minimum sRNA abundance of two, a p-value cut-off of 0.0001 and a length of 21 nt. Using the small RNA workbench from the University of East Anglia (http://srna-workbench.cmp.uea.ac.uk/; version 4.3.1 Alpha) (Stocks et al., 2012), the 21-mers ta-siRNA loci were identified using the TA-SI Prediction tool (p-value cut-off <0.0001). The algorithm for ta-siRNA enabled us to calculate phasing probability based on hypergeometric distribution (Chen, Li & Wu, 2007). Coding potential was estimated at each of the identified sites using the Coding Potential Calculator (CPC) (http://cpc.gao-lab.org/) (Fig. 3).
Detection of ta-siRNA loci (88 with sRNA from latex healthy trees, 0 from latex of TPD-affected trees and 12 from in juvenile plants  Annotation of the miRNA biogenesis and regulatory pathway and phylogenetic analysis

RESULTS
The analyses of miRNAs and of degradomes from data processing to functional annotation by gene ontology are depicted in Fig. 1. The first step was to annotate the stem-loop structures and the resulting miRNAs. The second step was to analyse the distribution of predicted miRNAs according to size by comparing three conditions in juvenile plants, and latex from healthy and TPD-affected trees. The third step was to analyse their cleavage activities either using the tissue-specific transcriptomes or a list of well-characterised rubber tree genes. Figure 3 shows the process of annotation of ta-siRNAs and their cleavage activities.

Annotation of MIR genes forming RNA hairpin structures and producing small RNAs
Hevea miRNAs and MIR genes (hairpin-forming sequences) were annotated based on high-throughput sequencing data obtained from the five sRNA libraries and from the resequenced genome of clone PB 260. The genome was mapped on the reference genome of clone Reyan7-33-97 (Table S1). An atlas of 1,042 MIR genes, associated with miRNA production, was generated ( Fig. 1, Tables S2 and S3), and annotated. Using both the MirBase and PNRD databases, only, respectively, 29 and 36 MIR sequences matched known MIR sequences from different plant species (Table S4). This study complements our previous analyses (Gébelin et al., 2012;Gébelin et al., 2013b), which showed only eight conserved MIR genes annotated from juvenile plants, and four conserved and 11 new MIR genes from latex of healthy and TPD-affected trees.
To perform a degradome analysis, Ahmed and co-workers (Ahmed et al., 2014) recommend taking miRNA and isomiRNA (miRNA variants) pools into account. A Hevea sRNAome containing 1,839 redundant sequences was constructed with 20-22 nt sized sRNA sequences originating from MIR genes, with or without similarity to known miRNAs ( Fig. 1, Table S5). This file was used as input to compare miRNA size class distribution across sRNA libraries, and to explore the degradome libraries to detect cleaved target transcripts/miRNA pairs.

Determination of miRNA size class distribution and base preference
The distribution of unique miRNA accessions identified in small RNA libraries from juvenile stressed plants, latex from healthy and TPD-affected tree tissues is shown in Fig. 4. The total lengths ranged from 17 to 35 nucleotides. The distribution of the size classes showed a major peak at 21 nt, which accounts for 24.4 to 25.2% of the total number of predicted miRNAs (Fig. 4). Interestingly, this proportion was seen to be globally conserved in the different libraries tested (Fig. 1, Tables S2 and S3). The 5 nucleotide identity of a miRNA guides the loading into AGO proteins (Rogers & Chen, 2012) (Frank et al., 2012Zha, Xia & Yuan, 2012;Yu, Jia & Chen, 2017). AGO1 is known to sort miRNA with a 5 uracil, and AGO2, 4, 7 and 10 also showed slicing activity with no sorting priorities (Rogers & Chen, 2012). Hevea predicted miRNAs 21-nt in size were further analysed for their 5 base preference (Fig. 5). In the latex from healthy and TPD-affected trees, the 5 nucleotide  of predicted miRNA was mainly uracil, in contrast to the guanidine found in the latex from juvenile plants.

Detection of 136 cleaved target transcript/miRNAs pairs in the degradome at the transcriptome level
The sequencing of the six degradome libraries, from six distinct tissues (leaf, bark, root, flowers, latex and embryos), generated a total of 16 million single reads. The highest and lowest numbers of cleaned reads were found in latex (134,611 reads) and embryo (9,283 reads) libraries, respectively. The distribution of the number of reads per library after the cleaning process is shown in Table S1.
To validate the quality of the degradome libraries, we looked for the position of miRNAdirected cleavage sites for three known target transcripts that had previously been identified experimentally (Gébelin et al., 2012). We successfully detected SBP/CL2120Contig2 targeted by miR156/acc_480780, ARF/CL6582Contig1 by miR160/acc_370 and chloroplastic CuZnSOD/CL4308Contig 2 by miR398/acc_420 (Fig. 6). Degradome products corresponding to these known targets were detected in degradome libraries only in certain tissues, suggesting post-transcriptional regulation could be tissue specific (Table S6). For instance, cleavage products of the chloroplastic CuZnSOD were found in latex and root libraries. Those of ARF were identified in libraries from bark, and those of SBP in leaf and bark. Moreover, we checked for the statistical classification of the category of degradome, which is based on coverage at the cleavage site compared to the total length of the transcript. The three known miRNA-mediated cleaved targets belonged to the worst degradome, i.e., category 4, consisting of only 1 read at the cleavage site (Fig. 6). This result  Table 1, Table S7). By contrast, the embryo degradome yielded the lowest number of targets with only 4 (2.9%) detected. Only 26 target transcripts (19.1%) were cleaved by known miRNAs (16 different accessions). The miRNAs involved in mRNA slicing exhibited a tissue-specific activity pattern with very little overlap between tissues, as shown in the Venn diagram (Fig. 7).
All the cleaved targets identified belonged to degradome categories 3 and 4 (meaning only one or very few reads starting at the cleavage site), except for six mRNAs which were in degradome category 0 (meaning the majority of reads started at the cleavage site). Among them, functional domains were predicted in the putative peptides of only two of the targets found in latex, CL5497Contig1, which included a leucine-rich repeat domain (cysteinecontaining subtype), and CL8280Contig1, which included an RNA-binding/CRM domain. Gene ontology annotation strongly suggests that a wide range of biological processes are under post-transcriptional regulation in different rubber tree tissues, including sugar biosynthesis in leaves, ion transport in roots, protein imports and sugar metabolism in the bark, and RNA metabolic processes in latex (Fig. 1, Table 1, Table S8).

Detection of 173 additional miRNA-mediated cleaved transcripts in curated genes involved in laticifer differentiation and natural rubber production
The occurrence of post-transcriptional regulation of genes involved in important biological pathways of Hevea was then investigated in more detail. The gene families have already been described by Duan et al. (2013), Kuswanhadi et al. (2010), Pirrello et Carbohydrate metabolic process,(1->3)-beta-D-glucan biosynthetic process, glycolytic process, glutamine biosynthetic process, nitrogen compound metabolic process, signal transduction, photosystem II assembly, malate transport, photosynthesis, protein import  (Fig. 1, Table 2, Tables S9 to S13). In the present study, a high proportion of transcripts were found to be cleaved by miRNA, including genes with functions related to ethylene biosynthesis (23 genes out of 34 tested including With respect to TPD-associated genes (Venkatachalam, Thulaseedharan & Raghothama, 2007;Venkatachalam, Thulaseedharan & Raghothama, 2009), nine out of 12 genes showed evidence of post-transcriptional regulation in all tissues tested, with the exception of embryo. Among them, transcripts produced by HbTOM20, which encoded a translocase of the outer mitochondrial membrane, were cleaved by both miR396 and miR395. Moreover, a gene similar to HbTCTP, which has multiple functions and is associated with stress/hormone responses and TPD, may be cleaved by two newly identified miRNAs (Pyoung96216 and health13399; Table S13). The latter completes our previous study in which HbTCTP transcripts were degraded by two distinct sRNAs, miRf12236-akr and miR1023b-3p (Deng et al., 2016).
Unlike reports in the literature (Chen, 2004), we found no evidence of AP2/ERFs transcripts targeted by miR172 in our experimental conditions. Instead, these transcripts were targeted by other known miRNAs (miR156, miR167, miR9386 and miR396) and by 97 newly identified miRNAs. Examples of the good quality of the detection of miRNAmediated cleaved transcription factors (HbMYC, HbRAV-03, HbERF-VIIa8 and HbERF-VIIa5), with a degradome category of 0, in a T-plot generated by CleaveLand, can be seen in Fig. 8.

Identification of putative ta-siRNA-producing loci and their targets
In addition to miRNAs, another type of 21-nt small RNA may also target transcripts to induce their degradation. Trans-acting siRNAs (ta-siRNAs) are generated through the processing of a double-stranded RNA precursor, transcribed from a TAS gene targeted by a single 22-nt miRNA or by two 21-nt miRNAs (Bologna & Voinnet, 2014), through cleavage sites spaced at 21-nt intervals along its sequence (Komiya, 2017). Putative Hevea loci producing ta-siRNAs were identified by mapping sequenced 21-nt sRNAs on the assembled PB 260 genomic sequences (Fig. 3). This step led to the localisation of 88 non-redundant putative loci with the ta-siRNA originating from the latex of healthy mature trees, none from TPD-affected mature trees, and 12 from juvenile plants (Table  S14, Fig. 3). Annotation of the loci showed that two loci correspond to sequences coding for calcium-dependent protein kinase 1 (CDPK1) gene, and that 46 match the intron of farnesyl diphosphate synthase (FINT1), a gene containing a transposon (Zhang, Leclercq & Montoro, 2017). Finally, up to four ta-siRNA-producing loci were found in repeat-rich regions such as putative retrotransposons () and unannotated genomic scaffolds (0 to  Degradome analysis showed that the number of cleaved targets (15) was lower than miRNAs (Fig. 3, Table 3). Most ta-siRNAs cleaved fewer targets than miRNAs. Most of these targets were found in the leaf (6) whereas only one target was found in latex (Table  3 and Table S15). An important accumulation of cleavage products was observed for the CL7Contig17 transcript, which could point to highly efficient ta-siRNA-mediated processing in the latex of TPD-affected trees (accession TPD_latex_scaffold0858_100404).
No targets corresponding to ta-siRNAs derived from coding sequences were identified. Annotation of the ta-siRNA targets (Fig. 2, Table 3 and Table S16) showed that these targets were either involved in different biological processes, linked to those in which miRNA targets played a role, or that their mode of post-transcriptional regulation differed, depending on the nature of the tissue. For instance, most transcripts involved in sugar metabolism and transmembrane transport were regulated by ta-siRNAs in latex, while in bark, they were regulated by miRNA (Table 3, Tables S15 and S16). Moreover, compared to miRNA targets, more ta-siRNA targets were involved in rubber synthesis, hormone signalling and laticifer differentiation (HbREF8, HbHMGS1, JAZ6, Topless, MED25, CTR2, and HbMYB gene families; Tables S17 to S21). Overall, most of the targets cleaved by ta-siRNAs differed from those cleaved by miRNAs. A target cleaved by both types of sRNA required two distinct sites.

Identification and regulation of genes involved in microRNA biogenesis and regulatory functions in latex cells
MiRNA production requires MIR gene multi-step processing, which involves several multigene families. To annotate them in Hevea, 27 protein sequences involved in miRNA biogenesis and activity in Arabidopsis (Naqvi et al., 2012) were used as queries. Among them, 24 unique loci encoding putative orthologues were identified in the Hevea reference genome (Tang et al., 2016) using a sequence similarity search, including members of the AGO, DCL, HST, HEN1, SDN and SE families. Analysis of AGO and DCL gene families showed that each of them included four paralogous clades, in accordance with previous reports (Fang & Qi, 2016;Fukudome & Fukuhara, 2017  AGO5 orthologues and two putative AGO6 orthologues; Fig. 9A, Fig. S1), comparable in size to that of the recently duplicated Populus trichocarpa genome (15 genes including three AGO5 and one AGO6 ) (Lei et al., 2012). This means this gene family is bigger and more diversified than in the genomes of Ricinus communis (eight AGO genes) and Arabidopsis thaliana (10 genes). It was difficult to assign AGO2/3 (two genes found in Hevea: mRNA00000107 and mRNA00000478) and AGO4/8/9 (only one gene found in Hevea: mRNA00000007 ). By contrast, the DCL gene family of Hevea (five genes, including two putative DCL2 orthologues; Fig. 9B, Fig. S2) was in the same size range as that of both Ricinus communis and Arabidopsis thaliana (four genes each), while being smaller than that of Populus trichocarpa (six genes including two DCL2 orthologues). The transcripts of the putative Hevea DCL1 and AGO1 orthologues (mRNA00000029 and mRNA00000323, respectively; Fig. 10) were identified, and their expression response in response to ethephon stimulation and TPD occurrence was analysed (Montoro et al., 2018) (Table S22). DCL1 expression was significantly downregulated as a result of ethephon treatment (1,591 vs. 659 reads), but not according to TPD. A very high level of expression was observed for AGO1 across all conditions tested (>19,000 reads). While putative members of the Serrate, Hasty and SDN gene families were found to have variable expression levels depending on the family member, no evidence was found for regulation in either stressed or TPD conditions (Fig. 10, Table S22). In latex, putative orthologous transcripts were detected for each gene involved in either miRNA biogenesis or its regulation. This suggests that the components of the miRNA machinery are functional in this tissue. However, whether these different genes play the same roles in rubber as in model plant species remains to be determined. Indeed, Hevea AGO and DCL gene families likely expanded due to whole genome duplication (Deng et al., 2016).

DISCUSSION
The aim of this study was to identify MIR genes and their corresponding miRNAs in genomic sequences of Hevea clone PB 260. In the second step, the targets of these miRNAs were identified by analysing the degradome in different tissues. In a previous study (Gébelin et al., 2012), we described a surge in the accumulation of 21-nt sRNAs in the latex of TPD-affected rubber trees. Thanks to our miRNA annotation, we now show that these 21-nt sRNAs were not miRNAs derived from MIR genes. Therefore, miRNAs do not contribute to the 21 nt surge in sRNAs observed in TPD-affected trees. These 21 nt sRNAs could be classified as siRNAs, and may correspond to epigenetically-activated siRNA (easiRNA) (Borges & Martienssen, 2015). These siRNAs may be responsible for target transcript degradation through a pathway that is not directly dependent on miRNA-based cleavage (Bologna & Voinnet, 2014). This result underlines the relevance of analysing Hevea sRNAs beyond only the miRNA-based regulation pathways, in particular those produced by transposable elements.
Degradome analyses are generally performed with pooled tissue samples or with a single tissue. The present work highlighted tissue-specific regulation thanks to the parallel analysis of six tissue degradome libraries. Negative co-regulation of miRNA expression and its target is an important step in demonstrating the post-transcriptional control mechanism. For example, in response to saline stress, cleavage of chloroplastic CuZnSOD transcripts was correlated with upregulation of miRNA398 expression in bark and in roots (Gébelin et al., 2013a), but not in leaves. This regulation was not found in relation with the other stresses tested (drought, wounding, flooding, cold, high light intensity (Gébelin et al., 2013a)). In the degradome analysis, the target chloroplastic CuZnSOD was identified only in roots and bark with a small number of transcripts at the cleavage site, in perfect agreement with the results of co-regulation analyses. Although preliminary, our study enabled us to validate the link between degradome analysis and co-regulation analysis. It also underlines the need to account for the spatio-temporal dimension as well as the physiological context involving a case-by-case study. We also demonstrated the partial conservation of sRNA machinery and targets in Hevea. Conservation of miRNA-target pairs between plant lineages is well described (Cuperus, Fahlgren & Carrington, 2011). However, in a previous study, we showed partial conservation between Hevea and Arabidopsis miRNA binding sites for genes encoding ROS-scavenging enzymes (Gébelin et al., 2012;Zhang et al., 2019). Another example of partial conservation between plant is the targeting of AP2/ERF transcripts by miR172, which is well-documented in different plant species (Aukerman & Sakai, 2003;Wang et al.,

CONCLUSION
This study provided several new insights into Hevea concerning miRNA biogenesis and activity, tissue-specific post-transcriptional regulation in response to abiotic stress, and the impact of Hevea speciation with the whole genome duplication shared with cassava. The divergence observed from annual plant species suggests enhanced adaptation to stress in Hevea, particularly in laticifers, which are subjected to recurrent combinatorial stresses. Concerning the post-transcriptional regulation observed for sucrose metabolism in bark and in each step of the natural rubber biosynthesis pathway (Fig. 11), intra-and interspecific comparison of rubber-producing plants is needed to decipher the impact of post-transcriptional regulations in response to stress in the context of rubber production.
• Shuangyang Wu conceived and designed the experiments, analyzed the data, prepared figures and/or tables, authored or reviewed drafts of the paper, and approved the final draft.
• Benoît Farinas, Hélène Vignes and Kuswanhadi Kuswanhadi performed the experiments, prepared figures and/or tables, and approved the final draft.
• Stéphanie Pointet, Bénédicte Favreau, Shenghan Gao and Gaëtan Droc analyzed the data, prepared figures and/or tables, and approved the final draft.
• Enrique Ortega-Abboud, Songnian Hu and Chaorong Tang analyzed the data, authored or reviewed drafts of the paper, and approved the final draft.
• Jean-François Dufayard analyzed the data, prepared figures and/or tables, authored or reviewed drafts of the paper, and approved the final draft.
• Pascal Montoro conceived and designed the experiments, authored or reviewed drafts of the paper, and approved the final draft.

DNA Deposition
The following information was supplied regarding the deposition of DNA sequences: Data is available at the GSA and BIG

Data Availability
The following information was supplied regarding data availability: Degradome data and tissue-transcriptome data are available at the BIG Data Center Genome Sequence Archive: subCRA001506.