Transcript Proling Provides Insights into the Molecular Mechanisms of Harvesting-activated Latex Regeneration in Virgin Rubber Trees

Background: Natural rubber, an important industrial raw material with wide applications, is harvested in the form of latex (cytoplasm of rubber-producing laticifers) from Hevea brasiliensis (para rubber tree) by the way of tapping, i.e. removing a slice of trunk bark by a special knife. In regularly tapped rubber trees, latex regeneration consists of one of the main yield-limiting factors for rubber productivity. Conspicuous stimulation on latex production for the rst few tappings makes virgin (untapped before) rubber trees an ideal model to investigate the regulatory mechanisms of latex regeneration. To understand the underlying mechanisms, genome-wide transcript proling was conducted with a silver-staining cDNA-AFLP technology against the latex samples for the rst ve tappings. Results: A total of 505 non-redundant differentially expressed (DE) transcript-derived fragments (TDFs) were identied, of which 217 were up-regulated, 180 down-regulated, and 108 bell type-regulated among the ve tappings. About 72.5% of these DE-TDFs were functionally annotated, and classied into 11 functional categories, which were discussed with reference to harvesting-stimulated latex regeneration. The importance of sugar metabolism and rubber biosynthesis was highlighted, due to the fact that most of the DE-TDFs annotated in sucrose transport, sugar catabolism, glycolysis, tricarboxylic acid cycle and pentose-phosphate pathway and nine of the ten rubber biosynthesis pathway DE-TDFs were up-regulated by the tapping treatment. More than one tenth of the total DE-TDFs were randomly selected for expression validation by semi-quantitative RT-PCR, and 83.8% showed patterns consistent with their original cDNA-AFLP gel proles. Moreover, quantitative RT-PCR analysis revealed an 89.7% consistency for the 29 latex-regeneration related DE-TDFs examined. Conclusions: In brief, our results indicate the tapping treatment incurs extensive physiological and molecular changes in the laticifers of virgin rubber trees. The vast numbers of tapping-responsive DE-TDFs identied here provide a basis for unravelling the gene regulatory network for latex regeneration in regularly harvested rubber trees. strand cDNA diluted qRT-PCR with HbYLS8 includes: template, μl for μM forward and reverse primers, 10 μl 2×SYBR ® Premix Ex Taq TM (cid:0) (Takara) and 7.4 μl ddH 2 O. Roche's LightCycler 2.0 system was used for qRT-PCR analysis with the program as follows: 95 ℃ 30 sec; 94 ℃ 5 sec, 60 ℃ 20 sec, 72 ℃ 20 sec, 45cycles. Three technical replicates were analyzed for each of the three biological samples. All the cycle threshold (Ct) values from one gene were determined at the same threshold uorescence value of 0.2 using the 2 -ΔΔCt method. The primers of target and reference genes were listed in Table S4. Statistical analysis was performed using Student’s t-test.


