Genome-wide analysis in Hevea brasiliensis laticifers revealed species-specific post-transcriptional regulations of several redox-related genes

MicroRNA-mediated post-transcriptional regulation has been reported on ROS production and scavenging systems. Although microRNAs first appeared highly conserved among plant species, several aspects of biogenesis, function and evolution of microRNAs were shown to differ. High throughput transcriptome and degradome analyses enable to identify small RNAs and their mRNA targets. A non-photosynthetic tissue particularly prone to redox reactions, laticifers from Hevea brasiliensis, revealed species-specific post-transcriptional regulations. This paper sets out to identify the 407 genes of the thirty main redox-related gene families harboured by the Hevea genome. There are 161 redox-related genes expressed in latex. Thirteen of these redox-related genes were targeted by 11 microRNAs. To our knowledge, this is the first report on a mutation in the miR398 binding site of the cytosolic CuZnSOD. A working model was proposed for transcriptional and post-transcriptional regulation with respect to the predicted subcellular localization of deduced proteins.

. Workflow diagram illustrating the main steps in the identification of redox-related genes in the Hevea reference genome sequence and transcriptome. Reference redox-related amino acid sequences were downloaded from the UniProt database. These sequences were blasted against the Hevea genome and transcriptome. Scaffolds harbouring Hevea redox-related genes were validated manually with ORF. Redoxrelated contigs were also identified using blastx and GO annotations of the Hevea transcriptome. The two lists of contigs were merged and blasted against the Hevea genome to identify unique contigs. They were then manually annotated with ORF and genome mapping. (2019) 9:5701 | https://doi.org/10.1038/s41598-019-42197-8 www.nature.com/scientificreports www.nature.com/scientificreports/ Comparative analysis of published latex transcriptomes. In order to identify redox-related genes expressed in latex, contigs or unigenes annotated as redox-related genes were extracted from the Supplemental Table 2 of recently published latex transcriptome analyses obtained by RNA sequencing technology [23][24][25][26][27][28] . For each publication, redox-related contigs or unigenes were assigned to one of the 30 gene families using their initial blastx annotation (Table 2). A small number of contigs (28, 30 and 12) was counted for three studies 23,24,28 compared to the total gene number found in this work (Table 1) and other transcriptome analyses (912, 77, 231) [25][26][27] . The transcriptome published by Wei and collaborators had the largest number of redox-related contigs (234) but a lower coverage (0.37 Gb for all samples) 27 . This transcriptome was obtained from trees of rubber clone RRIM 600 with long-term latex flow. For several gene families, the number of contigs was larger than the gene number counted in the reference genome. Tang and co-workers published transcriptome data for a mixture of several tissues including latex. Thus, the RNAseq dataset from clone PB 260 26 was adopted for further analysis for the following reasons: high coverage (6 Gb per sample), largest number of redox-related contigs (912), representation of all gene families, good statistical design with the use of 3 biological replicates, and data from a comparison of latex from healthy and TPD-affected trees.

