Transcriptomic, Proteomic, and Metabolic Profiles of Catalpa bungei Tension Wood Reveal New Insight Into Lignin Biosynthesis Involving Transcription Factor Regulation

Lignin is a complex polymer in plant cell walls whose proportion is second only to that of cellulose and plays an important role in the mechanical properties of wood and stress resistance of plants. Here, we induced tension wood (TW) formation in Catalpa bungei by artificial bending and analyzed the lignin metabolism of the TW. LC-MS analysis showed that a significantly higher content of coniferyl aldehyde was observed in the TW cell wall than in the opposite wood (OW) and normal wood (NW) cell walls. TW had significantly lower contents of coniferyl alcohol than OW and NW. Raman spectroscopy results indicated that TW had lower total lignin than OW and NW. The transcription and translation levels of most of the differentially expressed genes (DEGs) involved in lignin monomer biosynthesis indicated upregulation in TW/OW and TW/NW. We found no significant difference in the transcription levels of three collision gases (CADs) between TW and OW or between NW, but their translation levels were significantly downregulated in TW, suggesting post-transcriptional control for CAD. We predicted and analyzed transcription factors that could target DEGs involved in lignin monomer biosynthesis in TW. Based on the analysis of the relationships of targeting and coexpression, we found that NAC (evm.model.group1.695) could potentially target 4CLs and CCoAOMT, that HD-Zip (evm.model.group7.1157) had potential targeting relationships with CCoAOMT, F5H, and CCR, and that their expression levels were significantly positive. It is speculated that the upregulation of NAC and HD-ZIP transcription factors activates the expression of downstream target genes, which leads to a significant increase in coniferyl aldehyde in TW. However, the decrease in total lignin in TW may be caused by the significant downregulation of CAD translation and the significant decrease in precursors (coniferyl alcohol). Whether the expression of CAD genes is regulated by post-transcriptional control and affects TW lignin metabolism needs further study.


INTRODUCTION
Wood, the secondary xylem of trees, is a multipurpose renewable material and is commonly used in papermaking, architecture, furniture manufacturing, and energy production due to its varied properties (Nakahama et al., 2018;Papadopoulos et al., 2019). Wood formation is considered secondary growth and involves the formation of vascular tissue, secondary cell wall (SCW), lignification, programmed cell death, and heartwood formation (Tian et al., 2007). Wood quality is mainly determined by cellulose, hemicellulose, and lignin contents (Nakahama et al., 2018). Lignification can increase the hardness and mechanical supporting force of the secondary xylem cell wall and enhance the ability of plants to transport water over long distances and resist biological and abiotic stresses (Vance et al., 1980;Zhao, 2016;Hu et al., 2019). Lignin is related to the properties of pulp and paper (Liang et al., 2020). Therefore, it is of great significance to study the distribution, synthesis and regulation of lignin in xylem for the genetic improvement of wood properties.
The lignin synthesis pathway has been clearly elucidated (Zhao, 2016). Generally, lignin in higher plants is formed by the oxidative polymerization of three kinds of lignin monomers (p-hydroxyphenyl lignin, guaiacyl lignin, and syringyl lignin). Studies have found that lignin monomer synthesis pathways in different plants are significantly different (Li et al., 2014). For example, caffeoyl shikimate esterase (CSE) is a key gene for lignin monomer synthesis in Arabidopsis thaliana, and mutation of this gene will greatly reduce the lignin content in A. thaliana (Shen et al., 2013;Vanholme et al., 2013). Overexpression of the F5H gene in A. thaliana with a defective COMT gene leads to enrichment in unusual lignin (5H lignin) (Weng et al., 2010). These results indicate the flexibility of the lignification process. This flexibility is the starting point for the artificial regulation of lignin biosynthesis and will increase the possible applications of many economically important plants. Research has found that reducing the expression of the 4CL gene during lignin synthesis in poplar will reduce Klason lignin by 30%, which reduces the elasticity modulus by 40% (Horvath et al., 2010). However, increasing 4CL gene expression in poplar will increase lignin content by approximately 40%, reduce other carbohydrate contents and significantly increase wood strength (Hu et al., 2019). Altering the expression of genes in the lignin synthesis pathway has been shown to significantly influence the chemical composition of cell walls. Wu et al. (2019) confirmed that lignin monomer biosynthesis was strictly regulated by genes in synergy. They found that downregulation of F5H in COMT-RNAi transgenic switchgrass inhibits S lignin biosynthesis and increases guaiacyl lignin units. In contrast, overexpression of F5H in COMT-RNAi transgenic plants will reduce the number of G units. Deficient lignin biosynthesis could be partially compensated for or completely recovered with different degrees of COMT downregulation. These factors make it possible to artificially and precisely regulate the synthesis of lignin monomers.
Wood formation is a complicated process under the control of a large set of transcription factors (Sun et al., 2021). Lignin biosynthesis and lignification are two of the steps in wood formation. Transcriptional regulation is important for lignin synthesis. Studies have indicated that the lignin metabolic pathway usually has the same transcriptional regulation mode as cellulose and hemicellulose biosynthesis (Zhao and Dixon, 2011;Nakano et al., 2015). NAC transcription factors are usually the first-layer managers in the regulation network of wood chemical composition biosynthesis. NACs can target a number of MYB genes (second-layer managers) and promote their direct or indirect regulation of the formation of secondary walls for plants, including lignin biosynthesis. In A. thaliana, MYB46, MYB58, and MYB63 are targeted by SND1, NST1, and NST2 (Zhou et al., 2009;Zhong et al., 2010). In atnst1-1 and atnst3-1-1 double mutants, secondary wall thickening was completely inhibited, and some lignin synthesis pathway genes were also inhibited (Mitsuda et al., 2005). Moreover, lignin biosynthesis and cellulose deposition genes were highly expressed in the myb75-1 mutant with thickening of the fiber cell wall, increasing the lignin content and changing the S/G ratio in the stem (Bhargava et al., 2010). The PtrWND genes could induce the expression of wood biosynthetic genes, including related transcription factors and structural genes, leading to ectopic deposition of lignin in poplars . Furthermore, several studies have shown that PtrMYB3 and PtrMYB20 can activate the lignin biosynthesis pathway (McCarthy et al., 2010). In addition to the NAC and MYB transcription factors, an increasing number of transcription factors have been found to be involved in wood formation. For example, homeobox genes are associated with lignification in bamboo (Xu et al., 2019), and PtrHB7 expression enhances the differentiation of cambial cells toward xylem cells in Populus (Zhu et al., 2013). HD-ZIP also regulates lignin biosynthesis in plants (Sun et al., 2020). These studies suggest a new level of complexity in lignin biosynthesis regulation.
Catalpa bungei is native to China and is an economically important ring-porous tree species with high-quality wood (Li et al., 2019). Traditional breeding for rapid growth of C. bungei has been carried out for years (Xiao et al., 2019). However, currently, the forest industry demands not only high yields but also high quality. Thus, it is necessary to study the mechanism of wood formation in C. bungei. Recently, with the development of high-throughput sequencing techniques, multiomics analysis has played an important role in the selection of functional genes and the construction of gene expression regulatory networks for wood formation (Shinya et al., 2016;Bygdell et al., 2017;Wang et al., 2018). We combined transcriptome, proteome, and confocal Raman imaging techniques to analyze the xylem of tension wood (TW), opposite wood (OW), and normal wood (NW) in C. bungei and identified the key genes that regulate lignin biosynthesis, providing new insight for wood property improvement by genetic engineering in C. bungei.