Introduction
Natural rubber (cis-1, 4-polyisoprene, NR) is an elastomer with superior properties that cannot be completely replaced by petroleum-derived synthetic rubber, and is used as an important industrial raw material related to national economy and people's livelihood. There are more than 2, 500 kinds of NRbearing plants. However, due to the advantages of good quality, high yield, and easiness for harvesting, Hevea brasiliensis (para rubber tree, Hevea thereafter) has become the sole commercial NR source [1]. Hevea trees need high temperature and humid climate conditions for normal growth and NR production yet vulnerable to typhoon, thus con ning its planting to restricted tropical areas [2].
Many factors are affecting Hevea rubber yield. From the physiological perspective, there are three main factors: duration of latex ow after tapping, the capability of latex regeneration between two consecutive tappings, and the ability of laticifer differentiation in bark cambium [3,4]. The number of laticifer rings in tapped Hevea bark is 1 to 3 times more than un-tapped, and tapping promotes laticifer differentiation [4]. Meanwhile, mechanical injury and jasmonic acid stimulation facilitate laticifer differentiation, latex regeneration and production [6,7]. Latex is the cytoplasm that ows out of the laticifers after tapping, of which dry weight more than 90% is rubber hydrocarbon, i.e. NR [8]. In regularly tapped Hevea trees, the expelled latex can be regenerated within 2 to 3 days, latex regeneration thus representing the main metabolic activity occurred in the laticifers of trunk bark where harvesting is conducted.
After tapping, the substance owing out of the laticifers includes rubber hydrocarbon, organelles, sugars, organic acids, nucleic acids and proteins/enzymes [3]. Tupý reported a positive correlation between sucrose content and latex yield [9,10,11]. A sucrose transporter, HbSUT3, is responsible for sucrose loading into laticifers, and its expression has been found to be induced by the treatments of ethylene and tapping, both bolstering the latex yield [12]. A kind of alkaline/neutral invertase cleaves sucrose in the latex into glucose and fructose that are then exploited in subsequent latex production [11]. Liu et al.
identi ed the responsible invertase gene, HbNIN2, and found its expression and enzymatic activity in the latex are positively correlated with latex yield [13]. NR biosynthesis occurs on rubber particles, a kind of organelle where the proteins of rubber elongation factor (REF)/small rubber particle protein (SRPP) are important components of the rubber biosynthesis machinery [14][15][16][17]. Priya et al. observed a positive correlation between REF mRNA abundance and the yielding levels of Hevea clones [18]. Amalou et al. revealed that a marked increase in transtonoplast ΔpH within Hevea laticifers, which consists of one of the major mechanisms of ethylene stimulation on latex yield [19]. Therefore, latex regeneration, with the biosynthesis of rubber hydrocarbon as a major activity, involves a complex regulatory network of gene expression, multi-enzyme reaction and physio-biochemical processes. Virgin Hevea trees produce very little latex at the rst tapping, and the latex yield increases signi cantly with subsequent tappings with regular intervals of 2 to 3 days, and reaches a relatively stable level after 7-10 tappings [12,20]. Adiwilaga revealed that in the latex the accumulation of farnesyl diphosphate synthase transcripts is induced by tapping [21]. Tang et al. concluded that in virgin Hevea trees the metabolic activity of laticifers is mediocre, and gets activated with the tapping treatment and maintains a high level of dynamic equilibrium for latex regeneration in regularly tapped trees [12]. Therefore, virgin Hevea trees could be an ideal material to study the mechanisms of latex regeneration and to identify the genes involved.
Genome-wide transcript pro ling techniques play an important role in large-scale screening and characterization of the genes involved in various biological processes. Since the advent of the cDNA-AFLP transcript-pro ling technique in 1996 [22], owing to its advantages of good repeatability, high sensitivity and high throughput, this technique has been successfully applied to the eld of plant biology, such as abiotic stress response [23], plant-microorganism interaction [24], hormone signal [25] and development [26].
In this study, the pro les of latex transcriptome in virgin Hevea trees were systematically compared for the rst ve tappings using a silver staining cDNA-AFLP protocol we previously established for Hevea latex transcriptome [27]. A total of 505 non-redundant TDFs were identi ed as differentially expressed (DE) with the ve tappings. These DE-TDFs were further veri ed for their expression patterns by semi-quantitative RT-PCR (sqRT-PCR) and quantitative RT-PCR (qRT-PCR). The results highlight the importance of sucrose transport and sugar metabolism as well as rubber biosynthesis in tapping-activated latex regeneration in virgin Hevea trees.