ROS production
Tocotrienol biosynthesis GDH1 and GDH2 were located on scaffold1364_78602 and scaffold1364_29743. The phylogenetic analyses revealed a recent duplication of the genes (Supplemental Figs 4,8 and 11). Subcellular localization of redox-related genes was performed using WoLF PSORT, CELLO2GO and Plant-mPLoc. The largest number of proteins was predicted in chloroplast. Given that laticifers are non-photosynthetic tissues, chloroplast and plastid predictions were assigned as plastidic proteins. Subcellular localization of latex proteins was predicted as follows: 82 in plastids, 70 in cytosol, 12 in nucleus, 7 in mitochondria, 2 in extracellular, 1 in vacuole, 2 in peroxisome and 7 non-predicted.
When exploring RNAseq data from latex 26 , sixty transcripts were abundant (>1000 reads), and twelve of them were very abundant (>5000 reads) for one or other of the conditions. Twenty-nine transcripts were induced and forty-eight repressed in response to ethephon in healthy trees. Nine transcripts were induced in response to ethephon in TPD-affected trees. Four of these genes (PPO2, PrxQ, TrxS12 and TrxS13) showed contrasting regulation: repressed in healthy and overexpressed in TPD-affected trees. For the clarity of this manuscript, gene expression data are presented in Fig. 2 (cf. discussion section).
small RNA-mediated post-transcriptional regulation of redox-related genes. Redox-related transcripts targeted by microRNAs and ta-siRNAs were searched using CLEAVELAND pipeline 13 in the degradome dataset obtained from various tissues (root, leaf, bark, latex, flowers and embryo) and the reference transcriptome for rubber tree clone PB 260. The degradome analysis did not revealed post-transcriptional regulations by ta-siRNA (data not shown). Of the 407 redox-related genes, 13 were targeted by 11 different microRNAs   Table 2. Annotation of Hevea latex redox-related genes from published latex transcriptomes. www.nature.com/scientificreports www.nature.com/scientificreports/ (Table 3). The degradome analysis revealed post-transcriptional regulation of these transcripts at spatial level (thanks to tissue-specific libraries) and cleaved transcript abundance level classified into degradome categories. The number of microRNA families was different for each tissue: 1 in roots, 7 in leaves, 7 in latex, 4 in bark and 2 in flowers. In latex, seven of these redox-related genes were targeted by 6 families of microRNAs. Three known families of microRNA (miR535, miR398b and miR394) targeted and cleaved transcripts from genes APX3 (ascorbate peroxidase 3), SOD2 (Cu/Zn superoxide dismutase 2), GR1 and GR2 (glutathione reductase 1 and 2), respectively. For transcripts from gene APX3, strong spatial regulation was observed with a greater abundance of miR535 in leaf (degradome category 0) compared to latex (degradome category 4). One transcript, encoding MPBQ/MSBQ methyltransferase 3 involved in tocotrienol biosynthesis, was also cleaved in bark and leaf by a new microRNA named miRNAn7. For miR398b, which cleaves chloroplastic Cu/Zn superoxide dismutase transcripts 29 , a low abundance of cleaved transcripts was found in latex (degradome category 3) and root (degradome category 4). Interestingly, the three cytosolic isoforms were not detected in the degradome libraries confirming the previous observation made by Gébelin and co-workers 11 . The miR398 binding site was further scanned and sequence variations were observed in the 5′ and 3′ seed regions but also between the very sensitive 10 th and 11 th nucleotides of the microRNA sequence targeting HbCuZnSOD1 and at the 12 th nucleotide in both HbCuZnSOD3 and HbCuZnSOD4 sequences (Table 4). Three microRNAs (miRNAn1 to 3), with cleavage activity in latex, were new microRNAs not yet annotated in the miRBase database (Table 3). Catalase 1, the unique cytoplasmic Cat gene showing the highest Cat expression in latex, was regulated by a new microRNA named miRNAn1. This mRNAn1 also targeted peroxisomal catalase 2. Cytosolic glutathione reductase 1 and plastidic glutathione reductase 2 were highly expressed in latex and targeted by miR394.
The expression of the 13 post-transcriptionally regulated genes was recalculated using the reads covering the cleavage site only, in order to check if the level of expression assessed by the number of reads describes the real functionality of mRNAs (Supplemental Tables 3 and 4). The expression of 8 of the 13 targeted transcripts were significantly affected by the new way of calculation. Significant fold changes observed in standard RNA sequencing for ethephon treatment or TPD occurrence disappeared for genes APX3, GR1, MDHAR2, Prx2C1, Px1, Px6 and VTE4 when using the number of reads covering the cleavage site to calculate the expression level. Finally, some effects of ethephon were maintained for GSTF1 and Prx2C1.