Plant Material
One-year-old C. bungei clone 8,402 was planted at the experimental base in Luoyang city, Henan Province. We used artificial bending to induce TW formation in July. The method was as follows: one end of twine was tied to the top of the plant, and the other end was tied to a heavy object. The object was then dragged to bend the plant. The plants were kept in a bent state with a consistent angle of 45 • in the same direction. Refer to the literature for a diagram of the bending treatment mode (Xiao et al., 2020). After 3 months of treatment, wood samples were collected at a height of 60 cm, and xylem was obtained to make tissue sections. The sections were double stained with saffronfast green, and the detailed staining steps are described in the literature (Xiao et al., 2020). The staining results were used to identify the TW and OW.

Ultrastructure Observation
Sample preparation was carried out essentially as previously described by Xiao et al. (2020). The xylem samples were cut using a razor blade and immediately fixed in 3% (w/v) glutaraldehyde. Samples were dehydrated in an ethanol series (30, 50, 70, and 100%) and embedded in Spurr resin at 60 • C. Transverse resin sections of xylem were cut using a diamond knife for ultrastructural observation under a transmission electron microscope (HITACHI HT-7700, Tokyo, Japan).

Confocal Raman Mapping
A LabRam Xplora confocal Raman microscope (Horiba Jobin-Yvon, Paris, France) with an Olympus confocal microscope (Olympus, Tokyo, Japan) and a linearly polarized 532 nm laser (Ventus VIS 532, Laser Quantum, Chester, United Kingdom) was used to measure the sample size. The power was set to 8 mW. Point-by-point scanning microprobe imaging was employed to acquire spectral data. The parameter settings were as follows: 1,200 mm −1 grating, 100 µm slit width, 300 µm numerical aperture, 0.7 µm scanning step, 1 s single-point spectral acquisition time, and 2 cm −1 spectral resolution. Every 15 spectra were acquired and averaged. Data were acquired and analyzed by LabSpec6 software. The entire experiment was carried out at a constant 25 • C.

