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

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.

Available data sets from five small-RNA-seq data were generated by (Gébelin et al., 2012; Gébelin et al., 2013b; Kanjanawattanawong et al., 2014; Lertpanyasampatha et al., 2012) comprising juvenile plants, latex of healthy and TPD-affected trees, young and fully developed leaves from mature trees. Complete tissue-specific (embryo, flower, latex, leaf, bark, root) transcriptomes were generated by Duan et al. (2013) and Piyatrakul et al. (2012).

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).

MIR gene identification

MIR genes were identified on the PB260 resequenced genome. MITP (https://sourceforge.net/projects/mitp/) was chosen with default parameters to identify hairpin structures and miRNAs in available small RNA datasets (Gébelin et al., 2012; Gébelin et al., 2013b; Kanjanawattanawong et al., 2014; Lertpanyasampatha et al., 2012). The small RNA–seq data described in Gébelin et al. are available under project number PRJCA001333 in the GSA (Wang et al., 2017a) and BIG Data Center (BIG Data Center Members, 2018). MITP uses RNAfold in the Vienna RNA Package (Piyatrakul et al., 2014) to detect hairpin structures. The longest hairpin sequences were further annotated using BLAST (Altschul et al., 1990) against MirBase (Kozomara & Griffiths-Jones, 2011) and plant non-coding RNA database (PNRD) (Yi et al., 2015) as references (Fig. 1).

Degradome libraries

Degradome libraries were constructed from a pool of total RNA from stressed and non-stressed plants, representing six different Hevea tissues (leaf, root, latex, bark, embryos and reproductive tissues) (see plant material), according to the protocol described in (German et al., 2009), with major modifications to account for the different sequencing techniques used (Miseq versus GAII in the 2009 publication). Briefly, libraries were produced from 150 µg of total RNA. The 5′ RNA-adapter sequence is RA5: 5′-GUUCAGAGUUCUACAGUCCGACGAUC- 3′. The 3′-adapter sequence is 5′-GTGACTGGAGTTCCTTGGCACCCGAGAATTCCATTTTTTTTTTTTTTTTTTV- 3′. The primer sequences for library amplification are 5′ Adapter Primer (RP1: 5′ AATGATACGGCGACCACCGAGATCTACACGTTCAGAGTTCTACAGTCCGA- 3′) with the 3′ Primer Adapter indexed (RPI1: 5′-CAAGCAGAAGACGGCATACGAGAT CGTGATGTGACTGGAGTTCCTTGGCACCCGAGAATTCCA- 3′, RPI2: 5′-CAAGCAGAA GACGGCATACGAGATACATCG GTGACTGGAGTTCCTTGGCACCCGAGAATTCC A-3′, RPI3: 5′-CAAGCAGAAGACGGCATACGAGATGCCTAA GTGACTGGAGTTCCT TGGCACCCGAGAATTCCA-3′, RPI4: 5′-CAAGCAGAAGACGGCATACGAGATTGGT CA GTGACTGGAGTTCCTTGGCACCCGAGAATTCCA-3′, RPI5: 5′-CAAGCAGAAGA CGGCATACGAGATCACTGT GTGACTGGAGTTCCTTGGCACCCGAGAATTCCA-3′, and RPI6: 5′-CAAGCAGAAGACGGCATACGAGATATTGGC GTGACTGGAGTTCCTT GGCACCCGAGAATTCCA-3′). The degradome data are available under project number PRJCA001333 in the GSA (Wang et al., 2017a) and BIG Data Center (BIG Data Center Members, 2018). The number of sequences produced and statistics are listed in Table S1.

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).

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).

Annotation of the miRNA biogenesis and regulatory pathway and phylogenetic analysis

A list of 95 Arabidopsis proteins involved in the RNA-mediated silencing machinery was extracted from the TAIR database (https://www.arabidopsis.org/) (Matzke & Mosher, 2014). The amino acid sequences were used to retrieve putative orthologous sequences using the SPALN pipeline with the parameters -O0 -Q4 -M5 -H180 (Iwata & Gotoh, 2012), from the PB 260 genomic sequences and the reference Hevea brasiliensis transcriptome (Piyatrakul et al., 2012). A set of 206 genes encoding putative AGO proteins from 23 plant species, and 75 genes and putative DCL proteins from 17 plant species, was used to locate candidate orthologous sequences in the rubber tree genome using the SPALN pipeline. The presence of signatures corresponding to conserved domains for either of these protein families, such as Piwi (PF02171), PAZ (PF02170) and Dicer dimerization (PF03368), was ascertained in the filtered genes using Interproscan. Putative H. brasiliensis AGO and DCL peptide sequences were classified through alignment against the corresponding full-length amino acid sequences from Populus trichocarpa, Ricinus communis and Arabidopsis thaliana (Daniel Rodríguez-Leal et al., 2016; Mukherjee, Campos & Kolaczkowski, 2013) using MUSCLE (version 3.8.31) (Edgar, 2004a; Edgar, 2004b). Cleaning was performed to remove gaps using trimAl (version 1.4) (Capella-Gutierrez, Silla-Martinez & Gabaldon, 2009), after which a phylogenetic tree was built using MrBayes (v3.2.6 x64).

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., 2012; Zha, 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 miRNA-directed 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 reveals that even miRNA-mediated cleaved targets predicted in weak degradome category can be biologically active.

At the transcriptome level, our analyses of the degradome products highlighted the prevalence of post-transcriptional regulatory events in latex cells and reproductive tissues in Hevea brasiliensis. Among the 136 transcript targets detected among all the tissues, 39 (28.7%) and 38 (27.9%) were detected in latex and reproductive tissues, respectively (Fig. 1, 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 (cysteine-containing 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 al. (2014), Piyatrakul et al. (2012), Putranto et al. (2015a), Putranto et al. (2015b), Tang et al. (2016), Venkatachalam, Thulaseedharan & Raghothama (2007) and Venkatachalam, Thulaseedharan & Raghothama (2009) (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 members of the ACO, ACS, ETR, CTR and EIN3 gene families), jasmonate biosynthesis (30 out of 50 genes tested, including a MYC transcription factor), and hormone signalling pathways (54 out of 149 including HbERF-IXb02; HbERF-VIIa07; HbERF-VIIa2; HbERF-VI05 HbERF-VIIa8, HbRAV-03; HbERF-VIIa5). In addition, among the 596 genes selected for their likely involvement in rubber biosynthesis pathway (Tang et al., 2016), 56 had degradation products compatible with post-transcriptional regulation (HbMVK1,2,3; HbMVD1,2; HbREF3,7,8; HbSRPP2; HbACAT3; HbHMGR1; HbHMGS1,2; HbPMD2 and HbGSP3 genes).

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 miRNA-mediated 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 38) (Table S14). No known TAS gene was identified in the Hevea genomic sequence. 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). Compared to other plant genomes, Hevea brasiliensis exhibits several AGO genes (13, including five putative 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).