Discussion
Apart from plant model species, this study is the most complete genome-wide analysis of ROS production and scavenging systems and antioxidant biosynthesis in a perennial crop. The main 30 redox-related gene families totalize 407 genes in Hevea. This is a larger number of genes compared to Arabidopsis especially due to the expansion of peroxidase genes in Hevea. Based on the RNAseq dataset, small RNA/target identification, and prediction of subcellular localization, a model of transcriptional and post-transcriptional regulations of the 161 redox-related genes expressed in latex was attempted for a rubber clone particularly prone to oxidative stress (Fig. 2). The redox-related proteins were predominantly localized in plastids (82 proteins) and cytosol (70 proteins). This comprehensive analysis highlighted critical steps of redox homeostasis in latex.
This study also revealed specific regulation of ROS-scavenging systems, which might be adapted to strong and steady ROS production in latex cells due to recurrent harvesting stress and latex regeneration between two tappings. Lutoids are polydispersed vacuoles with lysosomal properties. Previous biochemical studies revealed that NADPH oxidase is the main source of ROS in laticifer cytosol, especially under stress 18 . The present study revealed that this enzyme was mostly encoded by the Rboh2 gene in latex cells, and its expression was enhanced by ethephon application. These results suggest that Rboh2 encodes the main enzyme generating ROS at the outer surface of lutoid membrane in contact with cytosol.
This production of O 2 − in cytosol requires a powerful detoxification system in this compartment. Superoxide dismutase is the enzyme involved in the first step of detoxification inducing the dismutation of the superoxide anions, produced by the lutoid NADPH oxidase, into hydrogen peroxide 22 . CuZnSOD1 transcripts were much more abundant compared to other genes encoding SOD. Unlike Arabidopsis, none of the Hevea cytosolic SOD isoforms was subjected to post-transcriptional regulation by miR398. A mutation in the binding site makes miR398 ineffective. The high expression of the CuZnSOD1 gene might then support the maintenance of SOD activity and a consequent high level of anion superoxide dismutation. To demonstrate the biological relevance of post-transcriptional regulations, the physiological context (type and duration of stress) in which the regulation occurs should be further identified case by case. For example, the cleavage of the chloroplastic CuZnSOD transcripts was correlated with the upregulation of miRNA398 expression in response to a salt stress specifically in bark and root 29 .
The second step deals with the decomposition of H 2 O 2 to H 2 O and O 2 through five hydrogen peroxide scavenging pathways coexisting in cytosol (peroxidase, ascorbate peroxidase, glutathione peroxidase, peroxiredoxin www.nature.com/scientificreports www.nature.com/scientificreports/ and catalase). High and steady ROS production in latex cells requires Cat activity, which generally comes into play under stress. A decrease in Cat activity was recorded in TPD-affected trees enabling the general oxidative stress in latex cells 35 . Cat1 gene was highly expressed in latex and might be the main gene related to the Cat activity. Although post-transcriptional regulation was shown by microRNA miRn1, this microRNA did not efficiently cleave Cat1 transcripts in the tested biological conditions (low number of read ends at the cleavage site in degradome data). For the genes encoding thioredoxins, TrxH5 had the highest level of expression out of the 161 genes expressed in latex. From our knowledge, there is no published information related to the potential role of Prx in latex and further characterization is required. The ascorbate/glutathione cycle, involving in its last lines APx and GPx, is essential in the reduction of H 2 O 2 to H 2 O and O 2 . Regeneration of the ASA and GSH forms reduced by the ascorbate-glutathione cycle involved several enzymes encoded by MDHAR2, DHAR2, GR1 and GR2. The ethephon treatment did not transcriptionally activate genes involved in the glutathione/ascorbate cycle. Although some post-transcriptional regulations appeared in the degradome analysis showing that both the GR1 and GR2 transcripts, miR394 did not significantly cleave GR transcripts. APx has a high affinity for H 2 O 2 and can reduce it to H 2 O in chloroplasts, cytosol, mitochondria and peroxisomes, as well as in the apoplastic space. Of the three genes encoding a cytoplasmic ascorbate peroxidase, the HbAPx1 and HbAPx5 transcripts were the most abundant. Considering the lower expression of these 3 APx genes compared to the plastidic APx4, the cytosolic ASA pathway might have a lower reducing capacity than the plastid pathway, which is obvious since the production of ROS is known to be high in plastids. Of the 23 Hevea genes encoding a GST, 21 were predicted as cytosolic GST. Among them, the GSTU3 and GSTF1 genes were actively expressed in latex cells. As GST plays a central role in the use of the reduction power of GSH to detoxify electrophiles, glutathione might be considered as the most important antioxidant in laticifers.
Glutathione, ascorbate and vitamin E isomers are the major antioxidants in latex 22 . The glutathione biosynthesis pathway involves two ATP-dependent enzymes: γ-glutamate cysteine ligase (GCL) and glutathione synthetase (GS). Of the two GS and GCL genes identified in the rubber genome, only one of each was encoded protein predicted to be expressed in latex cytosol (GS1 and GCL2), one GS (GS2) and the two GCL (GCL1 and GCL2) being expressed in plastids. The genes encoding GS2 and GCL2 were significantly over-expressed in response to ethephon. There are four routes for ASA biosynthesis in plant: the L-galactose pathway, the myo-inositol oxygenase pathway, the salvage pathway via L-galactonate, and the L-gulose-pathway. Of these four routes, L-galactose is the major pathway in many plants 36,37 . The L-galactono-1 4-lactone (L-GalL) biosynthesis pathway occurs in cytosol, which consists of five enzymes (VTC1, GME, VTC2, VTC4 and GDH). All genes encoding these enzymes have homologues expressed in latex cytosol.
www.nature.com/scientificreports www.nature.com/scientificreports/ its targeting by miRNAn7. As γ-tocotrienol is the most abundant vitamin E isomer, its accumulation might be fostered by the weak capacity to produce α-tocotrienol and tocopherol.
To conclude, this study reveals new insights into small RNA-mediated post-transcriptional regulations of ROS-scavenging systems. To our knowledge, this is the first report on a mutation in the miR398 binding site of the CuZnSOD altering the post-transcriptional regulation described in model species. In addition, the literature mentioned microRNA-mediated post-transcriptional regulation on ROS production and scavenging systems. This work paves the way to the study of adaptive mechanisms. Besides, several genetic studies have revealed the involvement of antioxidant compounds in complex traits of several species [40][41][42][43] . In Hevea, the 161 redox-related genes expressed in latex represent candidate genes for the identification of allelic variability. The development of molecular markers and the analysis of genetic variability of antioxidants should support breeding programmes, especially for traits relative to environmental stress.