Liquid Chromatography-Mass Spectrometry (LC-MS) Determination of Metabolites in the Lignin Monomer Biosynthetic Pathway
A total of 100 mg of sample powder was weighed and extracted overnight at 4 • C with 1.0 mL 70% aqueous methanol. Following centrifugation at 10,000 × g for 10 min, the extracts were subjected to SPE (CNWBOND Carbon-GCB solid phase extract (SPE) cartridge, 250 mg, 3 mL; ANPEL, Shanghai, China), and the eluates were filtered (SCAA-104, 0.22 µm pore size; ANPEL) before LC-MS analysis.
Linear ion trap (LIT) and triple quadrupole (QQQ) scans were acquired on a QQQ-LIT mass spectrometer (Q TRAP), API 6500 Q TRAP LC-MS/MS system, equipped with an ESI Turbo Ion-Spray interface, operating in positive ion mode and controlled by Analyst 1.6 software (AB Sciex). The ESI source operation parameters were as follows: ion source, turbo spray; source temperature, 500 • C; ion-spray voltage (IS), 5,500 V; ion source gas I (GSI), 55 psi; gas II (GSII), 60 psi; curtain gas (CUR), 25.0 psi; and collision gas (CAD), high. Instrument tuning and mass calibration were performed with 10 and 100 µmol/l polypropylene glycol solutions in the 6500 Q TRAP and LIT modes, respectively. Q TRAP scans were acquired as multiple reaction monitoring (MRM) experiments with the CAD (nitrogen) set to 5 psi. The declustering potential (DP) and collision energy (CE) for individual MRM transitions were determined with further DP and CE optimization. A specific set of MRM transitions was monitored for each period according to the metabolites eluted within this period.

Transcriptome Sequencing and Transcript Quantification
Three biological replicates of TW (pool of five trees), OW (pool of five trees), and NW (pool of five trees) were used to carry out multiomic sequencing. Total RNA from the TW, OW, and NW of each sample was isolated using an RNA reagent kit (DP441; Tiangen Biotech, Beijing, China) following the manufacturer's protocol. RNA quality was assessed on an Agilent 2100 bioanalyzer (Agilent Technologies, Palo Alto, CA, United States) and checked using RNase-free agarose gel electrophoresis. In total, nine cDNA libraries (three libraries each for TW, OW, and NW) were constructed for mRNA and small RNA high-throughput sequencing using an Illumina HiSeq TM 4000 platform (Illumina, United States) by Gene Denovo Biotechnology Co. (Guangzhou, China). Each base in the reads was assigned a quality score (Q) with a Phred-like algorithm using SOAPnuke software to assess the quality of the RNA-sequencing (RNA-seq) data (Li et al., 2012).
High-quality clean reads were obtained by removing adaptor sequences. RSEM software was used to quantify transcripts (Li and Dewey, 2011). Fragments per kilobase of transcripts per million mapped reads (FPKM) were used to normalize transcript expression. The transcripts with a fold change ≥ 1.5, a p-value < 0.05, and an FPKM > 1 in a comparison were significant differentially expressed genes (DEGs).

Protein Extraction and Labeling
Xylem samples of TW, OW, and NW were ground in liquid nitrogen into powder and transferred to 2 mL of lysis buffer [8 M urea, 2% SDS, 1 × Protease Inhibitor Cocktail (Roche Ltd., Basel, Switzerland)] and sonicated on ice for 30 min, and the homogenate was centrifuged for 30 min at 13,000 rpm and 4 • C.
The supernatant was transferred to a fresh tube. For each sample, proteins were precipitated with ice-cold acetone at −20 • C overnight. The precipitates were washed with acetone three times and redissolved in 8 M urea by sonication on ice. Protein quality was examined by SDS-PAGE (Supplementary Figure 1).
Then, 100 mL of 100 mM triethylammonium bicarbonate (TEAB) and 5 µL of 200 mM Tris(2-carboxyethyl)phosphine (TCEP) were added per 100 µg of protein. The mixture was incubated at 55 • C for 1 h. Then, 5 µL of 375 mM iodoacetamide was added and incubated for another 30 min. Ice-cold acetone was added to precipitate proteins, which were then redissolved in 100 µL TEAB. Proteins were digested with trypsin (Promega, Madison, WI, United States), and the peptide mixture was labeled with tandem mass tags (TMT). The peptide mixture was redissolved in 20 mM ammonium formate solution, pH 10.0, and then an Ultimate 3000 system (Thermo Fisher Scientific, Waltham, MA, United States) was used to fractionate the mixture by high-pH separation via a linear gradient. The column was reequilibrated under initial conditions for 15 min. The column flow rate was maintained at 1 mL/min, and the column temperature was maintained at 30 • C. See the literature (Xiao et al., 2020) for more details.

Protein Identification and Quantification
Tandem mass spectra were extracted, charge state deconvoluted and deisotoped by Mascot Distiller version 2.6. Then, the mass spectrometry data were transformed into mascot generic file (MGF) files with Proteome Discovery 1.2 (Thermo Fisher Scientific, Pittsburgh, PA, United States) and analyzed using the Mascot search engine (Matrix Science, London, United Kingdom; version 2.3.2). The Mascot database was set up for protein identification using the C. bungei reference transcriptome assuming the digestion enzyme trypsin with one missed cleavage allowed. Mascot was searched with a fragment ion mass tolerance of 0.050 Da and a parent ion tolerance of 20.0 ppm. Carbamidomethyl of cysteine and iTRAQ8plex of lysine and the n-terminus were specified in Mascot as fixed modifications. Deamidation of asparagine and glutamine and oxidation of methionine and acetyl groups at the n-terminus were specified in Mascot as variable modifications.
Protein identifications were accepted if they could achieve an false discovery rate (FDR) <1.0% by the scaffold local FDR algorithm. Proteins were classified to satisfy the principles of parsimony if they contained similar peptides and could not be differentiated based on MS/MS analysis alone. In total, 17,548 unique peptides corresponding to 5,366 proteins were quantified. Proteins identified with unique spectra ≥2 in all the samples were used for quantification. The ratios of reporter ions reflect the relative abundances of peptides to ensure relative protein quantification. Medians were used to average and quantify the Mascot search results. Proteins with a fold change >1.2 or <0.83 in a comparison and with an unadjusted significance level P < 0.05 were significant differentially expressed proteins (DEPs).