Screening of DE-TDFs
To extensively identify the differentially expressed TDFs (DE-TDFs) responding to the tapping treatment, all the 128 selective primer combinations of Apo /Mse enzyme restriction system were screened in the latex for the rst ve tappings in virgin Hevea trees (Table 1). On average, about 70 TDFs greater than 100 bp were discernable on the silver-stained polyacrylamide gels for each pair of primer combination ( Fig. 1). Therefore, nearly 9, 000 TDFs were expected to be pro led for each latex RNA sample when all the 128 Apo /Mse selective primer pairs are examined. In total, 651 DE-TDFs were identi ed and successfully cloned and sequenced. Taking the gene expression level at the rst tapping as a reference, the DE-TDFs identi ed were classi ed into three types according to their patterns of expression along with the ve successive tappings: up-regulation, down-regulation and bell-type regulation (Fig. 1). The upregulation type includes three subtypes: (i) increase successively; (ii) increase rst and then stabilize; (iii) increase rst and then decrease, but still higher than the rst tapping. The down-regulation type also includes three subtypes: (i) decrease successively; (ii) decrease rst and then stabilizes; (iii) decrease rst and then increase, but still lower than the rst tapping. The bell-type regulation includes two subtypes: (i) increase rst, reaching a high threshold, and then decrease to a level lower than the rst tapping (upwardbell); (ii) decrease rst, reaching a low threshold, and then increase to a level higher than the rst tapping (downward-bell).

