# Transcriptome analysis in Hevea brasiliensis latex revealed changes in hormone signalling pathways during ethephon stimulation and consequent Tapping Panel Dryness

Pascal Montoro, Shuangyang Wu, Bénédicte Favreau, Eva Herlinawati, Cécile Labrune, Marie-Laure Magniette, Stéphanie Pointet, Maryannick Rio, Julie Leclercq, Sigit Ismawanto, et al.

<figure><img src="image_1.png" /></figure>

## ► To cite this version:

Pascal Montoro, Shuangyang Wu, Bénédicte Favreau, Eva Herlinawati, Cécile Labrune, et al. Transcriptome analysis in Hevea brasiliensis latex revealed changes in hormone signalling pathways during ethephon stimulation and consequent Tapping Panel Dryness. Scientific Reports, 2018, 8, pp.8483. 10.1038/s41598-018-26854-y. hal-02624945

HAL Id: hal-02624945

https://hal.inrae.fr/hal-02624945

Submitted on 26 May 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
SCIENTIFIC REPORTS

<figure><img src="image_3.png" /></figure>

OPEN

# Transcriptome analysis in *Hevea brasiliensis* latex revealed changes in hormone signalling pathways during ethephon stimulation and consequent Tapping Panel Dryness

Pascal Montoro<sup>1,2</sup>, Shuangyang Wu<sup>1,2,3,4</sup>, Bénédicte Favreau<sup>1,2</sup>, Eva Herlinawati<sup>5</sup>, Cécile Labrune<sup>6</sup>, Marie-Laure Martin-Magniette<sup>6,7,8</sup>, Stéphanie Pointet<sup>1,2</sup>, Maryannick Rio<sup>1,2</sup>, Julie Leclercq<sup>1,2</sup>, Sigit Ismawanto<sup>5</sup> & Kuswanhadi<sup>5</sup>

Tapping Panel Dryness (TPD) affects latex production in *Hevea brasiliensis*. This physiological syndrome involves the agglutination of rubber particles, which leads to partial or complete cessation of latex flow. Latex harvesting consists in tapping soft bark. Ethephon can be applied to stimulate latex flow and its regeneration in laticifers. Several studies have reported transcriptome changes in bark tissues. This study is the first report on deep RNA sequencing of latex to compare the effect of ethephon stimulation and TPD severity. Trees were carefully selected for paired-end sequencing using an Illumina HiSeq 2000. In all, 43 to 60 million reads were sequenced for each treatment in three biological replicates (slight TPD trees without ethephon stimulation, and slight and severe TPD trees with ethephon treatment). Differentially expressed genes were identified and annotated, giving 8,111 and 728 in response to ethephon in slight TPD trees and in ethephon-induced severe TPD trees, respectively. A biological network of responses to ethephon and TPD highlighted the major influence of metabolic processes and the response to stimulus, especially wounding and jasmonate depression in TPD-affected trees induced by ethephon stimulation.

Received: 28 June 2017

Accepted: 21 May 2018

Published online: 31 May 2018

Natural rubber produced by *Hevea brasiliensis* accounts for more than 40% of rubber consumption worldwide. With more than 10 million hectares, rubber cultivation and its industry play an important socio-economic role in rubber producing countries. Natural rubber is synthesized in laticifers, which form an articulated network of anastomosed latex cells. Latex containing natural rubber is collected by tapping the bark of rubber trees. Latex is the cytoplasm of laticifers. Ethylene induces several metabolic pathways in latex. For certain *Hevea* clones, the application of ethephon, an ethylene releaser, stimulates latex regeneration between two tappings and latex flow after tapping, hence latex production. The growing demand for natural rubber prompts farmers to use quick-starter clones and intensive harvesting systems, which have led to increasingly large losses in rubber production. The physiological syndrome called Tapping Panel Dryness (TPD) is responsible for an estimated 12 to 20% loss in rubber production<fnref n="1" />, which could increase in a context of climate change.

Environmental and harvesting stresses can lead to oxidative stress in laticifers and consequently to TPD induction². A lutoid is a single-membrane microvacuole with lysosomal characteristics. Lutoidic NAD(P)H oxidase is a major source of superoxide anion species in latex (Chrestin *et al.* 1984). This production of reactive

<footnote n="1">¹CIRAD, UMR AGAP, F-34398, Montpellier, France. ²AGAP, Univ Montpellier, CIRAD, INRA, Montpellier SupAgro, Montpellier, France. ³CAS Key Laboratory of Genome Sciences and Information, Beijing Institute of Genomics, Chinese Academy of Sciences, Beijing, 100101, China. ⁴University of Chinese Academy of Sciences, Beijing, 100049, China. ⁵Indonesian Rubber Research Institute, Sembawa Research Centre, Palembang, Indonesia. ⁶Institute of Plant Sciences Paris Saclay IPS2, CNRS, INRA, Université Paris-Sud, Université Evry, Université Paris-Saclay, Bâtiment 630, 91405, Orsay, France. ⁷Institute of Plant Sciences Paris-Saclay IPS2, Paris Diderot, Sorbonne Paris-Cité, Bâtiment 630, 91405, Orsay, France. ⁸UMR MIA-Paris, AgroParisTech, INRA, Université Paris-Saclay, 75005, Paris, France. Correspondence and requests for materials should be addressed to P.M. (email: pascal.montoro@cirad.fr)</footnote>
oxygen species (ROS) induces peroxidatic degradation of the unsaturated lipids of the membrane. Such lutoid instability leads to the release of lutoid proteins in cytosol and promotes rubber particle agglutination³. Under exacerbated conditions, an irreversible type of TPD, called brown bast (BB), occurs. This form is related to deformation of the bark due to thylosoid formation, lignified gum, and abnormal division of parenchyma cells, and may be related to a cyanogenesis process⁴,⁵. These two types of TPD were recently called ROS-TPD and BB-TPD⁶. TPD, trunk phloem necrosis, brown bast, dry bark, etc. show the same anatomical symptom⁶. However, these syndromes might be related to different stress origins such as rootstock-scion incompatibility⁷, soil compaction, and environmental or harvesting stress⁸.

Numerous molecular studies have led to the identification of TPD-regulated genes, such as *HbMyb1*, *HbTOM20* and *HbTCTP1*¹,⁹⁻¹⁶. Differential transcript abundance was observed for twenty-seven candidate genes, related to TPD occurrence, in latex and phloem tissues for ROS-scavenging, ethylene biosynthesis and signalling genes⁶. The technologies used in these previous studies led to partial identification of genes associated with TPD occurrence. More recently, the technology of RNA sequencing has opened up new perspectives for a full comprehensive analysis of TPD through a transcriptome comparison between TPD-free and TPD-affected trees. Given the stoppage in latex flow, deep sequencing of RNA has focused solely on bark tissues of severe TPD-affected trees. These analyses notably revealed a down-regulation of genes related to rubber biosynthesis¹⁷ as well as jasmonate production¹⁸. A few studies have attempted to analyse latex from early TPD-affected trees. The first one revealed post-transcriptional regulation with both RNA degradation and a shift in miRNA biogenesis¹⁹. The second suggested differential expression of proteins related to rubber biosynthesis, the reactive oxygen metabolism and cell apoptosis²⁰.