Prediction of Key Transcription Factors
The DEGs related to lignin monomer synthesis between TW and OW or between TW and NW were selected. TBtools  was used to extract 2,000-bp sequences before the CDSs of key genes as potential promoter sequences based on the C. bungei genome data. The promoter sequences of functional genes were imported into the online database PlantTFDB, 1 and the potential transcription factors of target genes were predicted and regulated based on the Arabidopsis genome. The best match between Arabidopsis and C. bungei was found by basic local alignment search tool (BLAST). Then, we analyzed the differential expression of TFs between TW and OW and between TW and NW. The significantly differentially expressed TFs were candidate TFs.

Statistical Analysis
SPSS 22.0 was used for variance analysis to test the significant differences in metabolite content and gene expression among different wood types. The Duncan method was used for multiple comparisons after variance analysis. GraphPad Prism 6.0 was used to draw scatter plots of gene expression and linear regression analysis. Pearson correlation analysis of the gene transcription and translation levels was carried out by SPSS 22.0, and P < 0.05 was the threshold for significance.

The Anatomy of Tension Wood, Opposite Wood, and Normal Wood
Observation of the transverse section morphology of the C. bungei xylem showed that the differentiation of TW vessels was greatly affected. The number of vessels was significantly reduced, and their size was smaller than those in OW and NW (Figure 1). TW had a thinner cell wall than OW and NW according to the ultrastructure of the cell wall and the statistical results for cell wall thickness in our previous study (Xiao et al., 2020). Slice staining showed that TW could be dyed with fast green (Figures 1B,C), indicating lower lignification of TW. These results also suggested that the formation of the TW cell wall was inhibited and that its chemical composition might be greatly changed.

Raman Spectroscopic Imaging and Microspectroscopy of Tension Wood, Opposite Wood, and Normal Wood
Confocal Raman microscopy was used to obtain more information about lignin in different cell wall regions of TW, OW, and NW. Confocal imaging showed that the total lignin signal intensity from OW and NW fibers was greater than that from TW. This indicated decreased lignin deposition in TW (Figure 2). Figure 3 shows that the Raman spectral intensities of the secondary layer (SL), middle layer (ML), and cell corner (CC) of TW were lower than those of OW and NW. According to previous literature (Table 1), characteristic Raman peaks at 1,210, 1,275, and 1,337 cm −1 reflect H, G, and S lignin, respectively. There was no obvious peak at 1,210 cm −1 from the fiber cell wall in this study. However, there was a weak peak at 1,275 cm −1 and a strong peak at 1,337 cm −1 , which indicated a high relative content of S lignin but no H lignin in C. bungei wood.
The peak at 1,607 cm −1 reflects the presence of aromatic lignin skeletal vibrations. In the SL, ML, and CC, the Raman intensity of this peak from OW was approximately 37.27, 14.42, and 15.64%, respectively, higher than that from TW. The Raman intensity of this peak from NW was also higher than that from TW, but the standard deviation in these intensities was large. This result suggests that the difference in lignin deposition between TW and OW or NW is mainly reflected in the S layer. Similarly, the characteristic peak at 1,275 cm −1 , which reflects G lignin in the ML layer, was not different between TW and OW, while in the S layer, there was a 28.31% higher value for OW than for TW. Although the intensity of the peak at 1,275 cm −1 in the ML was different between TW and NW, there was a slightly large error for NW; thus, we could not accurately assess the differences. The results above further confirmed that the lignin metabolism difference in TW fibers mainly occurred in the SL. The Raman intensities of each functional group FIGURE 2 | Confocal Raman spectral imaging of tension wood (TW), opposite wood (OW), and normal wood (NW). Confocal Raman images of G-lignin monomer at 1,275 wavenumber/cm −1 , S-lignin monomer at 1,337 wavenumber/cm −1 , and total lignin at 1,607 wavenumber/cm −1 . SL, secondary wall layer; ML, middle lamella; CC, cell corner.
in lignin molecules in the CC layer were generally higher than those in the other layers, indicating that the lignin deposition in TW, OW, and NW began in the CC layer, and these peaks for OW were all greater than those for TW. Raman confocal imaging also clearly showed the strongest Raman spectral signals at 1,250-1,290 cm −1 , 1,320-1,338 cm −1 , and 1,590-1,620 cm −1 in CC. This implied that lignin synthesis and deposition in TW were inhibited from the beginning ( Figure 3C).

