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.

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 were used for subsequent experiments. Polyadenylated mRNA was enriched using oligo (dT) magnetic beads.

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₂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).

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).

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₁₁H₁₂O₄ 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).

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 (Lertpanyasampatha et al., 2014; Pramoolkit et al., 2014) and Ubiquitin (Li et al., 2011; 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 genes (Schmittgen and Livak, 2008). Three technical replicates per sample were analyzed to ensure reproducibility and reliability.

To evaluate whether the RNA-seq data were sufficient for further analysis, we first assessed their global quality. The RNA-seq experiment generated a total of 125,970,348 (TW), 127,449,704 (OW), and 126,820,416 (NW) raw reads. After trimming, 123,519,666 (TW), 124,605,772 (OW), and 124,379,132 (NW) clean reads remained (Table 2). Among the total clean reads, 115,335,701 (TW), 116,188,884 (OW), and 115,545,396 (NW) were mapped to the rubber tree genome with mapping ratios of 93.37% (TW), 93.26% (OW), and 92.90% (NW) (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₂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.

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.

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 down-regulated; 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 up-regulated 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 G- and H-lignin (Figure 4), potentially leading to higher lignin content in OW than NW. Furthermore, three out of six transcripts encoding peroxidases were significantly up-regulated in OW, while the other three were significantly up-regulated 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.

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.

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.

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 non-coding RNAs and the candidate genes identified through our RNA-seq analysis in rubber tree reaction wood formation.

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 (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.

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 (Sun et al., 2021). 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.