This paper sets out to provide new insight into the molecular regulatory mechanism of ethephon-induced TPD in latex and to identify a set of differentially expressed genes associated with ethephon-induced TPD occurrence. The study identified a set of differentially expressed genes associated with TPD occurrence. Several earlier publications tried to study TPD, but at the last stage of TPD. They therefore only obtained 100% TPD-affected trees and had to study bark tissue. The novelty of our study is to carefully monitor the onset of the syndrome and manage to select trees with slight and severe TPD symptoms still able to produce some latex. Regular monitoring of the dry cut length (DCL) on the tapping panel revealed that all the trees had a variation in DCL value⁶. High latex viscosity in a localized part of the bark reduces latex flow. This transient symptom concerns a small proportion of the tapping cut. In contrast, the ROS-TPD type resulting from the physiological fatigue of laticifers concerns a larger share of the tapping cut, which increases regularly over the tapping period and requires a tapping rest of several months for total recovery. Ethephon treatment accelerates the latex metabolism and induces earlier and greater TPD occurrence, especially for rubber clone PB 260 with a high latex metabolism⁶. Firstly, the effect of ethephon was analysed on slight TPD-affected trees. Secondly, slight and severe TPD-affected trees subjected to ethephon stimulation were compared. RNA sequencing was carried out in three biological replicates and a differential expression analysis was carried out. This study identified 728 differentially expressed genes associated with ethephon-induced TPD, 272 of which were specifically regulated by TPD and not by ethephon treatment. When ethephon-induced TPD became established, an overall decrease in gene expression was observed, except for genes related to the GO "stimulus response" term, particularly those related to the wounding and jasmonic acid pathways.

## Results

**RNA sequencing and mapping on the reference transcriptome.** Latex flows out from laticifer vessels after tapping soft bark. Latex drops appear instantly on the tapping cut. The dry cut length (DCL) was recorded just after tapping and before latex flowed along the cut. Most trees showed a partial DCL ranging from 0 to 100% (Fig. 1). Slight TPD-affected trees (DCL < 25%, Fig. 1B) were found in trees that were treated and not treated with ethephon. By contrast, severe TPD-affected trees (DCL > 75%, Fig. 1D) were identified only in ethephon-treated trees. The DCL was not strictly stable for each tree due to localized and partial TPD occurrence on the tapping panel (Suppl. Data 1). For that reason, trees were classified as slight or severe TPD depending on the mean DCL over the last 12 months of tapping before sampling (Table 1).

RNAs from the latex of non ethephon-treated trees (control), and both slight and severe ethephon-induced TPD trees were paired-end sequenced using an Illumina HiSeq 2000. In all, 43 to 60 million reads were obtained for each oriented TRUSeq library (Suppl. Data 2). Trimmed reads were mapped on the rubber clone PB 260 reference transcriptome²¹. Of the 86,941 contigs, 72,910 were previously annotated by BLASTX. Reads were counted for the 86,941 contigs.

**Identification of differentially expressed genes in response to ethephon and TPD severity.** A differential analysis was performed for the 86,941 contigs (Suppl. Data 3). Given the importance of small changes in RNA expression²², differentially expressed genes (DEGs) were selected when their FDR adjusted *p*-value was lower than 0.05. In latex, 8,111 DEGs were detected in response to ethephon stimulation with a similar number of over-expressed (3883) and under-expressed (4,228) genes (Suppl. Data 4, Table 2).

Ethephon-induced TPD-affected trees were identified in the set of ethephon-treated trees only. An analysis of slight and severe TPD-affected trees from the same set of ethephon-treated trees revealed 728 DEGs consisting of 97 over- and 631 under-expressed genes (Suppl. Data 5). Unlike DEGs identified in response to ethephon, 6 times more under-expressed genes than over-expressed genes were identified in trees with severe TPD. The expression of 456 of these ethephon DEGs (over- and under-expressed) was accentuated in response to TPD from ethephon to TPD for both over-expressed and under-expressed genes by ethephon were more highly expressed under TPD as well as for under-expressed genes. Of the 728 DEGs associated with severe TPD, expression of 272 genes was not significantly changed by ethephon treatment in slight TPD trees (Suppl. Data 6). The 272 DEGs were mostly under-expressed (244) rather than over-expressed (28).
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Gradient of Tapping Panel Dryness (TPD) symptoms in *Hevea brasiliensis*. (A) Normal latex flow along the entire tapping cut. (B) Tapping cut of a healthy tree. (C) Tapping cut of a slight TPD-affected tree. (D) Tapping cut of severe TPD-affected tree.</figcaption></figure>

<table><caption>**Table 1.** Phenotype of trees selected for the RNA sequencing experiment. Trees were tapped every two days for 3 years. Ethephon stimulation was applied 12 times a year. The dry cut length (DCL) was monitored monthly. Trees were considered as slightly affected when DCL < 25% and severely affected by TPD when DCL > 75%.</caption><thead><tr><th>Class of trees</th><th>Ethephon treatment (%)</th><th>DCL (%)</th><th>Tree phenotype</th></tr></thead><tbody><tr><td rowspan="3">Control, no-ET</td><td>0</td><td>4</td><td>Slight TPD</td></tr><tr><td>0</td><td>4</td><td>Slight TPD</td></tr><tr><td>0</td><td>4</td><td>Slight TPD</td></tr><tr><td rowspan="3">Slight TPD, ET-treated</td><td>2.5</td><td>6</td><td>Slight TPD</td></tr><tr><td>2.5</td><td>15</td><td>Slight TPD</td></tr><tr><td>2.5</td><td>16</td><td>Slight TPD</td></tr><tr><td rowspan="3">Severe TPD, ET-treated</td><td>2.5</td><td>95</td><td>Severe TPD</td></tr><tr><td>2.5</td><td>80</td><td>Severe TPD</td></tr><tr><td>2.5</td><td>81</td><td>Severe TPD</td></tr></tbody></table>