LC-MS of Metabolites in the Lignin Monomer Biosynthesis Pathway in Tension Wood, Opposite Wood, and Normal Wood
To reveal the change in lignin metabolism during TW formation in C. bungei, we used LC-MS to determine the contents of metabolites in the monolignol synthesis pathways in TW, OW, and NW (Figure 4). The results showed that the content of p-coumaryl alcohol was very low in these three types of wood tissues, but coniferyl alcohol and sinapyl aldehyde in TW, OW, and NW were more abundant than p-coumaryl alcohol. This suggested that H lignin was almost absent from the fiber cell wall of C. bungei, and S lignin was dominant, which was consistent with the Raman spectroscopy results. Cinnamic acid and p-coumaric acid in TW were significantly more abundant than those in OW and NW. Interestingly, we found that both p-coumaryl alcohol and coniferyl alcohol were less abundant in TW than in OW, and the content of coniferyl alcohol in TW was significantly lower (41% of that in OW) than that in OW. Unfortunately, this study did not determine the amount of sinapyl alcohol, but we inferred that its content in TW would be lower than that in OW. The main reason why the total lignin content of TW was significantly lower than that of OW was the reduced biosynthesis of cinnamyl alcohol in TW.   Figure 2A). On the other hand, DEPs in TW vs. OW were also significantly enriched to phenylpropanoid biosynthesis (Q-value = 0.014) (Supplementary Figure 2C). The rich factor represents the proportion of DEGs in a metabolic pathway to the total number of genes, and the larger the value is, the higher the enrichment degree is. The rich factor of phenylpropanoid biosynthesis of DEGs was relatively high in TW vs. NW (0.06) and TW vs. OW (0.07) (Supplementary Figure 2B). And the rich factor of phenylpropanoid biosynthesis of DEPs in TW vs. OW was 0.13, it was higher than that in TW vs. NW (0.06) (Supplementary Figure 2D). This indicated the phenylpropanoid biosynthesis was more important for the difference between TW and OW. All of the above results showed the change of phenylpropanoid biosynthesis would influence the TW information.

Transcription and Translation Levels of Genes Involved in Lignin Monomer
Biosynthesis in Tension Wood, Opposite Wood, and Normal Wood The expression of genes and proteins involved in lignin biosynthesis was detected by transcriptome and proteome analyses. Detailed functional annotations of these genes can be found in Supplementary Table 1. Comparison of the expression levels of these genes or proteins between TW and OW and between TW and NW revealed that most of these genes had high levels of both transcription and translation in TW (Figure 5) Figure 5A). Among them, 4CL and PAL had high levels of translation in TW according to the proteome ( Figure 5B). Notably, most CADs in TW had higher transcription than those in OW and NW, but the translation of CADs in TW was lower than that in OW and NW. Correlation analysis of the gene expression levels showed that the transcription and translation levels of the two 4CLs and CCoAoMT genes were both significantly positively correlated (0.816-0.977). The transcription levels of 4CL (evm.model.group6.1597) and CCoAOMT (evm.model.group4.1455) showed an obvious linear correlation with their translation levels ( Figure 5C). In contrast, there was no obvious linear relationship between the transcription and translation levels of the three CADs, and FIGURE 5 | Transcription and translation levels of genes involved in lignin monomer biosynthesis in tension wood (TW), opposite wood (OW), and normal wood (NW). The numbers in panels (A,B) represent the biological replicates. (A) Transcription levels in TW, OW, and NW and differential expression of genes involved in the lignin metabolism pathway between groups, (B) translation levels in TW, OW, and NW and differential expression of genes involved in the lignin metabolism pathway between groups, (C) correlation analysis of the transcription and translation levels of candidate genes. Pearson's r indicates the Pearson correlation coefficient. ** indicates 1% significance. Gene transcription levels and protein expression levels showed in heatmap were normalized.
the Pearson correlation of their transcription and translation levels was also weak.
Quantitative Real-time polymerase chain reaction (qRT-PCR) was used to validate the results of RNA-seq ( Figure 6A). The transcription level of 4CL (evm.model.group6.1597) in TW was significantly higher than that in OW and NW. This result was consistent with the transcription level of CCoAOMT (evm.model.group4.1455). Correlation analysis showed that the results of the qRT-PCR of CCoAOMT were positively correlated (0.810) with its FPKM value from the transcriptome. On the other hand, the transcription levels of the three CAD genes were not significantly different among TW, OW and NW. However, their qRT-PCR results were weakly correlated with the transcriptome sequencing results (Figure 6B), which may FIGURE 6 | qRT-PCR validation of transcripts involved in lignin monomer biosynthesis in tension wood (TW), opposite wood (OW), and normal wood (NW). (A) qRT-PCR of candidate genes. (B) Correlation analysis of qRT-PCR and transcriptome results for candidate genes. Duncan's method was used for multiple comparisons. Different lowercase letters indicate a 5% significant difference. ns, not significant. Pearson's r indicates the Pearson correlation coefficient. ** indicates 1% significance. The expression of NW was as the control in qRT-PCR experiment.
have been due to the specificity of primers or the difference in RNA samples used in this experiment. In contrast, two 4CLs (evm.model.group13.1024, evm.model.group6.1597) and one CCR (evm.model.group6.1507) showed significantly higher transcription and translation levels in TW than in OW and NW. The translation level of 4CL (evm.model.group6.1597) in TW was 1.50 times that in OW and 1.14 times that in NW. However, the translation level of the three CADs FIGURE 7 | The expression of candidate proteins involved in lignin monomer biosynthesis in tension wood (TW), opposite wood (OW), and normal wood (NW). Duncan's method was used for multiple comparisons. Different lowercase letters indicate a 5% significant difference. One of the biological repetitions of OW serves as a reference for calculating protein relative expression.
in TW was significantly lower than that in OW, and the expression level of CAD (evm.model.group10.214) in TW was significantly lower than those in OW and NW (Figure 7). The translation levels of the three CADs in OW were approximately 16.38-19.20% higher than those in TW (Figure 7 and Supplementary Table 2).
Proteins are better for explaining biological problems as the basis for their function in live organisms. Therefore, the significantly low level of CAD enzymes in TW was probably the main factor limiting monolignol biosynthesis, which eventually led to a decrease in lignin deposition. We further mapped the 17 DEGs whose transcription and translation levels were both significantly different in the lignin monomer metabolic pathways (Supplementary Figure 3). This figure visually shows the transcription and translation levels of genes involved in each chemical reaction in the lignin monomer synthesis pathway during TW formation, increasing the understanding of the mechanism of lignin monomer biosynthesis in C. bungei.

