Wood Transcriptome Profiling Identifies Critical Pathway Genes of Secondary Wall Biosynthesis and Novel Regulators for Vascular Cambium Development in Populus

Wood, the most abundant biomass on Earth, is composed of secondary xylem differentiated from vascular cambium. However, the underlying molecular mechanisms of wood formation remain largely unclear. To gain insight into wood formation, we performed a series of wood-forming tissue-specific transcriptome analyses from a hybrid poplar (Populus alba × P. glandulosa, clone BH) using RNA-seq. Together with shoot apex and leaf tissue, cambium and xylem tissues were isolated from vertical stem segments representing a gradient of secondary growth developmental stages (i.e., immature, intermediate, and mature stem). In a comparative transcriptome analysis of the ‘developing xylem’ and ‘leaf’ tissue, we could identify critical players catalyzing each biosynthetic step of secondary wall components (e.g., cellulose, xylan, and lignin). Several candidate genes involved in the initiation of vascular cambium formation were found via a co-expression network analysis using abundantly expressed genes in the ‘intermediate stem-derived cambium’ tissue. We found that transgenic Arabidopsis plants overexpressing the PtrHAM4-1, a GRAS family transcription factor, resulted in a significant increase of vascular cambium development. This phenotype was successfully reproduced in the transgenic poplars overexpressing the PtrHAM4-1. Taken together, our results may serve as a springboard for further research to unravel the molecular mechanism of wood formation, one of the most important biological processes on this planet.


