RNA Sequencing Reveals Phenylpropanoid Biosynthesis Genes and Transcription Factors for Hevea brasiliensis Reaction Wood Formation

Given the importance of wood in many industrial applications, much research has focused on wood formation, especially lignin biosynthesis. However, the mechanisms governing the regulation of lignin biosynthesis in the rubber tree (Hevea brasiliensis) remain to be elucidated. Here, we gained insight into the mechanisms of rubber tree lignin biosynthesis using reaction wood (wood with abnormal tissue structure induced by gravity or artificial mechanical treatment) as an experimental model. We performed transcriptome analysis of rubber tree mature xylem from tension wood (TW), opposite wood (OW), and normal wood (NW) using RNA sequencing (RNA-seq). A total of 214, 1,280, and 32 differentially expressed genes (DEGs) were identified in TW vs. NW, OW vs. NW, and TW vs. OW, respectively. GO and KEGG enrichment analysis of DEGs from different comparison groups showed that zeatin biosynthesis, plant hormone signal transduction, phenylpropanoid biosynthesis, and plant–pathogen interaction pathways may play important roles in reaction wood formation. Sixteen transcripts involved in phenylpropanoid biosynthesis and 129 transcripts encoding transcription factors (TFs) were used to construct a TF–gene regulatory network for rubber tree lignin biosynthesis. Among them, MYB, C2H2, and NAC TFs could regulate all the DEGs involved in phenylpropanoid biosynthesis. Overall, this study identified candidate genes and TFs likely involved in phenylpropanoid biosynthesis and provides novel insights into the mechanisms regulating rubber tree lignin biosynthesis.


INTRODUCTION
The rubber tree (Hevea brasiliensis) is a deciduous perennial tropical tree native to the Amazon basin that produces natural rubber as well as rubber wood. Rubber wood can be used to make a variety of products, such as pulpwood for paper production and rubber wood-based panels, furniture, and joinery products (Jahan et al., 2011;Teoh et al., 2011;Severo et al., 2016). The economic importance of the rubber tree and increasing demand have promoted its widespread domestication (Priyadarshan, 2017). Given the increasing economic and application value of rubber wood, improving wood productivity and quality has become the focus of rubber tree breeding (Priyadarshan, 2017).
Wood mainly consists of cellulose, hemicellulose, and lignin. Lignin is a macromolecular biopolymer similar to cellulose, and is widely distributed in higher plants (Wang et al., 2018). During plant growth and development, lignin content increases, causing cell wall thickening and plant tissue lignification. As a crucial secondary metabolite of the phenylpropanoid biosynthesis pathway, lignin had profound effects on plant growth and development (Liu et al., 2018;Wang et al., 2019). It is extensively involved in the material transport of vascular bundles, cell wall structural integrity, stem strength, mechanical support, response to pathogens and other environmental stresses (You et al., 2013;Wang et al., 2019). Thus, the regulation of lignin synthesis affects not only lignin accumulation, but ultimately the growth and development of the entire plant.
The genes for many key enzymes in lignin biosynthesis have been identified, with some also shown to influence environmental responses (Liu et al., 2018;Wang et al., 2018;Wang et al., 2019). For instance, up-regulation of POD, PAL, C4H, and 4CL promotes lignin biosynthesis and increases muskmelon resistance to black spot disease (Han et al., 2018;Yan et al., 2019). Moreover, suppression of Os4CL expression leads to reduced lignin content in rice (Oryza sativa; Gui et al., 2011). Similarly, transgenic tobacco over-expressing the peroxidase gene SWPA4 show high POD activity and lignin accumulation (Kim et al., 2008). Several transcription factors regulate the expression of lignin biosynthetic genes (Liu et al., 2018). In different plant species, TFs from the MYB (Ohtani and Demura, 2019), NAC (Ohtani and Demura, 2019), WRKY (Gallego-Giraldo et al., 2016), MADS  and HSF (Zeng et al., 2016) families control lignin biosynthesis by regulating genes related to cell wall synthesis. These findings provided additional insight into the mechanism of lignin biosynthesis and accumulation, that is TFs could participate in lignin biosynthesis.
In this study, transcriptome sequencing of H. brasiliensis reaction wood uncovered genes involved in regulating lignin synthesis was firstly reported in rubber tree. Identifying transcriptome-level changes of lignin synthesis-related genes contributes to our understanding of rubber wood formation and may help elucidate the regulatory mechanisms of rubber tree lignin biosynthesis. Overall, our findings identify the candidate genes and serve as an important theoretical basis to develop breeding strategies that improve rubber tree quality.