Transcriptional Regulatory Network of Lignin Monomer Biosynthesis in Tension Wood
Transcription factors are some of the most important molecules that regulate gene expression. We used the PlantTFDB online database to predict the transcription factors that potentially regulate 17 candidate genes involved in lignin monomer synthesis. On the basis of the targeted relationship, we selected the differentially expressed transcription factors in TW and OW or in TW and NW as candidate regulatory genes ( Table 2). In this study, 14 transcription factors were identified. Unfortunately, these 14 proteins were not identified in the proteome. Therefore, only their transcriptional levels are shown here. HD-ZIP (evm.model.group7.1157), LBD (evm.model.scaffold464.1), MYB (evm.model.group2.1745), MYB_related (evm.model.group2.388), and NAC (evm.model.group1.695) were significantly differentially expressed in TW and OW, and the remaining transcription factors were significantly differentially expressed in TW and NW. Figure 8A shows the regulatory relationship between the candidate transcription factor and 17 target genes. Each transcription factor has multiple target genes with a node degree between 3 and 11 ( Table 2). The node degree could represent the importance of transcription factors in the metabolic pathway. The node degrees of bHLH (evm.model.group1.1335) and HD-ZIP (evm.model.group7.1157) were 11 and 10, ranking first and second, respectively. Thus, they are considered important transcription factors. They all had targeting relationships with the CCoAOMT, F5H, CCR, and CAD genes (Figures 8B,C). Second, we found that several genes that had a potential targeting relationship with NAC (evm.model.group1.695), including 4CL, CCoAOMT, and CCR, were more directly related to cinnamaldehyde biosynthesis. Therefore, we believe that NAC is also an important regulatory factor affecting lignin monomer synthesis ( Figure 8D). The matching sequence information of the three key transcription factors and target genes is shown in Supplementary Table 3.
The correlation analysis of gene expression (Figure 9) showed that although bHLHs (evm.model.group1.1335) had a high node degree, the correlation between their expression level and target gene transcription and translation levels was weak. Therefore, bHLHs were not considered to have a critical regulatory role in the synthesis of lignin monomers. In contrast, the transcription level of HD-ZIP (evm.model.group7.1157) was significantly positively correlated with the transcription and translation levels of CCoAOMT (evm.model.group4.1455), CCR (evm.model.group6.1507) and F5H (evm.model.group3.991), and the correlation coefficient for the former two groups was >0.8. Similarly, NAC (evm.model.group1.695) was significantly positively correlated with the transcription and translation levels of CCoAOMT (evm.model.group4.1455) and two 4CLs (correlation coefficient > 0.7) (Figure 9). HD-ZIP (evm.model.group7.1157) had a higher transcription level in TW (FPKM was 46.08), which was approximately twice that of OW and NW. It is speculated that this transcription factor is the most important transcription factor regulating lignin metabolism in C. bungei TW.

Regulation Pattern Mediated by Transcription Factors in Lignin Biosynthesis in Tension Wood
We preliminarily proposed the transcriptional regulation mechanism of lignin monomer synthesis of C. bungei TW according to the contents of metabolites in the lignin synthesis pathway, the transcription and translation of key genes and the potential regulatory relationship between candidate transcription factors and genes (Figure 10). When the stems of trees grew under stress, the expression levels of the NAC (evm.model.group1.695) and HD-ZIP (evm.model.group7.1157) transcription factors were significantly increased in the stress area. They could bind 4CL, CCoAOMT or F5H genes and activate their expression at the transcription level, which may promote the production of feruloyl-CoA and sinapoyl-CoA. Then, HD-ZIP further promoted the expression of the CCR gene and enhanced the conversion of coenzyme A to cinnamaldehyde, significantly increasing the content of coniferyl aldehyde in TW. However, we speculated that the significantly lower protein content of the CAD gene in TW may inhibit the conversion of coniferyl aldehyde and sinapaldehydes into coniferyl alcohol and sinapyl alcohol, resulting in a significant decrease in lignin monomer and total lignin in TW. The normal transcription level of CAD and the significantly lower translation level in TW suggested that CAD may be regulated at the post-transcriptional level and affect the translation of CAD. However, the specific regulatory mechanism remains to be further explored.