<table><caption>**Table 2.** Effect of ethephon and severity of TPD occurrence on the number of differentially expressed genes (DEGs).</caption><thead><tr><th rowspan="2">Level of expression</th><th colspan="3">Number of DEGs</th></tr><tr><th>Ethephon vs control (slight TPD trees)</th><th>Severe vs slight TPD (ethephon-treated trees)</th><th>TPD-specific (ethephon-treated trees)</th></tr></thead><tbody><tr><td>Over-expressed</td><td>3883</td><td>97</td><td>28</td></tr><tr><td>Under-expressed</td><td>4228</td><td>631</td><td>244</td></tr><tr><td>Total</td><td>8111</td><td>728</td><td>272</td></tr></tbody></table>

A principal component analysis (PCA) of the read data was carried out for 8,111 DEGs. The PCA highlighted the variation between the three types of samples (control, slight TPD trees treated with ethephon, severe TPD-affected trees treated with ethephon) (Suppl. Data 7). The first two components explained 92% of the variability between samples, discriminating between slight and severe TPD trees on the first axis (62%), and
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Heatmap hierarchical clustering of the 8,300 contigs based on the count distribution of the nine samples (three biological replicates × three treatments). The treatments are latex from control, slight TPD trees with ethephon stimulation, and latex from severe TPD-affected trees with ethephon stimulation.</figcaption></figure>

slight TPD-affected trees treated or not with ethephon (control), on the second axis (30%). Heatmap hierarchical clustering of the DEGs was used to visualize count distribution in the nine samples (Fig. 2). The samples were clustered by treatment, with two main clusters representing severe and slight TPD trees, the latter being sub-clustered according to ethephon treatment. The heatmap highlighted specific gene expression patterns for each of the three treatments. Latex samples from control trees without ethephon treatment had a large number of over-expressed genes, while ethephon-treated latex samples revealed overall up-regulation of genes, especially for genes under-expressed in the treatment without ethephon. Conversely, gene expression in latex from ethephon-induced severe TPD trees was mainly under-expressed.

**Gene ontology, classification and biological network of response to ethephon and TPD severity.**

In order to improve the annotation of DEGs, contig sequences were blasted using scaffolds from published genomes as the target²³. An InterProScan was performed with corresponding peptide sequences from the rubber reference genome. InterProScan GOs were merged to annotation in order to run Blast2GO analyses (Suppl. Data 8). The first two levels of gene ontology (GO) classification were presented for up- and under-expressed genes in response to ethephon and TPD severity (Table 3). A large number of ethephon-regulated DEGs was identified (>500) for the cellular process, metabolic process, binding and catalytic activity GO terms. More GO terms were counted for ethephon over-expressed genes compared to under-expressed genes, except for the "response to stimulus", "single organism process", "membrane" and "membrane part" terms. For all the ethephon-regulated genes, most of the GO terms had more than 100 DEGs except for "response to stimulus", "membrane part", "organelle part", "nucleic acid binding transcription factor activity" and "transporter activity". Ethephon under-expressed genes had fewer GO terms. For ethephon-induced TPD-related DEGs, GO annotation showed that most of the ethephon over-expressed genes participated in the biological process and molecular function. Level 2 GO terms belonged to "biological regulation", "cellular process", "metabolic process", "regulation of biological process", "binding and catalytic activity". TPD under-expressed genes were related to most of the GO terms except "nucleic acid binding transcription factor activity" and "transporter activity" in the molecular function group.

A network analysis was carried out for the four DEG groups. For DEGs over-expressed by ethephon, biological processes including "metabolic process" (primary metabolism, macromolecule, protein, etc.), and "response to stimulus" (hormone, abiotic, biotic, oxidative stress) were overrepresented (Fig. 3, Suppl. Data 9). For DEGs under-expressed by ethephon, similar functions were reversely regulated (Fig. 4, Suppl. Data 10). For DEGs over-expressed in severe ethephon-induced TPD, a much smaller number of genes was significantly annotated (Fig. 5, Suppl. Data 11). A few DEGs, related to response stimulus with, in particular, jasmonic acid, chitin, wounding and fungus, were still activated. Some genes related to the regulation of the biological process were also activated, in particular the regulation of response to stimulus and metabolic processes. For the last ethephon-induced TPD under-expressed DEG group, a very small number of DEGs was classified: GO terms were related to "response to stimulus" (stress and abiotic), "developmental process" and "secondary metabolic process" (phenylpropanoid metabolic process) (Fig. 6, Suppl. Data 12).
<table><caption>**Table 3.** Number of Differentially Expressed Genes (DEGs) for the first two levels of Gene Ontology (GO) for over- and under-expressed genes, in response to ethephon and severity of ethephon-induced TPD.</caption><thead><tr><th rowspan="2">GO Name</th><th rowspan="2">Level 2</th><th colspan="4">Number of DEGs</th></tr><tr><th colspan="2">Ethephon vs control</th><th colspan="2">Severe vs slight TPD</th></tr><tr><th>Level 1</th><th></th><th>Over-expressed</th><th>Under-expressed</th><th>Over-expressed</th><th>Under-expressed</th></tr></thead><tbody><tr><td rowspan="7">Biological process</td><td>biological regulation</td><td>187</td><td>135</td><td>11</td><td>37</td></tr><tr><td>cellular process</td><td>740</td><td>472</td><td>24</td><td>110</td></tr><tr><td>localization</td><td>111</td><td>99</td><td>4</td><td>35</td></tr><tr><td>metabolic process</td><td>768</td><td>571</td><td>28</td><td>143</td></tr><tr><td>regulation of biological process</td><td>182</td><td>129</td><td>11</td><td>37</td></tr><tr><td>response to stimulus</td><td>0</td><td>86</td><td>4</td><td>27</td></tr><tr><td>single-organism process</td><td>305</td><td>317</td><td>9</td><td>84</td></tr><tr><td rowspan="7">Cellular component</td><td>cell part</td><td>355</td><td>158</td><td>2</td><td>48</td></tr><tr><td>cell</td><td>355</td><td>158</td><td>2</td><td>48</td></tr><tr><td>membrane part</td><td>72</td><td>95</td><td>0</td><td>22</td></tr><tr><td>membrane</td><td>121</td><td>149</td><td>4</td><td>39</td></tr><tr><td>macromolecular complex</td><td>209</td><td>48</td><td>0</td><td>14</td></tr><tr><td>organelle</td><td>257</td><td>99</td><td>2</td><td>34</td></tr><tr><td>organelle part</td><td>96</td><td>35</td><td>0</td><td>10</td></tr><tr><td rowspan="4">Molecular function</td><td>binding</td><td>912</td><td>768</td><td>27</td><td>183</td></tr><tr><td>catalytic activity</td><td>589</td><td>583</td><td>23</td><td>143</td></tr><tr><td>nucleic acid binding transcription factor activity</td><td>0</td><td>0</td><td>6</td><td>0</td></tr><tr><td>transporter activity</td><td>0</td><td>0</td><td>3</td><td>0</td></tr></tbody></table>