DE-TDFs annotation and redundancy removal
The DE-TDFs were made clean by removing the sequences of vector, primer and adaptor at both ends and then annotated by Blastx online searching (http://blast.ncbi.nlm.nih.gov/Blast.cgi) against the NCBI nonredundant protein database (nr), with a threshold of E value <10 -4 and score >50. According to the blasting against the Hevea latex transcriptome database [17], the DE-TDFs belonging to the same transcript and sharing similar expression pattern in the cDNA-AFLP analysis were clustered together, and only the longest TDF was retained. As a result, a total of 505 non-redundant DE-TDFs were obtained, including 217 (43.0%) up-regulated (Table S1), 180 (35.6%) down-regulated (Table S2) and 108 (21.4%) bell-regulated (Table S3).

Functional classi cation of non-redundant DE-TDFs
According to the results of Blastx searching, the 505 non-redundant DE-TDFs were divided into four categories: (i) Proteins with clear functional annotation; (ii) Unclassi ed proteins, the functional annotation of the proteins being multiple; (iii) Predicted protein, showing high homology with a predicted protein in the database; (iv) No hit, no homologous sequence in the database. Most (366, 72.5%) of these DE-TDFs had homology to genes with known functions, whereas 30 (5.9%) and 59 (11.7%), respectively, belonged to unclassi ed proteins and predicted proteins and the remaining 50 (9.9%) were with no hit ( Table 2).
With reference to the functional categories of plant genes de ned by Bevan et al. [28], the 366 DE-TDFs with known functions were classi ed into 11 functional categories, of which a new category, rubber biosynthesis, had been singled out from "secondary metabolism" (Fig. 2 and Table 2). Of these categories, ve, i.e. cell growth and division, protein degradation and storage, cellular structure, secondary metabolism, and rubber biosynthesis had a higher portion of up-regulated DE-TDFs than that of the down-and bell-regulated DE-TDFs together ( Table 2). Strikingly, nine of the ten DE-TDFs implicated in rubber biosynthesis were up-regulated, and the remaining one is upward bell-regulated.

Validation of expression pattern by sqRT-PCR
To determine the reliability of the cDNA-AFLP results, 80 DE-TDFs, covering >10% of the DE-TDFs we identi ed and sequenced, were randomly selected from each functional category and subjected to sqRT-PCR analysis using speci c primers for TDFs with 18S rRNA as the reference gene (Table 3). About 84% (67 TDFs) showed the expression pattern consistent with their cDNA-AFLP gel pro les, indicating the high reliability of the cDNA-AFLP screening.

qRT-PCR analysis of latex regeneration-related DE-TDFs
According to the metabolism pathways reported to be involved in latex regeneration, 29 latex-regeneration related DE-TDFs were further investigated by qRT-PCR analysis for their expression patterns across the ve successive tappings (Table S4). About 90% of the qRT-PCR results were consistent with their original cDNA-AFLP expression pro les ( Fig. 3; Table 4). The genes of these DE-TDFs are putatively involved in the pathways of primary metabolism, rubber biosynthesis and regulation, transporters and intracellular transport. Of the ten rubber biosynthesis pathway DE-TDFs, nine revealed qRT-PCR patterns similar to their cDNA-AFLP results (Table 4 and Fig. 3).

Discussion
Functional categories with reference to tapping-stimulated latex regeneration Tapping can stimulate the regeneration of latex, especially in virgin Hevea trees [12,20]. A number of early studies have shown that the rst few tappings greatly stimulate the metabolism of laticifers, accompanied by the enhanced expression of several speci c genes involved in latex regeneration [18,21,[29][30]. The latex ows out of laticifers after tapping, and in order to compensate for the loss of cytoplasm (latex) and maintain the balance of intracellular metabolism, the laticifers require large amounts of RNA and proteins to be synthesized before the next tapping. Of the 366 DE-TDFs identi ed with known function (Table 2), 26.2% were classi ed into the functional category of transcription and protein synthesis (Fig. 2), representing the largest category, 42.7% of which were up-regulated by the tapping treatment. These results indicated that tapping signi cantly affects the ways of laticifers to synthesize RNA and proteins, thus laying a foundation for their subsequent physiological response to the tapping treatment [3]. Laticifers are believed to be a defense system for Hevea trees to cope with biological and abiotic stresses, and the latex exuded after bark wounding has been found to play roles in resisting pathogen infection, insect feeding and abiotic stress [31]. The tapping itself is a kind of abiotic stress upon Hevea trees. Therefore, it is reasonable that "stress and defense" also accounted for a large portion of the functional DE-TDFs identi ed responsive to tapping ( Table 2; Fig. 2), ranking the second place in functional categories. The harvesting stress response has been suggested to be one of the key factors affecting latex regeneration and rubber productivity in Hevea trees [32].
The category of transporters and intracellular was the third largest among the 11 functional categories, accounting for 12.3% of the total functional DE-TDFs (Fig. 2). This corresponds well to the sink effect caused by the large loss of latex after tapping. The process of regenerating the expelled latex involves the synthesis, transport, loading and subcellular localization of a large number of organelles, proteins, nucleic acids, sugars, etc., all of which require the active involvement of transporters and intracellular transportrelated proteins [3,33]. DE-TDFs involved in signal transduction were also highly represented, accounting for a proportion of 11.5% for the total DE-TDFs (Fig. 2). A variety of signaling pathways within Hevea laticifers, including ethylene, jasmonic acid and wound signaling, have been reported to be extensively participate in latex regeneration and regulation [4,17,29,32,[34][35]. The proportions for the two categories, protein degradation and storage and primary metabolism were also high, covering, respectively, 9.3% and 8.2% of the total functional DE-TDFs. Their high representation suggested that with the progress of tapping, in order to meet the balance of supply and demand of all substances in latex regeneration, protein turnover rate becomes faster and primary metabolism gets active. In a word, these results indicated that the latex regeneration regulated by tapping involves a complex multi-gene regulatory network, as well as a physiological and biochemical response process.

Sugar metabolism and rubber biosynthesis in tappingstimulated latex production
In regularly tapped Hevea trees, the main metabolic activity of the laticifers is latex regeneration, which centers on the biosynthesis of rubber hydrocarbon that accounts for about 90% of the dry weight of fresh latex [3]. Sucrose has been identi ed as the precursor material for rubber biosynthesis in laticifers, providing the carbon skeleton and energy required for latex regeneration [11,36]. In Hevea trees normally tapped at intervals of 2-4 days, each tree produces dozens to hundreds of milliliters of fresh latex, and the removed latex could be effectively recovered before the next tapping to ensure the sustained productivity of the tree [3,12]. Therefore, the laticifers are an active carbon sink, and the effective supply of sucrose is a key factor to determine the yield of latex [10,37]. In this study, the genes of a sucrose transporter and a sugar transporter were among the DE-TDFs identi ed, both of which were signi cantly up-regulated with the increase of tappings (Table S1). Interestingly, the former is just the sucrose transporter HbSUT3 we previously identi ed to be critical in sucrose uptake into laticifers and rubber production in exploited Hevea trees [12]. The up-regulation of these two transporters indicated an active involvement of sucrose and sugar transport in tapping-stimulated latex regeneration. Sucrose catabolism and the following pathways of glycolysis, tricarboxylic acid cycle and pentose phosphate provide essential components, i.e. the carbon skeleton (acetyl CoA), the reducing power (NADPH) and the energy (ATP) for the nal rubber biosynthesis pathway [3,11]. Therefore, sugar metabolism becomes one of the core metabolic pathways contributing to latex regeneration and rubber biosynthesis in Hevea [3,11,38,[40][41]. This study identi ed multiple DE-TDFs involved in sucrose cleavage and the three above mentioned sugar metabolism pathways (Table 4; Tables S1 to 3), including those encoding neutral/alkaline invertase, fructokinase, phosphofructokinase, glyceraldehyde 3-phosphate dehydrogenase, pyruvate kinase, pyruvate dehydrogenase, and glucose-6-phosphate dehydrogenase, etc. Most of them were up-regulated in the latex for the rst few tappings (Table S1; Fig. 3). It is worth noting that the up-regulated DE-TDF (M16-A7-1) as identi ed by both cDNA-AFLP (Table S1) and qRT-PCR (Fig. 3) turned out to be HbNIN2, the neutral/alkaline invertase that is responsible for sucrose catabolism in Hevea laticifers [13].
There are 20 gene families directly involved in the rubber hydrocarbon biosynthesis and termed as rubber biosynthesis genes [17,31]. The DE-TDFs identi ed in this study involved six of these families, including cis-prenyltransferase, hydroxymethylglutaryl coenzyme A synthase, 3-hydroxy-3-methylglutaryl-coenzyme A reductase, farnesyl diphosphate synthase, rubber elongation factor and small rubber particle protein (Table 4). Among the nine DE-TDFs, eight were demonstrated by qRT-PCR to be up-regulated with the tappings (Fig. 3). In addition, a DE-TDF (M8-A5-6) annotated as inorganic pyrophosphatase was also bolstered by the tapping treatment (Fig. 3). A vacuolar type of inorganic pyrophosphatase has been found to locate on rubber particles and essential for the incorporation of IPP monomers into elongating rubber molecules [40].

Strength and weakness of the cDNA -AFLP technique
The cDNA-AFLP technique has been widely applied in various eukaryotes including the Hevea tree for transcript pro ling due to its advantages of stringency, reproducibility, cost-effectiveness, genome-wide coverage and the ability to distinguish among highly homologous genes [24,39,[41][42][43]. In this study, as determined by sqRT-PCR and qRT-PCR, about 84% and 90%, respectively, of the selected DE-TDFs were veri ed for their cDNA-AFLP pro les (Tables 3 & 4; Fig. 3), re ecting a high reliability of this technique in screening tapping-responsive DE-TDFs in Hevea latex. According to a previous in-silico estimation [27], about 84% of the genes expressed in Hevea latex could be visualized using the silver-staining cDNA-AFLP technique with the restriction enzyme pair of Apo I and Mse I exploited here. The sucrose transporter HbSUT3 [12] and the neutral/alkaline invertase HbNIN2 [13] that have been reported to be up-regulated in the latex of virgin Hevea trees by the tapping treatment were among the DE-TDFs identi ed in this study (Table 4; Fig. 3; Table S1), demonstrating a high transcript coverage of this technique. However, compared with the currently popularly used next generation RNA-sequencing technique that relies on expensive DNA sequencers and specialized bioinformatics [44], the cDNA-AFLP is labor-intensive. Nevertheless, the cDNA-AFLP technique still have its niche among the various transcript pro ling techniques, and can be readily established in a mediocrely equipped and stringently funded lab to ful ll its customized transcript pro ling task.

Conclusions
A genome-wide cDNA-AFLP transcript pro ling identi ed a total of 505 tapping-responsive DE-TDFs in the rubber-producing laticifers of virgin Hevea trees. According to the 366 DE-TDFs with de nite functional annotation, the tapping treatment brought about extensive physiological and molecular changes in laticifers. The integration of these changes upgraded the mediocre level of laticifer metabolism in virgin trees to a high dynamic equilibrium of latex regeneration in regularly tapped trees. Further combined studies of transcriptomics, proteomics and metabolomics will bene t a deeper insight into the exact relationships (synergy or antagonism) among the vast number of biological pathways implicated in tapping stimulated latex production.

Plant materials
Hevea trees of Reyan7-33-97 clone, planted for eight years in the experimental eld of Chinese Academy of Tropical Agricultural Sciences (Danzhou, Hainan), were rst brought to tapping with a system of half spiral, every 3 days, no ethylene stimulation (S/2, d/3).

Extraction of latex total RNA
Five rubber trees attaining the tapping standard (trunk girth >= 50 cm at 1 m above the ground) were selected for latex collection. Twenty seconds after tapping, about ve mL of latex was allowed to ow into a centrifuge tube containing 5ml 2×RNA extraction buffer (0.3 M LiCl, 10 mM EDTA, 10% SDS, 100 mM Tris-HCl, pH8.0). The collected latex was placed in ice box and brought to laboratory for RNA extraction as described in Tang et al. [45]. Electrophoresis on a 1.5% formaldehyde denaturing agarose gel was used to detect the integrity of RNA samples.

Semi-quantitative reverse transcription PCR (sqRT-PCR)
The rst strand of cDNA was synthesized by reverse transcriptase kit (RevertAid TM First Srtand cDNA Synthesis Kit, Thermo), and then diluted ten times as the template for sqRT-PCR with 18S rRNA used as the reference. The amount of cDNA samples used in sqRT-PCR was adjusted to be the same for the ve tappings based on the level of 18S rRNA expression.

cDNA-AFLP analysis
The manipulations were conducted according to the procedure which we previously established for transcript pro ling in Hevea latex [27]. A total of 50 μg latex total RNA taken from each of the ve tapping samples was subjected to cDNA-AFLP analysis. The synthesized double-stranded cDNA was cut by the restriction enzymes of Apo I and Mse I, and all the 128 possible selective primer combinations with 8 Apo I primers and 16 Mse I primers ( the same as the pre-ampli cation in cDNA-AFLP analysis [27]. PCR products were fractionated by 1.2% agarose gel electrophoresis, and the target band was sliced and puri ed using AxyPrep TM DNA Gel Extraction Kit ( AxyGen).

DE-TDFs cloning and sequencing
The puri ed DE-TDFs were ligated with the T-vectors using the pMD18-T Vector Kit (Takara) in accordance with the manufacture's manual. The ligation mixture was used to transform E. coli JM109 competent cells and the transformants were sent to BGI Genomics Co., Ltd for sequencing.

Bioinformatics analysis
For DE-TDFs analysis, sequences of vectors and adaptors were rst trimmed off by using the VecScreen program on the NCBI website (https://www.ncbi.nlm.nih.gov/tools/VecScreen). Then, the clean DE-TDF sequences were subjected for homology analysis to publicly available GenBank non-redundant sequences databases (http://www.ncbi.nlm.nih.gov) using the BLASTX program. Also, the Gene Ontology (http://amigo1.geneontology.org/cgi-bin/amigo/go.cgi) database was used to investigate the molecular function of each DE-TDF in the cell, which was used as the basis for functional classi cation.
Quantitative RT-PCR (qRT-PCR) The expression pattern of candidate genes was detected by qRT-PCR. The rst strand of cDNA was diluted 20 times as the template for qRT-PCR with HbYLS8 as the reference gene as recommended in our previous study [46].

Availability of data and materials
All data generated or analyzed during this study are included in this published article and its supplementary information les. The datasets measured and analyzed during the current study are available from the corresponding authors upon reasonable request.