Plant Material and Sample Collection
The rubber tree clones (Reyan 7-33-97) used in this study were obtained by in vitro tissue culture and planted in the Hainan University experimental greenhouse (Danzhou,Hainan,China;109°29′25″ E,19°30′40″ N) at the end of June 2016. To explore the genes involved in reaction wood formation, rubber tree trunks of similar diameter (approximately 2 cm) were selected as experimental materials. Three three-year-old rubber trees were bent at a 30°angle for 300 days to induce reaction wood formation. The bending treatment started at 9 a.m. on August 17th, 2020, and ended at 9 a.m. on June 13th, 2021. The xylem samples were collected immediately after completion of the bending treatment. Mature xylem tissue samples from NW (normal wood), TW (tension wood), and OW (opposite wood) were isolated from the same individual to enable comparisons in the same genetic background. The bark was removed from the sampling area, after which TW (upper side of the branch) and OW (lower side of the branch) were collected from the same branch using a sharp chisel, as described in Li et al. (2013). NW, which represents the control for stem xylem tissue, was isolated from the same side of the tree at breast height, approximately 1 m above the ground. All samples were approximately 2 × 1 cm and 4-5 mm deep. Samples were collected in the morning, immediately frozen in liquid nitrogen, and stored at −80°C for RNA isolation.