**Identification and regulation of DEGs encoding transcription factors in response to ethephon and severity of TPD.** Transcription factors are crucial regulators of gene expression. They were extracted from the two classes of ethephon-induced TPD treatments in Supplementary Data 13 and 14, respectively. They were classified according to their membership of the 13 most represented families of transcription factors (Table 4). The top five families in response to ethephon treatment were Ethylene Response Factors (ERF; 39 contigs), followed by Auxin Response Factor (ARF; 33 contigs), WRKY (31 contigs), Homeobox superfamily (31 contigs) and bHLH (29 contigs). In addition, WRKY (4 contigs), Myb (2 contigs), bHLH (2 contigs), ERF (1 contig) and MYC (1 contig) were the top ten most over-expressed genes in response to ethephon treatment (Supplementary Data 13). For each family of transcription factors, some genes were over- or under-expressed at a similar level, except for ARF, homeobox and MADS box, which had more under-expressed genes (24/33, 34/31 and 8/12, respectively).

The comparison of slight and severe ethephon-induced TPD trees revealed a dramatic drop in the number of regulated transcription factors from 316 to 50 (Table 4). The under/over-expressed transcription factor ratio increased from 1.38 (183/133) in slight TPD-affected to 3.54 (39/11) in severe TPD-affected trees, highlighting a general decrease in plant signalling. Most of the TF families were under-expressed, except the GRAS family, which belongs to the SCARECROW (SCR) subfamily. They were identified in latex but were reversely regulated for slight TPD-affected (3 under-expressed) and severe TPD trees (2 over-expressed).

**Discussion**

This study reports on the first deep RNA sequencing in the latex of ethephon-induced TPD-affected trees. It was made possible by the selection of slight and severe TPD-affected trees producing latex, based on the monitoring of TPD occurrence in a large panel of trees for 3 years. Latex from TPD trees was collected just before the tapping cut became totally dry. Several other studies were previously carried out on latex from healthy trees only, notably to understand molecular mechanisms related to rubber production<sup>21,24-28</sup>. With regard to TPD, some gene expression markers, such as *HbMyb1*, *HbTOM20* and *HbTCTP1*, have been identified through differential display RT-PCR, microarrays and suppression subtractive hybridization<sup>9-11,15,29</sup>. The expression of six *MYB* genes was changed in ethephon-induced TPD-affected trees (one over-expressed and five under-expressed). In contrast, none of the other candidate genes was associated with TPD, neither *HbTOM20* (2 contigs (CL1183Contig1 and CL1183Contig4) not expressed in latex) nor *HbTCTP1*, despite the high expression of one of the two contigs, CL1Contig21026 (107-142 reads) and CL1Contig7138 (41881-52466 reads). This variability in gene expression markers reveals that some candidate genes are associated with either clonal response or certain environmental or experimental conditions, and are not directly related to the biological process studied. In addition, these technologies provided a partial and a-priori view of transcriptome expression. More recently, high-throughput sequencing was used to describe global gene expression changes in bark<sup>17,18</sup>. Ethephon-induced TPD resulted in an overall decrease in gene expression, which was less in latex (631 DEGs) than in bark (8,137 and 22,577 DEGs)<sup>17,18</sup>. The down-regulation of some genes related to rubber biosynthesis in bark<sup>17,18</sup> was confirmed at protein level<sup>20</sup>. Bark tissue consists of several cellular types (phloem, laticifer, vascular rays, sclereid). The larger number of over-expressed DEGs in bark, compared to latex, suggests that either at least one bark tissue type other
<figure><img src="image_6.png" /><figcaption>**Figure 3.** BINGO network analysis using Cytoscape 3.4.0 for DEGs over-expressed by ethephon. The colour of the nodes represents the significance of the over-representation, with a colour scale ranging from yellow at a p-value of 0.05 to dark orange at a p-value corresponding to 5 orders of magnitude smaller than the significance level (10⁻⁵ × 0.05).</figcaption></figure>

than laticifers might maintain a certain level of metabolic activity and gene expression, or the greater stringency in this study led to fewer DEGs being selected.

Ethephon is applied to the tapping cut of certain rubber clones to stimulate latex metabolism and prolong latex flow, and consequently latex production³⁰. Rubber clone PB 260, studied in this paper, has an intrinsic high latex metabolism⁶. For this clone, ethephon application can accelerate TPD occurrence: more than 60% of ethephon-treated trees were severely affected by TPD after 15 months of stimulation, while a small number of trees had a partial dry cut without ethephon treatment⁶. Several reports have revealed the dramatic change in gene expression in latex from ethephon-treated trees using RNA sequencing (8111 DEGs for clone PB 260 (this study), 10216 DEGs for clone PR 107²⁴, 509 DEGs for clone Reyan7-33-97²³), and cDNA microarrays (163 DEGs³¹). A proteomics analysis of ethylene-stimulated rubber latex identified 143 and 404 ethylene responsive latex proteins using mass spectrometry from 2-DE and DIGE gels, and iTRAQ, respectively³². The present study actually revealed only 3,883 over-expressed DEGs, mostly involved in “cellular process”, “metabolic process”, “binding” and “catalytic activity”. Wang and co-workers concluded that ethylene improves carbohydrate catabolism and energy production but not rubber biosynthesis itself directly³², confirming rapid acceleration of the glycolytic pathway supplying precursors for the biosynthesis of IPP, as suggested by Liu and co-workers too²⁴.

Ethylene biosynthesis and signalling pathways are particularly affected by ethephon stimulation. Following the first isolation of ethylene biosynthesis genes³³,³⁴, a genome-wide analysis showed that the *Hevea* genome activated 8 genes encoding S-adenosyl-L-methionine synthase (SAMS), 14 genes encoding 1-aminocyclopropane-1-carboxylic acid (ACC) synthase (ACS) and 16 genes encoding ACC oxidase (ACO)²⁵. Although 7 ACS and 8 ACO contigs were identified in the present RNA sequencing analysis, a small number of counts confirmed previous assumptions of low ethylene production in latex²³,³³. By contrast, latex cells are well-equipped for ethylene perception and signalling. A large number of ethylene receptors were expressed: 7 contigs in clone PB 260 (present work, Suppl. Data 4) and 4 genes in clone Reyan7-33-97²³. The reference genome has 181 Ethylene Response
<figure><img src="image_7.png" /><figcaption>**Figure 4.** BINGO network analysis using Cytoscape 3.4.0 for DEGs under-expressed by ethephon. The colour of the nodes represents the significance of the over-representation, with a colour scale ranging from yellow at a p-value of 0.05 to dark orange at a p-value corresponding to 5 orders of magnitude smaller than the significance level (10⁻⁵ × 0.05).</figcaption></figure>