Differences in Lignin Metabolism Among Tension Wood, Opposite Wood, and Normal Wood of Catalpa bungei
Tension wood is wood with unique anatomical and physical characteristics. Its properties have created a starting point for the study of wood formation mechanisms. Lignin is one of the main chemical components in wood. Changes in its content and structure will affect the quality and utilization of wood. Methods of directionally changing lignin are always popular research topics. In this study, the characteristics of lignin metabolism were evaluated by various methods. Histochemical staining showed that the lignin content of TW was significantly lower than that of OW, which agreed with previous studies (Jin and Kwon, 2009;Brereton et al., 2011;Mizrachi et al., 2015). The presence of a gelatinous layer (G layer) with no or an extremely low lignin content is considered the result of decreased lignification of TW (Pilate et al., 2004). Our previous work showed that the TW of C. bungei does not have a G-layer (Xiao et al., 2020); therefore, there may be a special regulation pattern for lignin metabolism during TW formation in C. bungei.
For a more comprehensive analysis of the changes in lignin metabolism during TW formation, LC-MS and confocal micro-Raman spectroscopy were used to further analyze the metabolic differences among TW, OW, and NW in the lignin monomer biosynthesis of C. bungei. Confocal micro-Raman spectrochemical imaging can reflect the differences in chemical components between samples more accurately without any chemical treatment and can also directly show the spatial distribution of lignin at the subcellular level (Gierlinger and Schwanninger, 2006). Analysis with this technique revealed an obvious difference in lignin metabolism between TW and OW or NW. In TW, not only was the total lignin content lower than that in OW and NW, but the contents of S lignin and G lignin were also lower than those in OW and NW. The Raman images of TW in poplar (Gierlinger and Schwanninger, 2006), maple, beech, and oak (Lehringer et al., 2008) were similar to our images. Notably, our results showed that TW had a significantly higher level of coniferyl aldehyde but a significantly lower level of coniferyl alcohol than OW and NW. Coniferyl alcohol, as one of the precursors of lignin monomers, had a great influence on the total lignin content. We speculate that the enzymes that catalyze this chemical process might have led to this result. In addition, higher contents of p-coumaraldehyde and sinapyl aldehyde and a lower content of p-coumaryl alcohol in TW than in OW and NW were also observed, although the difference did not reach the level of statistical significance. The conversion of cinnamaldehyde to cinnamyl alcohol is the main inhibitor of TW lignin biosynthesis. Unfortunately, sinapyl alcohol was not measured in this study, which is important evidence for our hypothesis.

Gene Expression Involved in Lignin Biosynthesis in the Tension Wood of Catalpa bungei
To test and verify our hypothesis, transcriptome and proteome sequencing of TW, OW, and NW were carried out. The results showed that most of the DEGs and proteins in the phenylpropanoid biosynthesis pathway between TW and OW were highly expressed in TW. The transcriptome of TW in Betula luminifera was similar to our transcriptome (Cai et al., 2018). However, the results were opposite for Populus, Eucalyptus, and Betula platyphylla, in which genes or proteins involved in lignin monomer synthesis were mostly downregulated during TW formation (Andersson-Gunneras et al., 2006;Chen et al., 2015;Mizrachi et al., 2015;Bygdell et al., 2017). This suggests that the expression of genes involved in lignin metabolism is highly dynamic during TW formation. In particular, we found that 4CL (evm.model.group6.1597) and CCoAOMT (evm.model.group4.1455) had high transcription and translation levels in TW, and their transcription and translation levels were significantly positively correlated. Thus, we believe they might be the key to the change in TW lignin biosynthesis in C. bungei. Their high expression could lead to high contents of feruloyl-CoA and sinapoyl-CoA in TW. On the other hand, the high transcription and translation of CCR (evm.model.group6.1507) may be the reason for the high coniferyl aldehyde content in TW. However, in this study, we found that the content of lignin was lower in TW than in OW and NW, contrary to the increase in the expression of most genes involved in lignin monomer biosynthesis in TW. Interestingly, we found that the translation levels of the three CADs were significantly decreased in TW compared with OW, even though their transcription levels were not different between TW and OW. Compared with mRNA, protein, as the final implementer of biological functions, plays a direct role in determining life processes. CAD is the key enzyme leading to the synthesis of coniferyl alcohol and sinapyl alcohol and catalyzes the final step of lignin monomer biosynthesis (Mansell et al., 1974). Therefore, the change in CAD translation levels may be a more important factor in the trait development of TW. Eudes et al. (2006) found that AtCAD1 partially complemented the A. thaliana CAD mutant and resulted in the restoration of lignin units. Baucher et al. (1996) showed that the total lignin content was not changed in transgenic poplar trees with antisense inhibition of CAD gene expression compared with normal poplar, but the cinnamaldehyde level was significantly increased in the transgenic plants. Furthermore, an increase in cinnamaldehyde and an unchanging lignin content were observed in transgenic tobacco with inhibited expression of CAD (Halpin et al., 1994). Natural CAD mutants of pine showed not only high cinnamaldehyde levels but also a 9% decreased lignin content (Mackay et al., 1997). All these studies indicated that CAD plays an important role in lignin biosynthesis. According to our proteomic and LC-MS results, we believe that the low expression of CAD protein in TW resulted in inhibition of lignin monomer biosynthesis and a lower content of lignin in TW. Similarly, Mauriat et al. (2015) found that the CAD enzyme was less abundant in the TW than in OW of poplar. However, more information is needed to determine the difference between the transcription level and translation level of CAD genes in TW.