Introduction
Woody biomass represents more than 90% of the total biomass produced within the earth's terrestrial ecosystems. Roughly 25% of the annual anthropogenic CO 2 emissions can be assimilated during woody biomass formation, suggesting that wood formation serves as one of earth's major long-term terrestrial carbon sinks [1,2]. Woody biomass has the potential to be renewable as well as carbon neutral with regard to its conversion into various forms of energy (e.g., electricity, gas, and liquid energy); thus, it has attracted attention in fields of sustainable energy [3][4][5][6][7][8][9][10].
Wood is primarily produced by woody perennials through a process called secondary growth (i.e., wood formation). Secondary growth is achieved by the vascular cambium, a cylindrical domain of pluripotent stem cells below the organ surface, forming wood (i.e., secondary xylem) inside and bast (i.e., secondary phloem) outside in a strictly bidirectional manner by coordinated cell division and differentiation [11]. However, the underlying molecular mechanisms of vascular Raw reads generated from RNA sequencing were cleaned by PRINSEQ-Lite 0.20.4 (http://prinseq. sourceforge.net/) using a Phred quality score of 20 or above with minimum lengths of 50 bp or greater. For genome-guided transcriptome assembly, a script in Trinity v2.2.0 was used [47,48]. In brief, the cleaned paired-end reads of each library were mapped to the P. trichocarpa genome v3.0 (http://www.phytozome.org) using Bowtie v. 1.2.2 [49]. The transcript abundance (e.g., read count) was determined by RSEM (RNA-Seq by Expectation Maximization, v. 1.3.0) [50] and represented by a FPKM (Fragments Per Kilobase of transcript per Million mapped reads) value. The edgeR was used for statistical analysis of the differential transcript abundance by input read counts [51] using the following script; $run_DE_analysis.pl -counts.matrix -method edgeR -output edgeR_dir -dispersion 0.1 -samples samples.txt. For transcript annotation, BLASTX (e-value 1e-5) was performed against A. thaliana (Athaliana_167_TAIR10) and P. trichocarpa (Ptrichocarpa_210_v3) protein datasets from Phytozome V12 (https://phytozome.jgi.doe.gov/pz/portal.html), respectively.

Quantitative Real-Time PCR (qRT-PCR) and Semi-Quantitative RT-PCR
One microgram of poplar total RNA was reverse transcribed to produce first-strand cDNA using a PrimeScript™ RT reagent kit (Takara, Otsu, Japan) following the manufacturer's instructions. The gene expression patterns were analyzed by quantitative real-time PCR (qRT-PCR) [52]. All primer sequences were designed using the Primer3 program (http://fokker.wi.mit.edu). Poplar Actin2 gene was used as the quantitative control [53]. The qRT-PCR was performed using a CFX96™ Real-time PCR Detection System (Bio-Rad, Hercules, CA, USA) with iQ™ SYBR ® Supermix (Bio-Rad). In Arabidopsis samples, the total RNA was extracted using TRIZOL reagent (Ambion, Austin, TX, USA), according to the suggested protocol, with slight modifications. One microgram of the total RNA was reverse transcribed using Superscript II reverse transcriptase (Invitrogen, Carlsbad, CA, USA) in a 20 µL reaction. RT-PCR was carried out using 1 µL of the two-fold diluted reaction product as a template. Amplified DNA fragments were separated on 1% agarose gel and visualized with ethidium bromide staining.

Vector Construction and Production of Transgenic Plants
The full-length cDNA encoding PtrHAM4-1 (Potri.005G125800) was amplified from hybrid poplar (clone BH) by PCR and inserted downstream of the 35S promoter in the pK2GW7 vector [54] using the Gateway cloning system in order to produce 35S::PtrHAM4-1 constructs. Vector constructs were introduced into Agrobacterium tumefaciens strain C58, which was used to transform Arabidopsis thaliana (Columbia; Col-0) and a hybrid poplar by the floral-dip method [55] and leaf disk transformation-regeneration method [56]. All of the constructs used in this study were verified by DNA sequencing.

Histological Analysis
Rosette level stems of 60-day-old Arabidopsis plants were used to obtain hand-cut cross sections and stained with 0.05% toluidine blue O for 1 min as described previously [56]. Interfascicular cambium-derived tissue (ICD) extensions were measured at all of the interfascicular regions from 60-day-old stem sections immediately above rosette level. In poplar, lengths of cambial layers were measured in a total of eight directions of each stem sections and at least five ramets of each transgenic line (3-month-grown in test tube after subculture) were used for statistical analysis. Images were captured using a microscope (CHB-213; Olympus, Tokyo, Japan) and analyzed by ImageJ software (National Institutes of Health; http://www.nih.gov/).

Generation of a Wood-Forming Tissue-Specific Transcriptome from a Hybrid Poplar
To explore the molecular mechanism of wood formation, we designed a series of tissue type-specific transcriptome analyses using a hybrid poplar (Populus alba × P. glandulosa, clone BH). Each wood-forming tissue was isolated from vertical stem segments representing a gradient of developmental stages with regard to secondary growth (i.e., immature, intermediate, and mature stems; Figure 1a), which was previously confirmed by microscopic observation of the stem cross sections [45]. Thus, a total of seven tissue types from at least 10 actively growing poplars (field grown, one-year-old) was collected and combined, including SL (Shoot apical meristem with Leaf primordia), IS (Immature Stem), IC (Intermediate stem-derived Cambium), IDX (Intermediate stem-derived Developing Xylem), MC (Mature stem-derived Cambium), and MDX (Mature stem-derived Developing Xylem). In addition, L (Leaf without major veins) was included as a negative control for wood formation (Figure 1a). SL includes shoot apex with leaf primordia; IS indicates immature stem tissue of the third to the fifth internodes containing procambium and primary xylem; IC and IDX were collected from intermediate stems (seventh to 10th internodes), while MC and MDX were collected from mature stems (20-25th internodes) with a scraping method after separating bark tissues [45] (see Materials and Methods).
The total RNA extracted from the aforementioned tissues was used directly for RNA-seq analysis (Hiseq 2000, Illumina, San Diego, CA, USA). Seven RNA-seq libraries were constructed and sequenced to produce 6.5 to 7.5 million reads per library, corresponding to more than 6 Gb per sample (Supplemental Figure S1). Using genome-guided assembly, about 78.5% of the sequences were mapped to the Populus trichocarpa genome (v3.0) (Supplemental Figure S2), with a total of 38,329 transcripts across all tissues after removal of redundant transcripts (from a total of 84,168 transcripts) (Supplemental Table S1). The total RNA extracted from the aforementioned tissues was used directly for RNA-seq analysis (Hiseq 2000, Illumina, San Diego, CA, USA). Seven RNA-seq libraries were constructed and sequenced to produce 6.5 to 7.5 million reads per library, corresponding to more than 6 Gb per sample (Supplemental Figure S1). Using genome-guided assembly, about 78.5% of the sequences were mapped to the Populus trichocarpa genome (v3.0) (Supplemental Figure S2), with a total of 38,329 transcripts across all tissues after removal of redundant transcripts (from a total of 84,168 transcripts) (Supplemental Table S1).

Reliability of Wood-Forming Tissue-Specific Transcriptome
Transcriptional interrelationships of the seven tissue libraries were evaluated by generating a sample correlation heatmap and matrix using the Trinity package [48] (Figure 1b,c). Interestingly, the seven tissues segregated into two groups: One with L, SL and IS; and the other with ID, MD, IC and MC (Figure 1b). In addition, MC and IC as well as SL and IS were closely located within the same clades. As expected, this result suggests that physiologically similar tissues have comparable transcriptomes. For example, in the correlation matrix corresponding to each tissue, MC is highly correlated with IC (Figure 1c), which is convincing, because both MC and IC consist primarily of cambial cells, and next comes MDX, which flanks MC. Interestingly, the next most correlative tissue to IC is IS-instead of IDX, which flanks IC-implying that IC is related to IS in containing immature meristem tissues. These results suggest that the tissue-specific transcriptome data is reliable.

Reliability of Wood-Forming Tissue-Specific Transcriptome
Transcriptional interrelationships of the seven tissue libraries were evaluated by generating a sample correlation heatmap and matrix using the Trinity package [48] (Figure 1b,c). Interestingly, the seven tissues segregated into two groups: One with L, SL and IS; and the other with ID, MD, IC and MC (Figure 1b). In addition, MC and IC as well as SL and IS were closely located within the same clades. As expected, this result suggests that physiologically similar tissues have comparable transcriptomes. For example, in the correlation matrix corresponding to each tissue, MC is highly correlated with IC (Figure 1c), which is convincing, because both MC and IC consist primarily of cambial cells, and next comes MDX, which flanks MC. Interestingly, the next most correlative tissue to IC is IS-instead of IDX, which flanks IC-implying that IC is related to IS in containing immature meristem tissues. These results suggest that the tissue-specific transcriptome data is reliable.

Identification of Genes Involved in Vascular Cambium Initiation and Maintenance
Next, we analyzed our transcriptome data to identify genes that were involved in vascular cambium formation, of which little information is available to date, especially for woody perennials. Firstly, we checked expression of homologs to Arabidopsis genes known to be involved in cambium formation (Table 1). Interestingly, most of these were highly expressed in either IC or MC, or both IC and MC. For example, poplar genes homologous to PIN4, LOG1, and AtHAM4 were highly expressed in IC; while PIN1, AHK4/CRE1/WOL1, AHP1, CLE41, CLE44, AtAUR1, MOL1, and PXY homologs were highly expressed in MC. In contrast, WOX4, TMO5-Like1, LHW, TMO3/CRF2, AtHB8, SHR, TMO6, and DOF5.6 were abundantly expressed in both IC and MC when compared to other tissues. These results suggest that poplar has a gene regulatory network of cambial development similar to that identified in Arabidopsis.
To identify novel genes involved in vascular cambium initiation and maintenance within poplar, we attempted to isolate a group of genes preferentially expressed in either IC or MC. To do so, we selected genes preferentially upregulated (>3-fold) in either IC or MC compared to other tissues, which resulted in identification of 885 and 798 transcripts, respectively (Supplemental Tables S3 and  S4). Among these, we initially focused on transcriptional regulators (64 and 58 transcripts) (Tables 2  and 3). Table 2 shows abundantly expressed transcriptional regulators in IC, including homologs to Arabidopsis genes involved in cytokinin signaling (e.g., ARR11 and ARR12), negative regulation of xylem vessel formation (e.g., ANAC083/VNI2 and ANAC104/XND1), and meristem differentiation (e.g., PNF, WOX1, AtHAM3, and AtHAM4). Table 3 shows abundantly expressed transcriptional regulators in MC, including genes related to meristem maintenance (e.g., KNAT6), asymmetric stem cell division (e.g., SCHIZORIZA), and cell proliferation (e.g., AINTEGUMENTA, AINTEGUMENTA-like 5, GRF5). Additionally, abaxial and adaxial cell patterning-related genes (e.g., KANADI1, Dof5.1) are listed. Thus, the transcripts listed in Tables 2 and 3 include many genes known to be involved in vascular cambium initiation and maintenance. These results suggest that this listing is reliable and can be used to identify novel regulators of cambium initiation and maintenance in poplar. For example, seven poplar genes homologous to either Arabidopsis OBP3 or OBP4 are listed in Table 2, but their functional significance in cambial development remains unknown. Previously, it was found that flowering and development of cambium are inversely correlated [57,58]. Interestingly, transcription factors involved in the negative regulation of flowering were also found (e.g., EFM, NTL8, MYR2) ( Table 2).
To discover potential novel regulators involved in the initiation of vascular cambium, we performed a co-expression network analysis using the AspWood website [37] by querying the abundantly expressed transcriptional regulators of IC in Table 2. We found that most of the genes were interconnected with high correlation (Figure 6a). Among them, several genes showed a high degree of prominent centrality based on number of neighbors [37] including Potri.005G125800.1. Both Potri.005G125800.1 and its closest homolog, Potri.007G029200.1, are preferentially expressed in IC (Figure 6b) and show a higher similarity to AtHAM4 compared to other AtHAMs, which were found from both the amino acid sequence alignment ( Figure S3a) and phylogenic analysis ( Figure S3b). AtHAM4 was recently reported as a candidate for regulating cambium initiation in Arabidopsis by interacting with WOX4 [19]. Thus, we designated Potri.005G125800.1 and Potri.007G029200.1 as PtrHAM4-1 and PtrHAM4-2, respectively. The high spatial-resolution wood formation data [37] further support their preferential expression in the phloem and cambial tissues (Figure 6c).      (Figure 6b) and show a higher similarity to AtHAM4 compared to other AtHAMs, which were found from both the amino acid sequence alignment ( Figure S3a) and phylogenic analysis ( Figure S3b). AtHAM4 was recently reported as a candidate for regulating cambium initiation in Arabidopsis by interacting with WOX4 [19]. Thus, we designated Potri.005G125800.1 and Potri.007G029200.1 as PtrHAM4-1 and PtrHAM4-2, respectively. The high spatial-resolution wood formation data [37] further support their preferential expression in the phloem and cambial tissues (Figure 6c).  Table 2. PtrHAM4-1 (Potri.005G125800) was relocated to emphasize. (b) Tissue-specific expression of both Potri.005G125800.1/PtrHAM4-1 and its closest homolog, Potri.007G029200.1/PtrHAM4-2. This diagram was reconstructed from our RNAseq data. (c) PtrHAM4-1 is highly expressed in the phloem and cambial tissues. To obtain a gene expression profile by exploiting the high spatial-resolution wood formation data [37] the list of genes from Table 2 was queried to the AspWood website. The resulting heatmap showed that PtrHAM4-1 (indicated by red letters) is highly expressed in the phloem and cambial tissues.  Table 2.
PtrHAM4-1 (Potri.005G125800) was relocated to emphasize. (b) Tissue-specific expression of both Potri.005G125800.1/PtrHAM4-1 and its closest homolog, Potri.007G029200.1/PtrHAM4-2. This diagram was reconstructed from our RNAseq data. (c) PtrHAM4-1 is highly expressed in the phloem and cambial tissues. To obtain a gene expression profile by exploiting the high spatial-resolution wood formation data [37] the list of genes from Table 2 was queried to the AspWood website. The resulting heatmap showed that PtrHAM4-1 (indicated by red letters) is highly expressed in the phloem and cambial tissues.

Overexpression of Ptrham4-1 Enhanced Vascular Cambium Development in Both Transgenic Arabidopsis and Poplar Plants
To test whether PtrHAM4-1 is functionally involved in cambium formation, we generated transgenic Arabidopsis plants overexpressing PtrHAM4-1 using the CaMV 35S promoter (i.e., 35S::PtrHAM4-1). Overall, the 35S::PtrHAM4-1 Arabidopsis plants grew normally, but with a dramatic increase of cambial cell proliferation in the inflorescent stems of all five independent T3 homozygous lines (Figure 7). We observed the secondary xylem vessels in the interfascicular region of the 35S::PtrHAM4-1 Arabidopsis plants, which are not found in wild type (WT) plants (Figure 7c). This fact indicates that the vascular cambium was developed in the 35S::PtrHAM4-1 Arabidopsis plants. Accordingly, the 35S::PtrHAM4-1 Arabidopsis exhibited significantly increased interfascicular cambium-derived tissues (ICD) extension of >60% compared to WT plants (Figure 7c,d).
homozygous lines (Figure 7). We observed the secondary xylem vessels in the interfascicular region of the 35S::PtrHAM4-1 Arabidopsis plants, which are not found in wild type (WT) plants (Figure 7c). This fact indicates that the vascular cambium was developed in the 35S::PtrHAM4-1 Arabidopsis plants. Accordingly, the 35S::PtrHAM4-1 Arabidopsis exhibited significantly increased interfascicular cambium-derived tissues (ICD) extension of >60% compared to WT plants (Figure 7c and d).  The same construct (35S::PtrHAM4-1) was introduced to a hybrid poplar to further confirm the phenotypic significance of 35S::PtrHAM4-1 Arabidopsis plants. We generated a total of 25 independent transgenic poplar lines and found the increased cambial development (e.g., ICD extension) in many lines compared to the wild-type BH clone (i.e., BH) (Figure 8a). Subsequent quantification and gene expression analysis further showed that the phenotypic significance is nicely correlated with the expression level of the PtrHAM4-1 gene (Figure 8). Taken together, our results suggest that PtrHAM4-1 may function as an important player in the initiation of vascular cambium formation in poplar.
independent transgenic poplar lines and found the increased cambial development (e.g., ICD extension) in many lines compared to the wild-type BH clone (i.e., BH) (Figure 8a). Subsequent quantification and gene expression analysis further showed that the phenotypic significance is nicely correlated with the expression level of the PtrHAM4-1 gene (Figure 8). Taken together, our results suggest that PtrHAM4-1 may function as an important player in the initiation of vascular cambium formation in poplar.

Discussion
Secondary growth is one of the most important biological processes on earth. However, our knowledge concerning the underlying molecular mechanisms of vascular cambium initiation/proliferation and vascular patterning is still fragmented, especially for woody perennials. We hypothesized that genes abundantly and specifically expressed in wood-forming tissues may be important players in wood formation. Recently, many extensive wood-forming tissue-type transcriptome analyses of Populus have been reported. Shi et al. [36] in their study used four tissues (shoot tip, leaf, xylem, and phloem) and two wood-forming cell types (fibers and vessels) of P. trichocarpa. Sundell et al. [37] presented high spatial-resolution transcriptome data spanning the phloem, vascular cambium, early/developing xylem, and mature xylem with successive longitudinal tangential cryosections of P. tremula stem. Chao et al. [38] used three tissues from a hybrid poplar (Populus deltoids x P. euramericana cv. 'Nanlin895'): Shoot apex; internodes 1-3 (IN1-3); and internodes [4][5].
Here, we emphasized the transcriptome analysis of cambium differentiation and xylem cell fate specification by designing a gradient of secondary growth developmental stages, including cambium and xylem cell types, within a total of seven tissues: SL; IS; IC; IDX; MC; MDX; and L ( Figure 1a). To overcome the limitations of our single run transcriptome data (i.e., no biological replication), we tried; (1) to minimize the biological variation of the samples by combining tissue samples from at least 10 poplars; (2) to check the reliability of our wood-forming tissue-specific transcriptome data by employing the sample correlation heatmap and matrix data ( Figure 1c); (3) to validate our transcriptome data using well-known tissue-specific marker genes ( Figure 2). We performed not only an extensive transcriptome analysis but also a proof of concept experiment to investigate the molecular mechanism(s) underlying wood formation in poplar.

Critical Pathway Genes of Secondary Wall Biosynthesis in Woody Perennials
Utilizing our series of poplar wood-forming tissue transcriptomes, we sought to uncover the key players for each step in the biosynthesis pathways of secondary wall components (i.e., cellulose, xylan, and lignin) (Figures 3-5). Synthesis of secondary wall components depends on the import of sucrose and its subsequent metabolism. However, our understanding of cell wall precursor biosynthesis in developing wood is still limited [24].
Cellulose is a major constituent of plant cell walls and provides load-bearing strength by forming scaffolds with other cell wall polymers, such as xylan and lignin. Arabidopsis has six SUS (SUCROSE SYNTHASE) genes; SUS is believed to be a main route of carbon entry from sucrose to cellulose via production of UDP-Glc, a substrate for the cellulose synthase [25]. Of the seven SUS genes within our transcriptome data, three are highly expressed in MDX compared to expression in L tissue (Figure 3), indicating that these genes are important for cellulose biosynthesis in xylem. Among them, Potri.006G136700, the closest homolog to Arabidopsis SUS4 (AT3G43190), can be regarded as a key player as it exhibits the highest expression and specificity in MDX. An alternative pathway to production of UDP-Glc from sucrose is through invertases and UGP (UDP-GLC PYROPHOSPHORYLASE). However, UGP is not known to be rate-limiting for cell wall synthesis [59]. We found the gene Potri.013G110800 within the 13 invertases in our data, which is a homolog of Arabidopsis CINV2 (CYTOSOLIC INVERTASE) (AT4G09510) with a strong specificity to xylem. It has been reported that loss of cytosolic invertase affects cell wall synthesis in Arabidopsis [60]. Accordingly, Rende et al. [61] suggested that Potri.013G110800 (reported as CIN12) is responsible for supplying UDP-Glc for cellulose biosynthesis in the development of wood of a hybrid aspen (Populus tremula × tremuloides). Cellulose is produced at the cell surface by CESAs (CELLULOSE SYNTHASE A) [62]. Out of the 13 CESAs in our data, 10 exhibited higher expression in MDX, including: Potri.002G257900.1; Potri.011G069600.1; and Potri.006G181900.1, which are homologous to the well-known secondary wall CESAs; CESA4 (AT5G44030), CESA8 (AT4G18780), and CESA7 (AT5G17420), respectively [63].
Xylans contain xylose subunits as a backbone and are one of the major hemicelluloses found in the secondary cell walls of poplar. Unlike cellulose, xylan has reducing end oligosaccharides, which may act as primers or terminators as well as a variety of side chains [24,26]. UGD (UDP-GLUCOSE DEHYDROGENASE) is responsible for directing xylan biosynthesis from UDP-Glc and for catalyzing the production of UDP-glucuronate [64] (Figure 3). Among the four UGDs within our transcriptome data, three were more highly expressed in MDX-including Potri.008G094300.2-which is highly specific to MDX. UXS (UDP-XYLOSE SYNTHASE) utilizes UDP-glucuronate as a substrate to produce UDP-xylose, a subunit of the xylan backbone. Both Potri.010G207200.5 and Potri.001G237200.3 have higher expression and specificity within MDX, thus, these genes are likely the major UXS genes catalyzing this step. The xylosyl backbone is synthesized by the IRX9, IRX10, and IRX14 genes [65]. Our data suggest that Potri.006G131000.1, Potri.001G068100.3, and Potri.001G067500.1 are key players and are homologous to IRX9, IRX10, and IRX14, respectively. In addition, the poplar genes homologous to those in Arabidopsis that synthesize the reducing end oligosaccharide of xylan, such as IRX7/FRA8, IRX8, and GATL/PARVUS (GALACTURONOSYL TRANSFERASE-like), and those responsible for xylan modification, such as GUX, DUF579/GXM, and RWA, are listed in Figure 3. Based on MDX expression and specificity compared to L tissue, we can predict the essential players of each step.
Lignin is a phenolic compound providing compression strength and hydrophobicity to cell walls, and is a highly heterogeneous and complex polymer composed of p-hydroxylphenyl (H), guaiacyl (G), and syringyl (S) units. The lignin monomer biosynthetic pathway has been extensively investigated and well described [27]. In this pathway, 4CL (4-COUMARATE:CoA LIGASE) catalyzes to produce p-coumaroyl-CoA, a precursor for the biosynthesis of all three of the aforementioned monomers. Potri.001G036900.1 (similar to Arabidopsis 4CL2, AT3G21240.1) seems to be a key gene for this step (Figure 4). S units are synthesized by the action of F5H (FERULATE 5-HYDROXYLASE) and COMT (CAFFEIC ACID O-METHYLTRANSFERASE). Potri.007G016400.1 and Potri.012G006400.1 (similar to Arabidopsis F5H, AT4G36220.1 and COMT, AT5G54160.1, respectively) are prominent within these steps in poplar (Figure 4). Recently, the CSE (CAFFEOYL SHIKIMATE ESTERASE, AT1G52760.1) enzyme has been added to the lignin pathway of Arabidopsis. CSE provides an alternative route to caffeoyl-CoA by catalyzing caffeoyl shikimate [6]. The poplar genome has two CSEs (Potri.001G175000.1 and Potri.003G059200.1) and both exhibit high expression in MDX ( Figure 4). Since Arabidopsis cse mutants deposit less lignin but display a four-fold increase in saccharification yield without pretreatment [6], CSE is a promising target for the development of improved lignocellulosic biomass. Recently, Saleme et al. [66] reported that silencing CSE in a hybrid poplar that is morphologically indistinguishable from WT poplar resulted in up to 25% and 62% reduced lignin deposition and increased glucose release, respectively, without pretreatment.
This analysis was able to identify all the essential genes within each step of the biosynthesis pathways for secondary wall components (Figures 3-5), with these genes being potential targets for biotechnological improvement for the purpose of ascertaining quality of woody biomass. Very recently, Wang et al. [67] reported a multi-omics quantitative integrative analysis of lignin biosynthesis by perturbing 21 pathway genes to advance the strategic engineering of wood utilization, and these 21 target genes are exactly matched to key players identified from our analysis, with the exception of CSE (Figure 5b). A similar approach could be applicable to other pathways (i.e., those for cellulose and xylan) in combination in order to produce desirable woody biomass.

Transcriptional Regulators Involved in Vascular Cambium Development
Previously, Pineau et al. [68] reported that the hca (high cambium activity) mutant showed premature and numerous cambial cell divisions in both the fascicular and interfascicular regions.
To determine the influence of the hca mutation on global gene expression they performed transcriptome profiling and found AtHAM4/SCL15 (At4g36710), the most upregulated transcription factor belonging to a member of the HAM family, which plays an essential role in shoot meristem maintenance in a non-cell-autonomous manner [20,21]. Zhou et al. [19] demonstrated that HAM family members act as conserved interacting cofactors with WUS/WOX (WUSCHEL/WUSCHEL-RELATED HOMEOBOX) proteins. WUS is a homeodomain transcription factor expressed in the rib meristem of the SAM (shoot apical meristem), and is a key regulatory factor controlling stem cell populations in Arabidopsis [69]. WOX4, expressed in Arabidopsis procambial cells, defines the vascular stem cell niche and regulates cambial cell proliferation [15,18,70]. In particular, Zhou et al. [19] showed that HAM4 and WOX4 physically interact in vivo and were tightly co-regulated in both a spatial and temporal manner. For example, HAM4 and WOX4 are co-expressed in provascular or procambial cell types of various tissues and, in the stem transverse section, HAM4 is expressed specifically in the procambium and overlaps with WOX4 expression. Recently, Kucukoglu et al. [71] suggested that WOX4-like genes regulate cambial cell division activity and secondary growth by using PttWOX4a/b RNAi in Populus trees.
In our co-expression network analysis (Figure 6a), we found several transcription factors highly expressed in IC tissue. Among them, we focused on PtrHAM4-1, a homolog of AtHAM4, which is specifically expressed in the cambial tissues ( Figure 6). Indeed, overexpression of PtrHAM4-1 in transgenic Arabidopsis plants resulted in a dramatic increase of the vascular cambium with clear differentiation of secondary xylem vessels, which were known to be produced only from the vascular cambium (Figure 7). This phenotype was further confirmed by overexpressing the PtrHAM4-1 in hybrid poplars (Figure 8). To our knowledge, this is the first evidence that PtrHAM4-1 acts positively in vascular cambium development in plants. Taken together with a report from Kucukoglu et al. [71], our results suggest that a HAM4-WOX4 regulatory module may be conserved in Populus to achieve cambial meristem initiation and proliferation during secondary growth. However, the underlying molecular regulatory mechanism remains unknown. Chromatin modifications, including histone acetylation, have been implicated in meristem activity. Thus, interaction of HAM4/SCL15 with HDA19 (HISTONE DEACETYLASE19) [72] may provide a clue for the regulation of downstream gene expression, which will be a priority for our future studies.
In summary, a comprehensive wood-forming tissue-specific transcriptome analysis from a hybrid poplar successfully pinpointed many essential genes involved in the biosynthetic pathways of secondary wall components. These genes could be focal points for the biotechnological improvement of wood properties within the production of biomaterials and/or biofuels. Furthermore, the transcriptional regulators involved in vascular cambium development were isolated and demonstrated their validity via functional characterization of PtrHAM4-1 using a heterologous expression. Thus, our results may offer insights for disentangling the complex mechanisms of wood formation, one of the most important biological processes on this planet.

Supplementary Materials:
The following are available online at http://www.mdpi.com/2073-4425/10/9/690/s1, Figure S1: RNA-seq raw data processing, Figure: S2. Mapping of the sequence reads to the genome of Populus trichocarpa, Figure S3: PtrHAM4-1 is homologous to AtHAM4, Table S1: A total of 38,330 transcripts across all tissues, after removing redundant transcripts from a total of 85,209 transcripts, Table S2: List of primers used in this study, Table S3: List of transcripts preferentially expressed in IC (885 transcripts), Table S4: List of transcripts preferentially expressed in MC (798 transcripts).