Factors (ERF)²³, while only 91 have been classified and functionally verified²¹. For rubber clone PB 260, 27 and 39 DEGs encoding ERF were identified in latex by Piyatrakul and co-workers²¹ and in the present study, respectively. In clone Reyan7-33-97, 21 DEGs were reported²³. These results suggest a regulation of ethylene responsive genes in latex by ERFs in these rubber clones. Several members of three ERF groups (I, VII, VIII) are highly expressed in latex²¹. HbERF-II and HbERF-VIIIa deduced proteins contained the repressor EAR motif. Their high expression under normal conditions suggests a certain negative control of ethylene response. By contrast, their lowering in TPD-affected trees is expected to activate ethylene responsive genes. This deregulation could be the source of latex cell dysfunction. Several members of HbERF-IX are assumed to be essential integrators of complex hormonal signalling pathways in *Hevea*. HbERF-IXc4 and HbERF-IXc5 are two orthologs of ERF1, which is at the cross-talk of the ethylene and jasmonate signalling pathways³⁵. Lastly, ethylene stimulation induced several families of microRNAs³⁶⁻³⁸. This also suggests post-transcriptional regulation by ethylene.

This study also revealed dramatic down-regulation of DEGs related to jasmonic acid (JA) stimulus, carbohydrate and wounding responses. Tapping induces mechanical wounding, which is responsible for jasmonate
<figure><img src="image_8.png" /><figcaption>**Figure 5.** BINGO network analysis using Cytoscape 3.4.0 for DEGs over-expressed in severe TPD. The colour of the nodes represents the significance of the over-representation, with a colour scale ranging from yellow at a p-value of 0.05 to dark orange at a p-value corresponding to 5 orders of magnitude smaller than the significance level (10⁻⁵ × 0.05).</figcaption></figure>

production. Several upstream regulators down-expressed during TPD occurrence may explain these results. For instance, MYC transcription factors involved in the activation of JA-responsive genes were strongly under-expressed by ethephon³⁹. Regulation of MYC gene expression by ethephon was confirmed by this study as well as in TPD-affected trees. Several ethylene signalling transcription factors are also known to drive the response to jasmonate. The EIN3 gene, which encodes a primary transcription factor of the ethylene signalling pathway, is induced by both methyl jasmonate and ethylene, and might control JA/ethylene-responsive genes³⁶,⁴⁰. ERFs from groups II and VIII can also more specifically control JA-responsive genes. Contig CL17512Contig1 corresponding to HbERF-IIb2, an ortholog of ORA57 involved in the control of JA biosynthesis³⁵, was expressed at a very low level. This confirms previous assumptions that the expression of most genes related to jasmonate synthesis is repressed in TPD-affected trees¹⁸.

The involvement of oxidative stress in TPD was recently reviewed². This study showed that ethephon stimulation and TPD severity respectively induced the under-expression of 113 and 34 genes, which are related to "oxido/reductase activity" (Suppl. Data 8). Other transcriptomic and proteomic analyses also highlighted that ROS metabolism, the ubiquitin-proteasome pathway and programmed cell death were associated with TPD occurrence²⁹,⁴¹. The involvement of reactive oxygen species in the latex metabolism and TPD syndrome, and the wounding tolerance of rubber trees, suggest that some HbERFs from group VII, involved in the response to hypoxia, might play a major role in latex production. Several members of three ERF groups (I, VII, VIII) are highly expressed in latex²¹. In TPD-affected trees, the consumption of oxygen by NADH-cytochrome-c-oxidoreductase from lutoids is particularly high⁴². According to these authors, NADH-dependent oxygen consumption is inhibited by superoxide dismutase activity⁴³. ERF-VII transcription factors are known to regulate hypoxia-responsive genes. HbERF-VIIa17, an ortholog gene of AtEBP, might play a role in the response to the accumulation of reactive oxygen species occurring during latex regeneration²¹.

To conclude, this study is the first high-throughput RNA sequencing in latex at the onset of ethephon-induced TPD. The overall differentially expressed genes associated with ethephon stimulation and the severity of TPD have been identified and a gene network has been proposed to address this issue. The depression in jasmonate
<figure><img src="image_9.png" /><figcaption>**Figure 6.** BINGO network analysis using Cytoscape 3.4.0 for DEGs under-expressed in severe TPD. The colour of the nodes represents the significance of the over-representation, with a colour scale ranging from yellow at a p-value of 0.05 to dark orange at a p-value corresponding to 5 orders of magnitude smaller than the significance level (10⁻⁵ × 0.05).</figcaption></figure>

<table><caption>**Table 4.** Main transcription factor families regulated in response to ethephon and severity of ethephon-induced TPD.</caption><thead><tr><th rowspan="2">Transcription factor family</th><th colspan="3">Ethephon vs control</th><th colspan="3">Severe vs slight TPD</th></tr><tr><th>Total</th><th>Up</th><th>Down</th><th>Total</th><th>Up</th><th>Down</th></tr></thead><tbody><tr><td>ERF</td><td>39</td><td>23</td><td>16</td><td>3</td><td>1</td><td>2</td></tr><tr><td>ARF</td><td>33</td><td>9</td><td>24</td><td>3</td><td>0</td><td>3</td></tr><tr><td>WRKY</td><td>31</td><td>16</td><td>15</td><td>4</td><td>4</td><td>0</td></tr><tr><td>Homeobox</td><td>31</td><td>7</td><td>24</td><td>5</td><td>0</td><td>5</td></tr><tr><td>bHLH</td><td>29</td><td>13</td><td>16</td><td>2</td><td>0</td><td>2</td></tr><tr><td>MYB</td><td>21</td><td>10</td><td>11</td><td>6</td><td>1</td><td>5</td></tr><tr><td>MADS box</td><td>12</td><td>4</td><td>8</td><td>1</td><td>0</td><td>1</td></tr><tr><td>TCP</td><td>12</td><td>5</td><td>7</td><td>2</td><td>0</td><td>2</td></tr><tr><td>GATA</td><td>11</td><td>5</td><td>6</td><td>1</td><td>0</td><td>1</td></tr><tr><td>MYC</td><td>9</td><td>4</td><td>5</td><td>1</td><td>0</td><td>1</td></tr><tr><td>HSF</td><td>7</td><td>4</td><td>3</td><td>2</td><td>0</td><td>2</td></tr><tr><td>NAC</td><td>4</td><td>1</td><td>3</td><td>1</td><td>0</td><td>1</td></tr><tr><td>GRAS</td><td>3</td><td>0</td><td>3</td><td>2</td><td>2</td><td>0</td></tr><tr><td>Other</td><td>74</td><td>32</td><td>42</td><td>17</td><td>3</td><td>14</td></tr><tr><td>Total</td><td>316</td><td>133</td><td>183</td><td>50</td><td>11</td><td>39</td></tr></tbody></table>