Methods
Identification and classification of redox-related genes in the Hevea brasiliensis genome and transcriptome. Redox-related genes were identified from both the Hevea reference genome and transcriptome (Fig. 1). An amino acid sequence dataset was created by downloading sequences of thirty redox-related gene families from the UniProt database (http://www.uniprot.org/) using published accession numbers mostly from Arabidopsis, except for the polyphenol oxidase (PPO) family, which is absent in Arabidopsis. Sixteen families were selected for ROS production and scavenging (Table 1). In addition, protein sequences of genes involved in the biosynthesis of three major antioxidants in latex (ascorbate, glutathione, and tocotrienol) were collected. This dataset was blasted against the published Hevea genome 25 and transcriptome 26 . Redox-related contigs were also identified using blastx and GO annotations of the Hevea transcriptome. The two lists of contigs were merged and blasted on the rubber genome to identify unique contigs. Redox-associated genes were classified for each gene family related to ROS production, ROS-scavenging and regulation, and antioxidant biosynthesis (ascorbate, glutathione and tocotrienols). phylogenetic analysis of redox-related genes. The full length amino acid sequences of Arabidopsis redox-related protein were aligned with the amino acid deduced sequences from Hevea clone Reyan 7-33-97 genome. Identities of proteins are provided in Supplemental Table 5. The polyphenol oxidase family being absent in Arabidopsis, we used the Populus PPO gene family. This alignment was made by Muscle via Mega 6 44 . Amino acid sequence of Arabidopsis actin 1 or Arabidopsis glutamate cysteine ligase was used as outgroup control. The phylogenetic trees were generated in Mega 6 by Bootstrap method with 500 replications after alignment. prediction of the subcellular localization of redox-related proteins. The subcellular location of redox-related genes was predicted with translated sequences using WoLF PSORT (http://www.genscript.com/ wolf-psort.html), CELLO2GO (http://cello.life.nctu.edu.tw/cello2go/) and Plant-mPLoc (http://www.csbio.sjtu. edu.cn/bioinf/plant-multi/). The 3 predictors were successfully tested on subcellular localization prediction 45 . The matching ratio between the prediction result and protein location was calculated according to Xiong's Supplemental Table 2. The matching ratios from these 3 predictors ranged from 50% to 80%. The prediction of subcellular localization was considered as acceptable when the matching ratio of merged results was above 90%.
Identification of small RNA and target mRNA couples. Degradome data for several Hevea tissues (latex, leaves, male and female flowers, seeds, root, bark and somatic and zygotic embryos) were obtained according to a protocol adapted from German 46 . Hevea microRNAs from small RNAseq data published by Gébelin and co-workers 11,29 were annotated by MITP (https://sourceforge.net/projects/mitp/files/). This pipeline complies to the recommendations set by Axtell and coll 47 , looking from hairpin structures, producing miRNA and miRNA* with up to 3 bulges or 6 unpaired bases between miRNA and miRNA*. The prediction was done with sequences  Table 4. Comparison of HbmiR398 (acc_420) cleavage site between cytosolic and chloroplastic CuZnSOD isoforms. Arrow indicated the cleavage site observed experimentally for HbCuZnSOD2 by miR398 (Gébelin et al. 2012) and in the degradome analysis. Sequence variations in cytosolic isoforms sequences compared to HbCuZnSOD2 are in bold and highlighted character.