Molecular Regulation of Lignin Biosynthesis in the Tension Wood of Catalpa bungei
Many transcription factors, such as MYBs and NACs, have been proven to regulate lignin biosynthesis, and this regulatory hierarchy is similar to that of cellulose (Zhou et al., 2009;Zhao and Dixon, 2011). In our study, we also found that NAC (evm.model.group1.695) could regulate the expression of 4CLs and CCoAOMTs to influence lignin monomer biosynthesis. Yang et al. (2019) found that overexpressed WOOD-ASSOCIATED NAC DOMAIN protein 3 (PdWND3A) in Populus deltoides can increase the lignin content and syringyl and guaiacyl lignin (S/G) ratio of transgenic plants, while the expression of F5H was increased. This result suggested that NAC may regulate the expression of the F5H gene to change lignin monomer synthesis in poplar. A recent study showed that the Eucalyptus NAC transcription factor can bind to the promoter of the CCoAOMT gene to activate its expression, and overexpression of this gene in A. thaliana will greatly promote xylem development and lignin synthesis (Sun et al., 2021). With the development of molecular biology techniques, more knowledge about the lignin synthesis mechanism, beyond the NAC-MYB transcriptional cascade, has been revealed (Zhao, 2016). In poplar, HD-ZIP transcription factors are highly expressed in the developing xylem (Schrader et al., 2004). Moreover, HD-ZIP transcription factors have been identified as important determinants of vascular bundle differentiation in A. thaliana (Zhong and Ye, 2004). With further research, Roc8 (one HD-ZIP TF) was found to be able to target the expression of LAC, a lignin synthesis gene, in rice, while knocking out OsRoc8 led to a significant decrease in the lignin content (Sun et al., 2020). These results supported our finding that there was a potential targeting relationship between HD-ZIP (evm.model.group7.1157) and several genes involved in lignin monomer synthesis, and the expression levels were significantly positively correlated with the target gene abundance. This implies that HD-ZIP plays a critical role in the wood formation of C. bungei. Interestingly, the high transcription level and low translation level of CAD in this study also suggest the existence of post-transcriptional regulation. And we guess maybe there were some miRNAs involved the regulation process, cause miRNAs were important to post-transcriptional regulation of genes. On the other hand, CAD is a key enzyme in the last step of the biosynthesis of lignin monomers, the low protein expression level of CAD is likely an important factor affecting the synthesis of lignin monomers when most lignin synthesis genes are highly expressed. Therefore, the hypothesis that miRNAs regulate CAD protein expression in translation inhibition and influence TW lignin metabolism is indeed worth exploring. Previous studies have shown that miRNAs are involved in the regulation of lignin metabolism; for instance, miR397 negatively regulates the LAC gene to change the lignin content of poplar (Lu et al., 2013). miR6433 influences the synthesis of S lignin by binding with the F5H gene (Fan et al., 2020). Overexpression of miR156 can reduce the S/G ratio in poplar (Rubinelli et al., 2013). All of the above results suggest a new level of complexity in lignin biosynthesis regulation.

CONCLUSION
In this study, we analyzed the differences in lignin biosynthesis in C. bungei at the transcription, translation, metabolism, and histochemistry levels between TW and OW or NW. We found that NAC and HD-ZIP transcription factors may be the key transcription factors regulating the biosynthesis of lignin monomers in C. bungei TW by analyzing the matching relationship between the transcription factors and the promoter sequences of key genes and the correlation of their transcription and translation levels. The significantly increased expression of NAC and HD-ZIP in TW may activate the expression of the 4CL, CCoAOMT, F5H, and CCR genes and enhance the production of coniferyl aldehyde in TW. In contrast, the significantly decreased translation of CAD genes in TW may be the reason for the significantly decreased content of coniferyl alcohol, which ultimately limits lignin monomer biosynthesis in TW. However, the regulatory mode that changes the translation level of the CAD gene in TW needs further research.

DATA AVAILABILITY STATEMENT
The RNA-seq datasets are available in Sequence Read Archive (SRA) in NCBI with the accession number PRJNA559964.

AUTHOR CONTRIBUTIONS
YX performed the experiments, processed, analyzed, and interpreted the data, and wrote the manuscript with contributions from the other authors. JL, FY, and TZ assisted with the interpretation of the data and the writing of the manuscript. WM, KZ, and HY prepared the plant material and assisted in the data analysis. NL and TZ carried out the manuscript revision. JW conceived the original research plan and assisted with the design and interpretation of the data. All authors contributed to the article and approved the submitted version.

ACKNOWLEDGMENTS
We would like to thank all colleagues and friends who have contributed to this study.