response and ROS-scavenging systems were a common feature observed in latex (this study) and in bark (previous studies). Allelic variation in the sequences of gene expression markers associated with TPD might be exploited on already mapped and phenotyped segregating populations in order to identify the genetic markers associated with phenotypes related to latex production and TPD susceptibility.
**Methods**

**Plant material.** Latex was collected from 9-year-old trees of clone PB 260 at the Sembawa Research Centre belonging to the Indonesian Rubber Research Institute (Palembang, Sumatra, Indonesia). Trees were tapped for three years using the S/2 d2 harvesting system (semi-spiral and tapping every two days) with or without application of 2.5% ethephon 12 times a year (ET 12/y). The dry cut length (DCL) was monitored after tapping for 3 years. Trees were categorized as slight- (<25% of DCL) or severe-TPD trees (>75%) according to the mean DCL for the last year of tapping. Nine trees were selected for further RNA sequencing: three slight-TPD trees under S/2 d2 (control), three slight-TPD trees under S/2 d2 ET 12/y, and three severe-TPD trees under S/2 d2 ET 12/y.

**Total RNA isolation.** For each selected tree, 6 mL of latex was collected in a 15-mL Falcon tube containing 6 ml of 2X alkaline fixing buffer (0.1 M Tris-HCl, 0.3 M LiCl, 10 mM EDTA, 10% SDS pH 9) with constant manual mixing, and immediately deep-frozen in liquid nitrogen. Total RNAs were isolated using the caesium chloride (CsCl) cushion method adapted from Sambrook⁴⁴ by Duan and coll.³⁶. Latex samples were thawed in a water bath at 50 °C and transferred to centrifuge tubes. After centrifugation for 30 min at 15,000 g at 20 °C, the aqueous fraction was transferred to a 13 mL-tube, then the volume was adjusted to 10 mL with extraction buffer containing 4 M guanidium isothiocyanate, 1% sarcosine, 1% polyvinylpyrrolidone and 1% β-mercapto-ethanol. After homogenization, the tubes were kept on ice, then centrifuged at 10,000 g at 4 °C for 30 minutes. The supernatant was transferred to a new tube containing 3 mL of 5.7 M CsCl. Ultracentrifugation in a swinging bucket was carried out at 82,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. Although DNA could not cross the caesium cushion for this centrifugation condition, a DNAse treatment was performed using TurboDNAse (Ambion, Life Technologies, Carlsbad, USA). DNA contamination was then checked by PCR amplification using primers of the *Actin* gene including the intron sequence. RNAs were conserved at -80 °C pending further analysis.

**RNA sequencing and count analysis.** The quality and quantity of total RNA samples were checked using the Agilent Bioanalyzer and Quant-iT RiboGreen RNA Assay Kit, respectively. The sequencing technology used was an Illumina Hiseq 2000 (thanks to IG-CNS for allowing us privileged access to carry out sequencing). RNA-seq libraries were produced by the TruSeq Stranded mRNA SamplePrep Guide_15031047_D protocol (Illumina®, California, U.S.A.) using 2 µg of total RNA. The RNA-seq samples were sequenced in paired-end (PE) mode with a sizing of 260 bp and a read length of 100 bases. Seven samples per Hiseq 2000 lane were analysed using individual bar-coded adapters, generating approximately 60 million PE reads per sample.

Each RNA-seq sample followed the same pipeline of bioinformatics analyses from trimming to the counting of transcript abundance as follows. Read pre-processing criteria included trimming library adapters and performing quality control. The raw data (fastq) were trimmed by the fastx toolkit for Phred Quality Score Qscore >20, read length >30 bases, and the ribosome sequences were removed with the sortMeRNA tool⁴⁵. The mapper Bowtie version 2⁴⁶ was used to align reads against the rubber clone PB 260 reference transcriptome²¹ (with-local option and other default parameters). Reads were counted for the 86,941 contigs. The abundance of each contig was calculated by a local script which parses SAM files and counts only paired-end reads for which both reads map unambiguously one contig, removing multi-hits. According to this procedure, around 95% of PE reads were associated with a contig.