RNA Extraction and Qualification
Nine samples (NW1, NW2, NW3; TW1, TW2, TW3; OW1, OW2, OW3) from 300 days rubber tree reaction wood were used for RNA extraction. A modified cetyltrimethyl ammonium bromide (CTAB) method was used to isolate total RNA according to Chang et al. (1993). Genomic DNA was eliminated by DNase treatment. RNA degradation and contamination were assessed on 1% agarose gels. A NanoPhotometer ® spectrophotometer (IMPLEN, CA, United States) and the RNA Nano 6000 Assay Kit of the Bioanalyzer 2,100 system (Agilent Technologies, CA, United States) were used to evaluate RNA quality and concentration for further analysis. Only RNA samples with absorption OD 260/280 ratios from 1.9 to 2.2, OD 260/230 ratios ≥2.0, and RNA integrity number (RIN) values greater than 6.8 For Illumina sequencing, fragmentation buffer was added to produce shorter mRNA strands. Single-stranded cDNA was synthesized from the mRNA using random hexamer primers. Double-stranded cDNA was synthesized by adding buffer, dNTPs, and DNA polymerase I. The double-stranded cDNA was purified using AMPure XP beads and subjected to end repair, addition of the poly-A tail, ligation of the sequencing linker, and fragment size selection. Finally, the nine cDNA libraries were subjected to PCR enrichment and sequenced on the Illumina HiSeq 2,500 platform. Clean reads were obtained by removing low-quality sequence fragments caused by instrument errors, reads with low overall quality, 3' ends with base 10 quality score (Q) < 20 (Q −10log error_ratio ), reads containing N blur, any adapter sequences, and any sequences with <20 nucleotides. The clean reads were aligned to the reference (ref) genome sequence (Liu et al., 2020). The read count of each gene was obtained by mapping the clean reads to the ref genome. The read counts were converted into fragments per kilobase of exon model per million mapped reads (FPKM) values.
DEGs were selected based on the following criteria: |log 2 FC| ≥ 1 and P adjust (P adj ) < 0.05. All DEGs were mapped to individual terms in the Gene Ontology (GO) database (http://www. geneontology.org/), and the number of genes per term was calculated. GO enrichment analysis was performed using GOseq software to identify significantly enriched terms. Analysis of gene regulatory pathways was conducted using the Kyoto Encyclopedia of Genes and Genomes (KEGG) Pathway database (http://www.genome.jp/kegg/pathway.html), MapMan software (v 3.6.0RC1; http://mapman.gabipd.org) was used for the functional pathway analysis (Thimm et al., 2004).

Correlation Networks and Promoter Analysis
Co-expression among genes and TFs was assessed based on the Pearson correlation coefficient calculated in R (version 4.0.1) (Supplementary Table S1). The TF-gene interaction networks were visualized using Cytoscape (version 3.7.2; Shannon et al., 2003). Rubber tree gene promoter sequences (2000 bp upstream of the transcription start site in most cases) were retrieved and analyzed. Promoter analysis and prediction of TF binding sites was performed using PlantPAN 3.0 (http://plantpan.itps.ncku. edu.tw/; Chow et al., 2019).

Lignin Quantification
Quantification of lignin content was performed as in Xiong et al. (2005). The sample (1 g) was crushed and then 0.1 g was added into a centrifuge tube. Subsequently, 10 ml of 1% acetic acid solution was added; the sample was mixed and centrifuged at 4,500 rpm for 5 min. The resulting precipitate was washed with 5 ml of 1% acetic acid, followed by addition of 3-4 ml of ethanol and ether mixture (1:1). The mixture was soaked at room temperature for 3 min, and the supernatant was discarded for 3 times. The precipitate in the centrifuge tube was evaporated in a boiling water bath, and 3 ml of 72% sulfuric acid was added to the precipitate. The precipitate was stirred well with a glass rod, and cellulose was dissolved by incubating the sample at room temperature for 16 h. Then, 10 ml distilled water was added to the test tube, stirred well with a glass rod, and placed in a boiling water bath for 5 min. After cooling, 5 ml of distilled water and 0.5 ml of 10% barium chloride solution were added, mixed well, and centrifuged. The precipitate was washed twice with distilled water, and 10 ml of 10% sulfuric acid and 10 ml of 0.1 mol/L potassium dichromate solution were added to the washed lignin precipitate. The test tube was placed in a boiling water bath for 15 min, stirring constantly during the process, and allowed to cool for later use.
All materials from the cooled test tube were transferred to a beaker for titration and the remaining material was washed with 15-20 ml distilled water. Then, 5 ml of 20% KI solution and 1 ml of 0.5% starch solution were added to the beaker and titrated with 0.2 mol/L sodium thiosulfate. Note: Separate titration and add 10 ml of 10% sulfuric acid and 10 ml of 0.1 mol/L potassium dichromate solution as a blank sample.
Lignin content was calculated using the following equation: % lignin k(a−b) n×48 , where 48 is the 1 mol of C 11 H 12 O 4 equivalent to sodium thiosulfate, k is the concentration of sodium thiosulfate (mol/L), a is the volume of sodium thiosulfate consumed in blank titration (ml); b is the volume of sodium thiosulfate consumed by the titration solution (ml); n is the mass of the sample (g).

Validation of Gene Expression by
Reverse-Transcription Quantitative PCR (RT-qPCR) cDNA was synthesized by reverse transcription using 1 μg total RNA from the nine samples (NW1, NW2, NW3; TW1, TW2, TW3; OW1, OW2 and OW3), respectively. Primer Premier v5 software was used to design primers specific to the selected genes ( Table 1). Five genes were randomly selected from the NW, TW, and OW of rubber tree reaction wood. For RT-qPCR analysis, TB Green Premix Ex Taq II (Tli RNaseH Plus; Takara, Beijing, China) was used following the manufacturer's recommendations. PCR amplification was performed in a preheated (94°C) thermal cycler and incubated at 94°C for 2 min, followed by 40 cycles of 95°C for 5 s and 60°C for 30 s. For melting curve analysis, samples were denatured at 95°C for 15 s, then cooled to 60°C (4°C per second). Fluorescence signals were collected at 520 nm continuously from 60°C to 95°C (0.2°C per second). Expression stability of the 40S Pramoolkit et al., 2014) and Ubiquitin Chao et al., 2016) rubber tree genes has been previously evaluated and both were confirmed as suitable internal control genes. Thus, these two genes served as internal controls for normalization. Expression levels of the DEGs were calculated using the 2 −△△Ct method against the internal control Frontiers in Genetics | www.frontiersin.org October 2021 | Volume 12 | Article 763841  TGCAGAACGAAGAGGGTCAG  TTGTTGGCGGAGTTCAACCT  gene-GH714_003817  GCTAATGTCACGGGCTCCAT  TGGGATCGACCCAGATAACG  gene-GH714_035841  ATTCTGGGTCTGGTCTCCCA  CCCCTCCCCAATTCGGTTTT  gene-GH714_043603  GCGAAATCGCTCCGATCAAC  CGAGTCACTCTGGGTTCCAC  gene-GH714_019256  GCTCATGTGGCTTTCATGGC  GAACCCTTCCCCTTGCCTAC  gene-GH714_027628 GGTGAGCACATGCCTTGTTG GCCCTGAGTGAACGAAGTGA  genes (Schmittgen and Livak, 2008). Three technical replicates per sample were analyzed to ensure reproducibility and reliability.

Global Transcriptome Analysis of RNA-Seq Data
To evaluate whether the RNA-seq data were sufficient for further analysis, we first assessed their global quality.  Table 2). Based on previous studies (Geraldes et al., 2011;Chen et al., 2015a), these results indicated that our RNA-seq results detected most expressed genes were sufficient for subsequent quantitative analysis. To measure changes in gene expression and find key genes involved in reaction wood formation, we further selected DEGs that met the following statistical significance criteria: |log 2 FC| ≥ 1 and P adj ≤ 0.05. In total, 214 (TW vs. NW; 173 up-regulated and 41 down-regulated), 1,280 (OW vs. NW; 527 up-regulated and 753 down-regulated), and 32 (TW vs. OW; 26 up-regulated and 6 down-regulated) DEGs were identified in rubber tree reaction wood ( Figures 1A-C). After removing repetitive genes, this analysis identified 1,347 genes that were significantly differentially expressed in the TW, OW and NW xylem tissues ( Figure 1D). The DEGs obtained from three comparison groups were used for subsequent analysis.

GO and KEGG Enrichment Analyses of DEGs That Participate in Reaction Wood Formation
RNA-seq provided an overview of genes that are differentially expressed during reaction wood formation. To better understand the function of these DEGs and their associated biological processes, GO and KEGG enrichment analyses were carried out for the different comparison groups shown in Figure 1D. First, for DEGs in TW vs. NW, GO analysis revealed enrichment for intramolecular oxidoreductase activity, cell wall macromolecule metabolic process, and sulfate transmembrane transporter activity (Figure 2A), while KEGG enrichment analysis revealed that zeatin biosynthesis and plant hormone signal transduction pathways were significantly enriched ( Figure 2B). Next, GO analysis of the DEGs from OW vs. NW showed that biological processes including cellulose synthase (UDP-forming) activity, cellulose synthase activity, and UDP-glucosyltransferase activity were significantly enriched ( Figure 2C), and KEGG analysis showed that phenylpropanoid biosynthesis and plant hormone signal transduction pathway were significantly enriched ( Figure 2D). Last, GO analysis of the DEGs from TW vs. OW showed that cell growth, cell wall assembly, and cellulose microfibril organization were significantly enriched ( Figure 2E). KEGG analysis of DEGs from TW vs. OW showed enrichment for plant-pathogen interaction pathway ( Figure 2F). Taken together, results of the GO and KEGG analyses of DEGs from different treatment groups suggest that zeatin biosynthesis, plant hormone signal transduction, phenylpropanoid biosynthesis, and plant-pathogen interaction pathway might function in the formation of reaction wood.

Expression Analysis of DEGs Involved in Phenylpropanoid Biosynthesis and Lignin Quantification
The phenylpropanoid biosynthesis pathway is responsible for lignin production and thus plays an important role in wood formation. KEGG enrichment analysis indicated that DEGs from OW vs. NW are involved in the phenylpropanoid biosynthesis pathway. A total of 16 transcripts of five genes from OW vs. NW are related to lignin synthesis ( Figure 3A). The 4CL gene was down-regulated in OW; 4CL converts coumarate, caffeate, ferulate, 5-hydroxyferulate, and sinapate to coumaroyl-CoA, caffeoyl-CoA, feruloyl-CoA, 5-hydroxyferuloyl-CoA, sinapoyl-CoA (Figure 4). The CAD gene was also downregulated; CAD converts coumaryl aldehyde, caffeyl aldehyde, coniferaldehyde, 5-hydroxy-coniferaldehyde, and sinapaldehyde to coumaryl alcohol, caffeyl alcohol, coniferyl alcohol, 5-hydroxyconiferyl alcohol, and sinapyl alcohol ( Figure 4). Down-regulation of 4CL and CAD in OW vs. NW may imply that these two conversion steps are suppressed in reaction wood formation. The COMT and CCoAOMT genes were significantly upregulated in OW vs. NW. These two genes are involved in eight steps of lignin synthesis, including in the conversion of caffeyl alcohol to coniferyl alcohol and 5-hydroxyconiferyl alcohol to sinapyl alcohol (Figure 4). Up-regulation of COMT and CCoAOMT would promote the synthesis of Gand H-lignin (Figure 4), potentially leading to higher lignin content in OW than NW. Furthermore, three out of six transcripts encoding peroxidases were significantly upregulated in OW, while the other three were significantly upregulated in NW (Figure 4). These results suggest that changes in POD gene expression may not affect lignin content.
To link the gene expression data with effects on lignin synthesis, lignin levels were quantified in NW, TW, and OW. The average lignin content of each tissue was 263.54 mg/g in OW, 248.68 mg/g in NW, and 240.06 mg/g in TW. There were no statistically significant differences detected in the lignin content of these different tissues (Table 3). These results are consistent with the changes in expression levels of DEGs in the phenylpropanoid biosynthesis pathway.

Transcription Factors Mediated Transcriptional Regulatory Networks Involved in Phenylpropanoid Biosynthesis
DEGs encoding TFs in each comparison group were identified to discover potential transcriptional regulatory networks in phenylpropanoid biosynthesis ( Figure 3B). In total, 129 TFs and 16 DEGs associated with phenylpropanoid biosynthesis from OW vs. NW were used to construct a co-expression network using Pearson correlation coefficients.
Significantly co-expressed TF-gene pairs (|cor| ≥ 0.9 and p < 0.05) were selected to construct the transcriptional regulatory network ( Figure 5). In the network, 27 TF families regulate phenylpropanoid biosynthesis-related DEGs. COMT, 4CL, POD, CCoAOMT, and CAD were regulated by 20, 20, 24, 10, and 12 TF families, respectively ( Figure 5). Notably, the MYB, C2H2, C3H, and NAC families could regulate all the DEGs involved in phenylpropanoid biosynthesis ( Figure 5). We evaluated the promoter sequence of genes co-expressed with TFs for potential TF binding sites. Promoter analysis showed that MYB, C2H2, and NAC TF families could bind to the promoter sequence of all DEGs involved in phenylpropanoid biosynthesis (Supplementary Table S2). These results suggest that these three TF families are likely key players in rubber tree reaction wood formation by regulating numerous downstream genes involved in phenylpropanoid biosynthesis.

Validation of RNA-Seq Gene Expression Data by RT-qPCR
The expression patterns of most genes in the NW, TW and OW showed similar trends between the high-throughput RNA-seq  data and RT-qPCR data ( Figure 6). Although the relative expression level calculated by sequencing did not exactly match the expression values detected by RT-qPCR, the expression profiles were mostly consistent for the genes tested. These results confirm the reliability of the gene expression values generated from the RNA-seq data.

A Model for Reaction Wood Formation Identifies Genes Involved in Wood Formation of Rubber Tree and Other Tree Species
To explore the mechanisms of rubber tree reaction wood formation, wood formation-related genes were identified and their expression patterns were investigated. In this study, we identified 214, 1,280, and 32 DEGs from three comparison groups (TW vs. NW, OW vs. NW, and TW vs. OW), respectively ( Figure 1). Among them, 62, 1,103, and 3 DEGs were unique to TW vs. NW, OW vs. NW, and TW vs. OW, respectively. There were no DEGs common to the three comparison groups ( Figure 1D), which might indicate that these tissue-specific DEGs function in wood formation of different xylem tissues. In addition, GO and KEGG enrichment analyses showed that the DEGs from different comparison groups have unique functions and participate in different pathways (Figure 2), further confirming that these DEGs perform distinct functions in different xylem tissues. Based on the GO and KEGG enrichment analyses, our subsequent work focused on genes involved in phenylpropanoid biosynthesis, as well as TFs identified in the list of DEGs from OW vs. NW ( Figure 2D, Figure 3). The heatmap and interaction network showed that phenylpropanoid biosynthesis-related genes and TFs had divergent expression patterns in OW vs. NW, which may reflect the diverse functions of these genes and the complexity of the regulatory network associated with wood formation (Figures 3, 5). These results further demonstrate that wood development is a complex biological process. Genes involved in wood formation may be effectively identified through our model of reaction wood formation.
Previous research has shown that many factors are involved in reaction wood formation, such as non-coding RNA, single nucleotide polymorphisms etc. Xiao et al. (2020) previously reported that long non-coding RNAs (lncRNAs) might participate in Catalpa bungei tension wood formation by regulating genes involved in indoleacetic acid (IAA) and abscisic acid (ABA) synthesis. Similarly, microRNA (miRNA)-and lncRNA-mediated regulatory networks widely participate in reaction wood formation of Populus tomentosa (Pto), and single nucleotide polymorphisms of miRNAs and their target genes also influence this process (Chen et al., 2015b;Chen et al., 2016). Moreover, Liu et al. (2021) induced tension wood in stems of black cottonwood (Populus trichocarpa) and established a transcriptional regulatory network in which PtrHSFB3-1 and PtrMYB092 directly activate eight and 11 monolignol genes that participate in reaction wood formation, respectively. Future research should focus on characterizing the roles of noncoding RNAs and the candidate genes identified through our RNA-seq analysis in rubber tree reaction wood formation.

Characterization of Genes Associated With Phenylpropanoid Biosynthesis in Rubber Tree by Transcriptome Analysis
Lignin synthesis is positively associated with the expression levels of 4CL, COMT, CCoAOMT, CAD, and POD (Do et al., 2007;Vanholme et al., 2008;Wagner et al., 2011;Hu et al., 2017;Chanoca et al., 2019). In this study, COMT and CCoAOMT were up-regulated in OW vs. NW, suggesting that the lignin biosynthetic process may be enhanced in this tissue. However, 4CL and CAD were down-regulated, indicating that the reaction steps involving 4CL and CAD may be inhibited in OW versus NW. Moreover, three peroxidase gene transcripts were up-regulated in OW, while three were up-regulated in NW. This indicates that peroxidase genes likely generate different transcripts in different rubber tree xylem tissues to maintain lignin biosynthesis.
Previous studies showed that repressing steps of the lignin biosynthetic pathway, from PAL to CAD, may lead to reduced lignin content (Chanoca et al., 2019). In addition, studies have shown that single or double Arabidopsis knockout mutants in several peroxidase genes typically exhibit minor but noticeable reductions in lignin content and/or altered lignin composition in inflorescence stems, suggesting that changes in peroxidase expression levels may not strongly affect lignin synthesis Frontiers in Genetics | www.frontiersin.org October 2021 | Volume 12 | Article 763841 (Herrero et al., 2013;Barros et al., 2015). Thus, changes in peroxidase gene levels observed in this work likely did not affect lignin content between OW and NW. Overall, the lack of a significant difference in lignin content among the different tissue types in this study may be explained by the opposite expression patterns of 4CL, COMT, CCoAOMT, CAD, and by the minor contribution of the peroxidase genes to lignin levels.

NAC and MYB Transcription Factor Families Participate in Rubber Tree Lignin Synthesis
Wood formation is an essential yet complex biological process arising from plant secondary growth. Important transcription factors involved in secondary growth have been identified, such as members of the MYB and NAC TF families (Demura and Fukuda, 2007;Zhong and Ye, 2009). In this study, 129 TFs from 27 TF families were responsive to reaction wood treatment (Figures 3, 5). Among them, TFs from the MYB and NAC families could regulate all DEGs involved in phenylpropanoid biosynthesis and might play pivotal roles in rubber tree reaction wood formation. Previous studies have shown that overexpression of MYB transcription factors PtoMYB92, PtoMYB216, and PtoMYB74 could up-regulate genes involved in lignin biosynthesis and promote the formation of additional xylem layers, thicker xylem cell walls, as well as ectopic lignin deposition; these overexpression plants accumulated 13-50% more lignin (Tian et al., 2013;Li et al., 2015;Li et al., 2018). Similarly, overexpression of NAC141 in Arabidopsis resulted in enhanced expression of lignin biosynthesis genes, stronger lignification, larger xylem, and higher lignin content compared with wild-type plants . These findings support the TF-gene regulatory network for rubber tree wood formation generated in this study, in which MYB and NAC TFs likely play important roles.

CONCLUSION
Identification of DEGs through RNA-seq analysis on mature xylem from TW, OW, and NW of rubber tree reaction wood provided insight into the molecular basis of lignin biosynthesis. Our work demonstrated that zeatin biosynthesis, plant hormone signal transduction, phenylpropanoid biosynthesis, and plant-pathogen interaction pathways likely influence reaction wood development. A transcriptional regulatory network analysis revealed that three TF families (MYB, C2H2, and NAC) could regulate all DEGs involved in phenylpropanoid biosynthesis and play important roles in lignin biosynthesis for rubber tree reaction wood formation. Taken together, the findings presented here advance our understanding of the regulation of phenylpropanoid biosynthesis, a critical pathway for lignin production in perennial trees. The publicly available transcriptomic dataset provides essential information for further transcriptomic, genomic, and functional genomics research in rubber tree. Overall, the results generated in this study will serve as an important resource for future studies on this economically important tropical tree crop.