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

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.


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 21 . 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 22 , 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 TPDfrom 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 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   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 23 . 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 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 21,[24][25][26][27][28] . 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 [9][10][11]15,29 . 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 17,18 . 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) 17,18 . The down-regulation of some genes related to rubber biosynthesis in bark 17,18 was confirmed at protein level 20 . 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  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 30 . Rubber clone PB 260, studied in this paper, has an intrinsic high latex metabolism 6 . 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 6 . 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 24 , 509 DEGs for clone Reyan7-33-97 23 ), and cDNA microarrays (163 DEGs 31 ). 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 32 . 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 32 , confirming rapid acceleration of the glycolytic pathway supplying precursors for the biosynthesis of IPP, as suggested by Liu and co-workers too 24 .
Ethylene biosynthesis and signalling pathways are particularly affected by ethephon stimulation. Following the first isolation of ethylene biosynthesis genes 33,34 , 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) 23 . 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 23,33 . 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 23 . The reference genome has 181 Ethylene Response DEGs encoding ERF were identified in latex by Piyatrakul and co-workers 21 and in the present study, respectively. In clone Reyan7-33-97, 21 DEGs were reported 23 . 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 21 . 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 crosstalk of the ethylene and jasmonate signalling pathways 35 . Lastly, ethylene stimulation induced several families of microRNAs [36][37][38] . 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 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 39 . 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 36,40 . 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 35 , 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 18 .
The involvement of oxidative stress in TPD was recently reviewed 2 . 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 29,41 . 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 21 . In TPD-affected trees, the consumption of oxygen by NADH-cytochrome-c-oxidoreductase from lutoids is particularly high 42 . According to these authors, NADH-dependent oxygen consumption is inhibited by superoxide dismutase activity 43 . 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 21 .
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 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.   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 45 . The mapper Bowtie version 2 46 was used to align reads against the rubber clone PB 260 reference transcriptome 21 (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 47 . All steps of the transcriptomic experiment, from growth conditions to bioinformatic analyses were stored and managed through the CATdb database 48 (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) 49 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 −5 × 0.05).