Statistical analyses were carried out with R software (version 3.3.1) and the edgeR package (version1.12.0). Library sizes were normalized using the TMM method. A differential expression analysis was performed using an exact test based on negative binomial distribution using an empirical Bayes estimation of the dispersion, where data were filtered by the default mode proposed by edgeR: genes that did not have at least 1 read after a count per million in at least half of the samples were discarded. To control the FDR, raw p-values were adjusted. We considered as being differentially expressed, those contigs with an adjusted p-value under 0.05, fulfilling criteria for a robust RNA-Seq analysis⁴⁷. All steps of the transcriptomic experiment, from growth conditions to bioinformatic analyses were stored and managed through the CATdb database⁴⁸ (http://tools.ips2.u-psud.fr/fr/CATdb/) Project_ID NGS2014_07_Hevea according to the “minimum information about a high-throughput sequencing experiment” standards.

A principal component analysis (PCA) and heatmap clustering analysis were performed using the R Mixomics and flashClust packages, respectively, after log transformation of the normalized count data. Data scaling and the Euclidean distance were calculated by the Pareto and Ward methods, respectively.

A BINGO network analysis was carried out using Cytoscape 3.4.0 for the DEG groups. DEGs were converted into their respective *Arabidopsis thaliana* orthologs selected according to the best e-values. The At DEGs were implemented to perform gene ontology (GO) enrichment using the BiNGO plugin (Cytoscape® software)⁴⁹ following the procedure by default. Briefly, over-representation of the biological process was assessed when compared to *Arabidopsis thaliana* whole annotation as the reference. A hypergeometric test was used as a statistical test, then multiple testing correction was applied using Benjamini & Hochberg False Discovery Rate (FDR) correction at the significant level of 0.05. The Bingo graph represents the over-represented GO category as nodes. Node size is proportional to the gene number related to the respective category. The colour of the nodes represents the significance of the over-representation, with a colour scale ranging from yellow at a p-value of 0.05 to dark orange at a p-value corresponding to 5 orders of magnitude smaller than the significance level (10⁻⁵ × 0.05).

**Data deposition.** The RNA sequencing data have been submitted to the international repository GEO (Gene Expression Omnibus, Edgard R. et al. 2002, http://www.ncbi.nlm.nih.gov/geo) ProjectID GSE101568.
## References

1. Okoma, K. M., Dian, K., Obouayeba, S., Elabo, A. A. E. & N'guetta, A. S. P. Seasonal variation of tapping panel dryness expression in rubber tree *Hevea brasiliensis* Muell.arg. in Côte d'Ivoire. *Agriculture and Biology Journal of North America* **2**, 559-569 (2011).

2. Zhang, Y., Leclercq, J. & Montoro, P. Reactive oxygen species in Hevea brasiliensis latex and relevance to Tapping Panel Dryness. *Tree Physiol*, https://doi.org/10.1093/treephys/tpw106 (2016).

3. Sookmark, U. et al. Characterization of polypeptides accumulated in the latex cytosol of rubber trees affected by the tapping panel dryness syndrome. *Plant Cell Physiol* **43**, 1323-1333 (2002).

4. de Faÿ, E. Histo- and cytopathology of trunk phloem necrosis, a form of rubber tree (Hevea brasiliensis Müll. Arg.) tapping panel dryness. *Australian Journal of Botany* **59**, 563-574 (2011).

5. de Faÿ, E., Moraes, L. A. C. & Moraes, V. H. F. Cyanogenesis and the onset of tapping panel dryness in rubber tree. *Pesquisa Agropecuária Brasileira* **45**, 1372-1380 (2010).

6. Putranto, R. A. et al. Involvement of Ethylene in the Latex Metabolism and Tapping Panel Dryness of Hevea brasiliensis. *International journal of molecular sciences* **16**, 17885-17908, https://doi.org/10.3390/ijms160817885 (2015).

7. Sobhana, P., Thomas, M. & Jacob, J. Genetic distance between the root stock and scion of healthy and TPD affected trees of Hevea brasiliensis. *Natural Rubber Research* **18**, 101-104 (2005).

8. Nandris, D., Moreau, R., Pellegrin, F. & Chrestin, C. Rubber Tree (Hevea brasiliensis) Bark Necrosis Syndrome: 2. First Comprehensive Report on Causal Stresses. *Plant Dis.* **88**, 1047 (2004).

9. Chen, S. et al. Association of decreased expression of a Myb transcription factor with the TPD (tapping panel dryness) syndrome in Hevea brasiliensis. *Plant Mol Biol* **51**, 51-58 (2002).

10. Venkatachalam, P., Thulaseedharan, A. & Raghothama, K. Identification of expression profiles of tapping panel dryness (TPD) associated genes from the latex of rubber tree (Hevea brasiliensis Muell. Arg.). *Planta* **226**, 499-515, https://doi.org/10.1007/s00425-007-0500-8 (2007).

11. Venkatachalam, P., Thulaseedharan, A. & Raghothama, K. Molecular identification and characterization of a gene associated with the onset of tapping panel dryness (TPD) syndrome in rubber tree (Hevea brasiliensis Muell.) by mRNA differential display. *Mol Biotechnol* **41**, 42-52, https://doi.org/10.1007/s12033-008-9095-y (2009).

12. Venkatachalam, P., Geetha, N., Priya, P. & Thulaseedharan, A. Identification of a differentially expressed thymidine kinase gene related to tapping panel dryness syndrome in the rubber tree (Hevea brasiliensis Muell. Arg.) by random amplified polymorphic DNA screening. *Int. J. Plant Biology* **1**, e7 (2010).

13. Venkatachalam, P. Molecular cloning and characterization of a novel down-regulated TOM20 gene associated with the onset of tapping panel dryness (TPD) syndrome from the inner bark tissues of rubber tree (Hevea brasiliensis Muell. Arg.) by mRNA differential display. *Abstr Pap Am Chem S* 240 (2010).

14. Qin, B., Liu, X.-H., Deng, Z. & Li, D.-J. Identification of genes associated with tapping panel dryness in Hevea brasiliensis using oligonucleotide microarrays. *Chin J Trop Crops* **33**, 296-301 (2012).

15. Li, D., Qin, B., Deng, Z. & Liu, X.-H. Identifying the genes related to Tapping by Hevea brasiliensis oligonucleotide microarrays. *Plant Physiology Journal* **48**, 901-908 (2012).

16. Deng, Z. et al. Expression Profiles, Characterization and Function of HbTCTP in Rubber Tree (Hevea brasiliensis). *Frontiers in plant science* **7**, 789 (2016).

17. Li, D. et al. Transcriptome analyses reveal molecular mechanism underlying tapping panel dryness of rubber tree (Hevea brasiliensis). *Sci Rep* **6**, 23540, https://doi.org/10.1038/srep23540 (2016).

18. Liu, J. P., Xia, Z. Q., Tian, X. Y. & Li, Y. J. Transcriptome sequencing and analysis of rubber tree (Hevea brasiliensis Muell.) to discover putative genes associated with tapping panel dryness (TPD). *BMC Genomics* **16**, 398, https://doi.org/10.1186/s12864-015-1562-9 (2015).

19. Gébelin, V. et al. The small RNA profile in latex from Hevea brasiliensis trees is affected by tapping panel dryness. *Tree Physiology*, https://doi.org/10.1093/treephys/tpt076 (2013).

20. Yuan, K., Zhou, X., Wang, Z. & Yang, L. Identification and analysis of latex rubber particle proteins related with tapping panel dryness in Hevea brasiliensis. *Journal of Nanjing Forestry University (Natural Sciences Edition)* **38**, 36-40 (2014).

21. Piyatrakul, P. et al. Sequence and Expression Analyses of Ethylene Response Factors Highly Expressed in Latex Cells from Hevea brasiliensis. *PLoS One* **9**, e99367, https://doi.org/10.1371/journal.pone.0099367 (2014).

22. St Laurent, G. et al. On the importance of small changes in RNA expression. *Methods* **63**, 18-24, https://doi.org/10.1016/j.ymeth.2013.03.027 (2013).

23. Tang, C. et al. The rubber tree genome reveals new insights into rubber production and species adaptation. *Nature Plants* **2**, 16073, https://doi.org/10.1038/nplants.2016.73 (2016).

24. Liu, J. P., Zhuang, Y. F., Guo, X. L. & Li, Y. J. Molecular mechanism of ethylene stimulation of latex yield in rubber tree (Hevea brasiliensis) revealed by de novo sequencing and transcriptome analysis. *BMC Genomics* **17**, 257, https://doi.org/10.1186/s12864-016-2587-4 (2016).

25. Li, H.-L. et al. Comparative Transcriptome Analysis of Latex Reveals Molecular Mechanisms Underlying Increased Rubber Yield in Hevea brasiliensis Self-Rooting Juvenile Clones. *Frontiers Plant Science* **7**, 1204 (2016).

26. Wei, F. et al. Transcriptome sequencing and comparative analysis reveal long-term flowing mechanisms in Hevea brasiliensis latex. *Gene* **556**, 153-162, https://doi.org/10.1016/j.gene.2014.11.048 (2015).

27. Li, D. et al. Next-generation sequencing, assembly, and comparative analyses of the latex transcriptomes from two elite *Hevea brasiliensis* varieties. *Tree Genetics & Genomes* **11**, 1-9 (2015).

28. Duan, C. et al. Identification of the *Hevea brasiliensis* AP2/ERF superfamily by RNA sequencing. *BMC Genomics* **14**, 30, https://doi.org/10.1186/1471-2164-14-30 (2013).

29. Li, D. et al. Identification and characterization of genes associated with tapping panel dryness from Hevea brasiliensis latex using suppression subtractive hybridization. *BMC Plant Biol* **10**, 140, https://doi.org/10.1186/1471-2229-10-140 (2010).

30. Tungngoen, K. et al. Hormonal treatment of the bark of rubber trees (Hevea brasiliensis) increases latex yield through latex dilution in relation with the differential expression of two aquaporin genes. *J Plant Physiol* **168**, 253-262 (2011).

31. Nie, Z. et al. Profiling Ethylene-Responsive Genes Expressed in the Latex of the Mature Virgin Rubber Trees Using cDNA Microarray. *PLoS One* **11**, e0152039, https://doi.org/10.1371/journal.pone.0152039 (2016).

32. Wang, X. et al. Comprehensive Proteomics Analysis of Laticifer Latex Reveals New Insights into Ethylene Stimulation of Natural Rubber Production. *Sci Rep* **5**, 13778, https://doi.org/10.1038/srep13778 (2015).

33. Kuswanhadi et al. Isolation of Three Members of The Multigene Family Encoding ACC Oxidases in *Hevea brasiliensis* and Inverstigation of Their Responses to Ethylene Stimulation and Wounding. *Journal of Rubber Research* **13**, 185-205 (2010).

34. Zhu, J. H., Xu, J., Chang, W. J. & Zhang, Z. L. Isolation and Molecular Characterization of 1-Aminocyclopropane-1-carboxylic Acid Synthase Genes in Hevea brasiliensis. *International journal of molecular sciences* **16**, 4136-4149, https://doi.org/10.3390/ijms16024136 (2015).

35. Putranto, R. A. et al. Ethylene Response Factors Are Controlled by Multiple Harvesting Stresses in Hevea brasiliensis. *PLoS One* **10**, e0123618, https://doi.org/10.1371/journal.pone.0123618 (2015).

36. Duan, C. et al. Gene expression pattern in response to wounding, methyl jasmonate and ethylene in the bark of Hevea brasiliensis. *Tree Physiol* **30**, 1349-1359 (2010).
37. Gébelin, V., Leclercq, J., Hu, S., Tang, C. & Montoro, P. Regulation of MIR genes in response to abiotic stress in Hevea brasiliensis. *International journal of molecular sciences* **14**, 19587–19604, https://doi.org/10.3390/ijms141019587 (2013).

38. Pramoolkit, P. et al. Involvement of ethylene-responsive microRNAs and their targets in increased latex yield in the rubber tree in response to ethylene treatment. *Plant Physiol Biochem* **84**, 203–212, https://doi.org/10.1016/j.plaphy.2014.09.016 (2014).

39. Pirrello, J. et al. Transcriptional and post-transcriptional regulation of the jasmonate signalling pathway in response to abiotic and harvesting stress in *Hevea brasiliensis*. *BMC Plant Biology* **14**, 17 p (2014).

40. Yang, Z.-P., Li, H.-L., Guo, D. & Peng, S.-Q. Identification and characterization of the EIN3/EIL gene family in Hevea brasiliensis. *Tree Genetics and Genomes* **11**, 11 p (2015).

41. Yuan, K., Wang, Z. H., Zhou, X. M., Zou, Z. & Yang, L. F. The identification of differentially expressed latex proteins in healthy and tapping panel dryness (TPD) Hevea brasiliensis trees by iTRAQ and 2D LC-MS/MS. *Acta Agriculturae Universitatis Jiangxiensis* **36**, 650–655 (2014).

42. Jacob, J.-L. & d'Auzac, J. In *Physiology of Rubber Tree Latex* (eds d'Auzac, J. Jacob, J. L. & Chrestin, H.) 59–96 (CRC press, Inc., 1989).

43. Chrestin, H. In *Physiology of rubber tree latex* (eds d'Auzac, J. Jacob, J. L. & Chrestin, H.) 431–458 (CRC Press Inc., 1989).

44. Sambrook, J., Fritsch, E. F. & Maniatis, T. *Molecular cloning, a laboratory manual*. (CHS Press, 1989).

45. Kopylova, E., Noe, L. & Touzet, H. SortMeRNA: fast and accurate filtering of ribosomal RNAs in metatranscriptomic data. *Bioinformatics* **28**, 3211–3217, https://doi.org/10.1093/bioinformatics/bts611 (2012).

46. Langmead, B. & Salzberg, S. L. Fast gapped-read alignment with Bowtie 2. *Nat Methods* **9**, 357–359, https://doi.org/10.1038/nnmeth.1923 (2012).

47. Zhang, Z. H. et al. A comparative study of techniques for differential expression analysis on RNA-Seq data. *PLoS One* **9**, e103207, https://doi.org/10.1371/journal.pone.0103207 (2014).

48. Gagnot, S. et al. CATdb: a public access to Arabidopsis transcriptome data from the URGV-CATMA platform. *Nucleic Acids Res* **36**, D986–990, https://doi.org/10.1093/nar/gkm757 (2008).

49. Kucera, M., Isserlin, R., Arkhangorodsky, A. & Bader, G. D. AutoAnnotate: A Cytoscape app for summarizing networks with semantic annotations. *F1000Research* **5**, 1717, https://doi.org/10.12688/f1000research.9090.1 (2016).

## Acknowledgements

This work was supported by CIRAD and IRRI (Indonesian Rubber Research Institute). The POPS platform benefited from the support of the LabEx Saclay Plant Sciences-SPS (ANR-10-LABX-0040-SPS). We are grateful to Fina Bintardi and Painah for technical help, and Ludivine Soubigou-Taconnat for fruitful discussion. The authors also thank Peter Biggins for his revision of the English.

## Author Contributions

P.M. coordinated the project, supervised data generation and analyses, and drafted the paper. K. supervised the field experiment. P.M. and K. supervised latex collection. E.H. performed the field experiment and dry cut length analysis. P.M., K., E.H. and S.I. collected latex samples. M.R. carried out RNA isolation. C.L. carried out RNA sequencing. M.L.M. carried out mapping and edgeR analysis. P.M., S.W., B.F., S.P. and J.L. analysed the data.

## Additional Information

Supplementary information accompanies this paper at https://doi.org/10.1038/s41598-018-26854-y.

**Competing Interests:** The authors declare no competing interests.

**Publisher's note:** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_10.png" /></figure>

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2018