Global Transcriptomic Profile Analysis of Genes Involved in Lignin Biosynthesis and Accumulation Induced by Boron Deficiency in Poplar Roots

To uncover the transcriptomic mechanism of lignin accumulation caused by boron deficiency (BD), Nanlin895 (Populus × euramericana “Nanlin895”) was subjected to control (CK, 0.25 mg·L−1) and BD (0 mg·L−1) treatments for 3 days. RNA-Seq was carried out to survey the expression patterns of the lignin-regulated biosynthetic genes in response to BD. The results showed that 5946 genes were identified as differentially expressed genes (DEGs), 2968 (44.2%) of which were upregulated and 3318 (55.8%) of which were downregulated in response to BD. Among them, the expression of lignin monomer biosynthetic (PAL, CCR, CAD, COMT, F5H, PER/LAC) and modulated genes, for example, transcription factors (MYBs) and hormone signal regulating genes (GIDs, histidine kinase 1, coronatine-insensitive protein 1), were upregulated, and some hormone signal regulating genes, such as AUXs and BR-related (sterol methyltransferases), were downregulated under BD treatment. There are also some genes that were screened as candidates for an association with wood formation, which will be used for the further analysis of the function of lignin formation. These results provide an important theoretical basis and reference data in plant for further research on the mechanism of lignin accumulation under BD.


Introduction
Boron is one of the essential elements required for normal growth in higher plants [1]. It plays important roles in the physiological processes of cell wall structure and function, phenol metabolism, the structure and function of membranes, sugar transportation, carbohydrates, and nucleic acid metabolism [2][3][4]. Studies have shown that boron and the pectin polysaccharide RGII combine to form a dipolymer, which is involved in cell wall formation [5,6]. Boron deficiency (BD) results in the accumulation of pectinolytic-related substances and lignin in mature cells, inhibiting the formation and development of new tissues.
BD is a common phenomenon on a global scale, which is mainly caused by the leaching of rains in the soil [7]. Sillanpaa et al. [8] studied 190 samples of agricultural soil, and the results indicated that BD led to a 31% decrease in the yield of crops. With the emphasis on the harm and damage caused by BD, many studies have focused on crops, such as Arabidopsis [2], wheat [9], maize [10] and citrus [11,12]. From these studies, some genes, such as EXP14, EXPB1, XTHs, CSLB5, PG, and PME, related to the synthesis of the cell wall have been discovered under BD conditions [3]. However, until now, few researchers have focused on the cell wall synthesis of poplars under BD.
BD also leads to a decrease in forestry productivity. Populus is one of the most widespread trees that is planted all over the world because of its characteristics of fast growth, high yield, and strong adaptability [13,14]. It is a kind of raw material for the paper and pulp industries [15]. Lignin, a kind of phenol polymer, is an important component of the cell wall in vascular plants, and the content of lignin is an essential factor that influences the properties of the production of paper with poplars. In recent years, many researchers have made substantial efforts to breed planting clones with a lower lignin content; for example, some genetically modified poplars with a lower lignin content were obtained [16,17]. Thus, studies on the molecular mechanism of lignin accumulation after abiotic stresses (such as BD) are needed to breed low lignin poplar.
Previous studies by our group revealed that BD increased the lignin in citrus vessel cell walls [18]. In this paper, we propose that BD will cause the accumulation of lignin in poplar. To understand the molecular mechanism of lignin accumulation in poplars, poplar saplings were hydroponically cultured in nutrient solution with BD. Digital gene expression was performed to dissect the gene reprogramming profiles underlying lignin accumulation under BD.

Plant Materials and Treatment
Nanlin 895 (P. × euramericana "Nanlin895") used in this study were from the poplar germplasm nursery of Shishou City in Hubei Province, China. In early July, two-month-old culture saplings approximately 30-40 cm in height with no pests, good growth, and uniform size were taken out from the pots, cleaned up, and precultured in nutrient solution (1/2 Hoagland's nutrient solution and iron minor elements) for two weeks and cultured with deionized water for three days before BD treatment. During the treatment, the plants were divided into two groups randomly; a total of 168 poplar saplings with 12 independent biological replicates for each stage were included in each group and treated for different times (0, 1, 3, 5, 7, 10, and 15 days). Saplings were cultured using Hoagland's nutrient solution containing either 0.25 mg·L −1 boron (control, CK) or 0 mg·L −1 boron (boron deficiency, BD). Poplar saplings were cultured in black plastic buckets ventilated for 15 min every two hours, and the nutrient solutions were replaced weekly. Root tips were collected to determine the content of lignin at different time points (0, 1, 3, 5, 7, 10, and 15 days) with twelve biological duplications. For the saplings, we used three plants to determine the lignin content, three for RNA-Seq, three for qRT-PCR, and others as additional materials. Then, the root tips with the significantly higher lignin content were selected to extract RNA using random selection method. All the nutrient solutions used in this experiment were prepared with deionized water.

Determination of Lignin Content
For lignin content determination, the fresh samples were heated for 20 min at 105°C, then heated at 65°C, and maintained at a constant weight. Then, the dried samples were ground into powder. Two methods (Klason and ultraviolet spectrophotometry (UV-6100, Mapada, Shanghai, China) were used to determine the lignin content. For the calculation of standard light absorption (A standard ), Klason methods [19] were used to estimate the lignin content of the untreated samples. Absorptivity was calculated as Absorptivity = Abs × liters/W samples × k lignin Astandard were calculated as average number of different Absorptivity of biological replicates. For the lignin content of different stages of root samples, ultraviolet spectrophotometry was carried out [20]. Briefly, ca. 6 mg of root powder samples were mixed with 5 mL of 25% (w/w) acetyl bromide acetic acid and 0.2 mL of perchloric acid, sequentially. The tubes were sealed and incubated at 70°C for 30 min, and then the mixture was mixed with 10 mL of 2 M NaOH and diluted with glacial acetic acid to 100 mL. The lignin contents in the samples were determined spectrophotometrically at 254 nm (the maximum absorption wavelength was obtained by scanning the absorption peak using ultraviolet-visible spectroscopy in the range of 200-600 nm). The calculation of the lignin contents was performed via the formula lignin % = Abs × liters × 100%/W samples × A standard .

Total RNA Extraction and Digital Gene Expression Profiles
Root tips (ca. 3.5 cm from the top) at the stage of day 3 were collected for the following RNA isolation. Total RNA was extracted from the root tips of CK and BD treatments with 3 biological replicates at the same phase. RNA extraction was performed according to the manufacturer's instructions of an RNA Plant Reagent (TIANGEN, Beijing, China).
RNA samples extracted from the root tips were sent to Shenzhen Huada Genomics Co (Huada, Shenzheng, China), Ltd. for quality and quantity evaluation, RNA-Seq library construction, and high throughput sequencing. In brief, RNA-Seq library construction was performed according to the operating instructions (MGIEasy RNA Library Prep Set, version 3.0 (Huada, Shenzheng, China). In brief, mRNA was enriched by using Oligo(dT) magnetic beads (Dynabeads®mRNA Purification Kit; Invitrogen, Cat. No. 61006, Carlsbad, CA, USA) from total RNA under high-temperature conditions. Then, the mRNA was broken into short fragments by the addition of fragmentation buffer. Then, the short fragments of mRNA were used as templates for the synthesis of first strand cDNA and second strand cDNA. The double string cDNA purified with magnetic beads was subjected to repair of the 3' end, the addition of a dNTP base, and the ligation of sequencing adapters. The fragments were enriched using PCR amplification, and the products were purified using agarose gel electrophoresis. The quantity and quality of the completed cDNA libraries were tested using an Agilent 2100 Bioanalyzer (Agilent, Santa Clara, CA, USA) and an ABI Step One Plus Real-Time PCR System (Thermo Fisher Scientific, Santa Clara, CA, USA). The qualified library was sequenced with Illumina HiSeq TM2000 (Illumina, San Diego, CA, USA).
The clean reads were mapped to reference sequences of P. trichocarpa (ftp://ftp.jgi-psf.org/pub/ compgen/phytozome/v9.0/Ptrichocarpa/annotation) SOAPaligner/SOAP2. Sequence saturation was used to analyze the number of genes detected that tended to saturate. The distribution of reads on the reference genes was used to evaluate the randomness.

Identification of Differentially Expressed Genes (DEGs)
The gene expression level was calculated by using the RPKM method (Reads Per kb per Million reads), as suggested by Wagner et al. [20].
The RPKM method can eliminate the influence of different gene sequencing discrepancies on the calculation of gene expression levels. Therefore, the RPKM values can be directly used for comparing the differences in gene expression among samples. The absolute value of log 2 (BD-RPKM/CK-RPKM) was used to represent the fold change of two genes with a threshold FDR ≤ 0.001, and the absolute value of |log 2 Ratio| ≥ 1 was used to judge the significance of the difference in gene expression.
For the analysis of function and pathway, all the DEGs were mapped to GO terms in the database (http://www.geneontology.org/) and pathway terms in the KEGG database (http://www.genome.jp/ kegg/), as described by Ye et al. [21] and Minoru et al. [22], respectively. GO term enrichment analysis was conducted in the agriGO database with the singular enrichment analysis (SEA) tool as described by Du et al. [23,24]. For functional categories, DEGs were submitted to MapMan software (version 3.6.0RC1, Golm, Germany) according to the standard protocol with minor modifications [24].

Quantitative Real-Time PCR Analysis
qRT-PCR analysis was used to verify the DEG results. Nine differentially expressed candidate genes involved in the biosynthetic pathway of lignin monomers and regulatory pathways were chosen for quantitative real-time PCR. For this analysis, gene-specific primers were designed according to the reference gene sequences using Primer 5.0 (Canada), and qRT-PCR was performed according to the ZHUANGMEN manufacturer specifications (Product code: Catalog # ZF101-102, Beijing, China). Briefly, a 20 µL reaction mixture included 10 µL 2× SYBR qPCR Mix, 0.4 µL Dye ROX, 0.5 µL cDNA temple, 0.4 µL forward primer (10 µM), 0.4 µL reverse primer (10 µM), and ddH 2 O to a final volume, and the reaction was performed on the 7500 Fast ABI Real-time PCR system (Applied Biosystems, Foster City, CA, USA) using a two-step method: 94 • C for 3 min, 94 • C for 15 s, and 60 • C for 1 min. The Poplar actin2 gene was used as an internal reference [25], and the relative expression levels of the genes were presented by 2 -∆∆CT . The experiment was conducted with three biological reduplicates and three technological repeats. The primers used in this study are displayed in Table S1.

Statistical Analysis
Statistical analysis was carried out with SPSS 17.0 (SPSS, Chicago, IL, USA). Data were tested for normal distribution, and one-way analysis of variance (ANOVA) was used to compare the differences between means. p-values lower than 0.05 were considered significant.

Effect of BD on the Lignin Content in Poplar Roots
The root lignin content of different stages of "Nanlin895" under treatment and nontreatment conditions were determined using ultraviolet spectroscopy at a 254 nm wavelength. The lignin content gradually increased in the early stage of the treatment, peaked on day 3, and tended to remain constant for the following days ( Figure 1). For the control, the lignin content was the lowest (28.42%) on day 3 and the highest (32.75%) on day 15 ( Figure 1). In contrast, it was the highest on day 3 (39.28%), and the average content was remarkably higher than that of the control under boron deficiency conditions ( Figure 1).  [25], and the relative expression levels of the genes were presented by 2 -ΔΔCT . The experiment was conducted with three biological reduplicates and three technological repeats. The primers used in this study are displayed in Table S1.

Statistical Analysis
Statistical analysis was carried out with SPSS 17.0 (SPSS, Chicago, IL, USA). Data were tested for normal distribution, and one-way analysis of variance (ANOVA) was used to compare the differences between means. p-values lower than 0.05 were considered significant.

Effect of BD on the Lignin Content in Poplar Roots
The root lignin content of different stages of "Nanlin895" under treatment and nontreatment conditions were determined using ultraviolet spectroscopy at a 254 nm wavelength. The lignin content gradually increased in the early stage of the treatment, peaked on day 3, and tended to remain constant for the following days ( Figure 1). For the control, the lignin content was the lowest (28.42%) on day 3 and the highest (32.75%) on day 15 ( Figure 1). In contrast, it was the highest on day 3 (39.28%), and the average content was remarkably higher than that of the control under boron deficiency conditions ( Figure 1).

Read Mapping and Analysis of DEGs
To investigate the changes in gene expression in transcription in poplars under boron deficiency, poplar saplings were cultured in Hoagland's nutrient solution containing either 0.25 mg·L −1 boron (control, CK) or 0 mg·L −1 boron (boron deficiency, BD). Root tips with the highest contents of lignin were used to extract total RNA and were further utilized to construct libraries (CK and BD); high-throughput sequencing was performed using Illumina HiSeq TM2000. More than 11 million clean reads, above 93% of the total raw reads, were obtained from each library ( Figure S1). More than 7 million clean reads from each library were mapped to reference genes, and approximately 78% of the clean reads were mapped to the reference genome (Table S2). Between BD and CK, 5944 genes had an FDR ≤ 0.001 and |log 2 Ratio| ≥ 1 ( Figure 2, Table S3). Among these significantly expressed genes, 2628 (45.1%) genes were upregulated and 3318 (54.9%) genes were downregulated in response to BD stress ( Figure 2, Table S3).

Read Mapping and Analysis of DEGs
To investigate the changes in gene expression in transcription in poplars under boron deficiency, poplar saplings were cultured in Hoagland's nutrient solution containing either 0.25 mg·L −1 boron (control, CK) or 0 mg·L −1 boron (boron deficiency, BD). Root tips with the highest contents of lignin were used to extract total RNA and were further utilized to construct libraries (CK and BD); highthroughput sequencing was performed using Illumina HiSeq TM2000. More than 11 million clean reads, above 93% of the total raw reads, were obtained from each library ( Figure S1). More than 7 million clean reads from each library were mapped to reference genes, and approximately 78% of the clean reads were mapped to the reference genome (Table S2). Between BD and CK, 5944 genes had an FDR ≤ 0.001 and |log2Ratio| ≥ 1 (Figure 2, Table S3). Among these significantly expressed genes, 2628 (45.1%) genes were upregulated and 3318 (54.9%) genes were downregulated in response to BD stress ( Figure 2, Table S3). Figure 2. Number of significantly differentially expressed genes in "Nanlin895" roots exposed to boron deficiency. GO term analysis was performed; DEGs were classified into three categories, i.e., biological process, cellular component, and molecular function, and significantly enriched in 13 GO terms ( Figure 3, Table S4). This analysis revealed five biological processes, including microtubule-based process (GO:0007017), microtubule-based movement (GO:0007018), metabolic process (GO:0008152), carbohydrate metabolic process (GO:0005975), and photosynthesis (GO:0015979); three molecular functions, including microtubule motor activity (GO:0003777), transcription factor activity (GO:0003700), motor activity (GO:0003774); and five cellular components, including chromosomal parts (GO:0044427), chromosome (GO:0005694), chromatin (GO:0000785), nucleosomes (GO:0000786), and protein-DNA complexes (GO:0032993), were significantly enriched in the 13 GO terms (Figure 3, Table S4). Number of significantly differentially expressed genes in "Nanlin895" roots exposed to boron deficiency. GO term analysis was performed; DEGs were classified into three categories, i.e., biological process, cellular component, and molecular function, and significantly enriched in 13 GO terms ( Figure 3, Table S4). This analysis revealed five biological processes, including microtubule-based process (GO:0007017), microtubule-based movement (GO:0007018), metabolic process (GO:0008152), carbohydrate metabolic process (GO:0005975), and photosynthesis (GO:0015979); three molecular functions, including microtubule motor activity (GO:0003777), transcription factor activity (GO:0003700), motor activity (GO:0003774); and five cellular components, including chromosomal parts (GO:0044427), chromosome (GO:0005694), chromatin (GO:0000785), nucleosomes (GO:0000786), and protein-DNA complexes (GO:0032993), were significantly enriched in the 13 GO terms (Figure 3, Table S4).  . The significantly enriched GO terms for differentially expressed genes in "Nanlin895" roots exposed to boron deficiency. The x-axis and y-axis indicate the names of clusters and gene number assigned to the cluster, respectively.
Pathway analysis of the DEGs based on the KEGG database was performed to characterize the functional consequences caused by BD. The results certified that there are 11 pathways significantly enriched in response to BD, in which genes participating in metabolic pathways (26.65%), plant hormone signal transduction pathways (8.92%), and starch and sucrose metabolism (3.47%) were the top three pathways enriched under BD (Table 1). Other pathways, such as nitrogen metabolism; glycine, serine, and threonine metabolism; and ascorbate and aldarate metabolism were also enriched in response to BD (Table 1). Table 1. Significantly enriched KEGG pathways of significantly differentially expressed genes in roots of "Nanlin895" exposed to boron deficiency.

Candidate Genes Associated with Lignin, Cellulose, and Flavonol Synthesis Under BD
Lignin is biosynthesized from the phenylalanine pathway. Here, the candidate genes involved in lignin biosynthesis were identified as flow and regulation pathways. The results of global gene expression showed that numerous candidate genes were affected in the roots of "Nanlin 895" under BD conditions (Tables 2,3,S5). Some were upregulated, and these genes included PAL, CCR, CAD, COMT, F5H, PER, and LAC ( Table 2,S5). Additionally, a total of six identified candidate transcription factors (MYB(2), TALE (2), and WRKY (2); Tables S3,S5) were upregulated. Some candidate genes also participated in hormone biosynthesis or signal transduction during cellular differentiation, including  (Table 1). Other pathways, such as nitrogen metabolism; glycine, serine, and threonine metabolism; and ascorbate and aldarate metabolism were also enriched in response to BD (Table 1). Table 1. Significantly enriched KEGG pathways of significantly differentially expressed genes in roots of "Nanlin895" exposed to boron deficiency.

Candidate Genes Associated with Lignin, Cellulose, and Flavonol Synthesis Under BD
Lignin is biosynthesized from the phenylalanine pathway. Here, the candidate genes involved in lignin biosynthesis were identified as flow and regulation pathways. The results of global gene expression showed that numerous candidate genes were affected in the roots of "Nanlin 895" under BD conditions (Tables 2 and 3 and Table S5). Some were upregulated, and these genes included PAL, CCR, CAD, COMT, F5H, PER, and LAC (Table 2 and Table S5). Additionally, a total of six identified candidate transcription factors (MYB(2), TALE(2), and WRKY (2); Tables S3 and S5) were upregulated. Some candidate genes also participated in hormone biosynthesis or signal transduction during cellular differentiation, including gibberellin receptors (GID1/2), gibberellin oxidases (gibberellin 20-oxidase/gibberellin 2-oxidase), sterol methyltransferase, chitinases (ethylene response factor), and some auxin-responsive proteins (Tables S3 and S5). These different expression patterns of transcription factors and signal transduction response genes are involved in the regulation of the dynamic accumulation of lignin during root development. Table 2. Selected significantly differentially expressed genes related to phenylpropanoid biosynthesis and cellulose biosynthesis in roots of "Nanlin895" exposed to boron deficiency.  Table 3. Significantly differentially expressed transcription factors and phytohormone related genes in roots of 'Nanlin895' exposed to boron deficiency. Flavonols were also synthesized from the phenylalanine pathway, and some of the major enzymes were identified. Interestingly, the CHS and F3'5'H genes identified were all downregulated ( Figure 4; Table 2). This suggests that the deduced flavonol synthesis pathway may promote the upregulation of lignin. Cellulose-associated genes were also identified in this study, and five genes (hexokinase (3) and cellulose synthase (2)) were also downregulated ( Table 2).

Gene ID
Flavonols were also synthesized from the phenylalanine pathway, and some of the major enzymes were identified. Interestingly, the CHS and F3'5'H genes identified were all downregulated ( Figure 4; Table 2). This suggests that the deduced flavonol synthesis pathway may promote the upregulation of lignin. Cellulose-associated genes were also identified in this study, and five genes (hexokinase (3) and cellulose synthase (2)) were also downregulated ( Table 2).  Table 2. Selected significantly differentially expressed genes related to phenylpropanoid biosynthesis and cellulose biosynthesis in roots of "Nanlin895" exposed to boron deficiency.

Validation of RNA-Seq Results Using qRT-PCR
To verify the results of the DEGs, nine DEGs were selected to analyze their expression levels using qRT-PCR in two libraries ( Figure 5). The selected genes contained lignin monomer synthesis genes and transcription factors (Table S1). The relative expression levels of CK and BD were comparable with those of the RNA-Seq data (Figure 4). Although the results of qRT-PCR were different from those of RNA-Seq, they all shared a similar direction of change between qRT-PCR and RNA-Seq (Figure 4). This comparison validated the results of RNA-Seq as accurate.

Validation of RNA-Seq Results Using qRT-PCR
To verify the results of the DEGs, nine DEGs were selected to analyze their expression levels using qRT-PCR in two libraries ( Figure 5). The selected genes contained lignin monomer synthesis genes and transcription factors (Table S1). The relative expression levels of CK and BD were comparable with those of the RNA-Seq data (Figure 4). Although the results of qRT-PCR were different from those of RNA-Seq, they all shared a similar direction of change between qRT-PCR and RNA-Seq (Figure 4). This comparison validated the results of RNA-Seq as accurate.  Table S1. Solid circles and histograms represent the results obtained from RNA-Seq and qRT-PCR, respectively. Bars indicate the mean ± SE (n = 3). Detailed information on the selected genes is available in Table S1.

BD Caused the Differential Expression of Lignin Biosynthetic Enzyme Genes
The phenylpropanoid pathway is responsible for the synthesis of a large variety of secondary metabolites, including phenol esters, coumarins, flavonoids, and lignin. Boron starvation has been found to significantly affect phenolic metabolism [26]. The expression levels of CHSs, CHI, and F3'5'H, which function as the rate-limiting enzymes in the flavonoid biosynthetic pathway (Figure 4, Table 2), were downregulated under BD, but this was not the case for the lignin monomer synthetic pathway. It can be hypothesized that the reduction of the flavonoid metabolic pathway compensated for the conversion of intermediate metabolites to the lignin biosynthetic pathway (Tables S2 and S5).
Phenylpropanoids, a key product in the biosynthesis of phenolic compounds, are catalyzed by phenylalanine ammonia-lyase (PAL). Previous reports suggest that PtrPAL2,4,5 is expressed principally in xylem and root tips and mainly determines the production of lignin, which differs from the PtrPAL1,3 function in producing the phenol metabolites [26]. In our investigation, the mRNA levels of PtrPAL4 and PtrPAL5 were upregulated more than 2 folds, which implies excess cinnamic acid promotes upregulation of the downstream monolignol biosynthesis genes and then contributes to the biosynthesis of monolignol in the early time period under BD [27]. The further mechanism still needs to be studied in the future.
Cinnamic acid is catalyzed by a series of enzymes and produces p-coumaroyl-CoA, which is then processed by cinnamoyl-CoA reductase (CCR) to coniferaldehyde, which in turn is converted to coniferyl alcohol by the action of CAD.
The downregulation of PtCAD1 results in a reduction of Klason lignin and an alteration of the proportion of monolignol. This result indicates that PtCAD1 is involved in lignin biosynthesis [28]. In this study, we found that BD induces the upregulation of Ptr CAD1, which may suggest that the accumulation of lignin is associated with its expression. The other CAD that showed upregulation may be induced by other stresses and, as described in a previous study [29], the function can be redundant, which requires additional verification. Although two CCRs were upregulated, we found that CCR2 [30], which was revealed as being involved in the biosynthesis of lignin, was not found, suggesting that the two CCRs may be involved in the lignin biosynthetic pathway.
COMT and F5H are both multifunctional enzymes. F5H functions in changing guaiacyl (G) monolignol to syringyl (S) monolignol and can utilize a large set of substrates, including ferulic acid, coniferaldehyde, and coniferyl alcohol. COMT is responsible for catalyzing O-methylation to a cluster of substrates, including caffeic acid, 5-hydroxy-coniferaldehyde, caffeoyl alcohol, 5-hydroxyl-feruloyl CoA, and 5-hydroxy-coniferyl alcohol [31]. Previously, the downregulation of COMT expression led to a decrease in lignin syringyl monomer content [32]. Consistently, the transcript levels of F5H1 and COMT1/2 were increased under BD, which implied that lignin accumulation was caused by increasing the intermediate product conversion of G and S lignin monomer biosynthesis. We also found that 4CL, C3H, and HCT, the key enzymes of the lignin monomer biosynthetic precursor, decreased during BD, which implied the reduction of G and S lignin monomers and perhaps supported the increase of H:S [33].
Laccases and PER, both oxidases, are responsible for the last step of lignin biosynthesis. Disrupting LAC4 and LAC17 or overexpressing the PRX gene AtPrx17 in Arabidopsis plants reduced or caused an increase in the lignin content, respectively [34,35]. The removal of hydrogen peroxide from the lignin-forming cell culture line of Norway spruce prevented or hindered the formation of lignin [36]. The increased expression levels of LACs and PERs, especially LAC20 [37], CWPO-C [38], and potri.011g120300 and potri.009g042500 homologs with AtLAC17 and AtLAC4 may play key roles in lignin accumulation under BD.

Transcription Factors Involved in the Cell Wall Metabolic Process Under BD
BD can affect the expression of many transcription factors, such as MYBs, WRKYs, and NACs [36]. Our data identified two PtrMYBs that shared a closer correlation with AtMYB58 and AtMYB63, which are the second layer of the regulators of lignin synthesis. These two genes were upregulated in response to BD and may activate downstream monolignin biosynthetic enzyme genes. Some genes (LAC4, LAC17, MYB46/83, and SND2) related to SCW formation coexpressed with BLH6, and the overexpression of BLH6 resulted in an accumulation of lignin deposition, and T-DNA mutants caused a decrease in lignin [33,39]. Our results are in accordance with these previously reported results. Our data demonstrated that two TALE genes (Potri004G159300 and Potri009G120800), which are homologous to BLH6, were upregulated in response to BD. Lignin biosynthesis is controlled by a regulatory cascade of upstream transcription factors controlling the formation of secondary walls by activating several other transcription factors, including MYB43, MYB46, MYB58, MYB63, MYB85, MYB103, SND3, KNAT7, BLH6, and WRKY12 [40][41][42]. Three WRKY genes (Potri001G328000, Potri002G193000, and Potri008G103300) ( Table 2) shared closer correlations with WRKY12, suggesting that these genes play critical roles in regulating the upregulation of monolignin biosynthetic enzyme genes, but the specific roles still need to be studied.

Hormones Participate in the Signal Transduction of Cell Wall Formation Under BD
Plant hormones have been inextricably linked with BD and cell wall formation [43][44][45]. In vascular plants, lignification occurs only in certain types of cells, such as tracheary elements and fiber cells.
To adjust themselves to the changing environment and thicken the cell wall, meristematic cells must receive signals that initiate the cell differentiating and cell wall thickening programs. Changes in phytohormone levels are linked with the synthetic and metabolic pathway genes [46]. Here, we have demonstrated that some GA20OXs and GIDs function as synthases and that the receptors of GA are upregulated; GA2OXs catalyze the deactivation of bioactive GA and are downregulated. Consistent with this, the overexpression of AtGA20ox in transgenic tobacco plants showed a higher content of lignin in the AtGA20ox plants compared to wild-type controls [47,48], implying that the increased levels of GA may cause the accumulation of lignin. In contrast, some results demonstrated that spraying winter wheat with exogenous GA 3 significantly decreased lignin accumulation compared to the control, suggesting that GA interacts with other hormones to function during secondary cell wall formation under boron deficiency [49].
Auxin was considered to take part in the cell expansion and the biosynthesis of the cell wall by repression of the KNOX gene [50,51]. The expression levels of the genes encoding auxin-responsive proteins were found to be downregulated, which means that auxin levels decreased under BD and agreed with the results of Bairu et al. [52]. The mutant iar4 leading to the functional loss of auxin could result in the accumulation of lignin [51]. Ban et al. [53] found that the application of exogenous 2,4-D in high concentrations can inhibit the expression of PAL, which may support the hypothesis that low auxin can oppose the inhibition of PAL expression in different directions. Some other reports showed that BD reduces the elongation of cells but not the division, and three receptor genes of cytokinin were found to be upregulated in the DEGs, which agrees with previous results [54].
Brassinosteroids (BRs) can also protect plants from environmental stresses [55]. Sterol synthetic mutants of sterol methyltransferase 1 showed a deficiency of BRs, cellulose, and the ectopic deposit of lignin but not other cell wall components, such as the content of neutral sugars and pectin [56], suggesting that the downregulation of two sterol methyltransferases decreased the BR content. The reduction of the cellulose content further inhibited cell elongation, which not only compensated for the lignin content but also increased the total content of the cell wall by increasing the superficial area of the cell wall [56].
Jasmonic acid (JA) has been shown to be a significant regulator of plant responses to environmental stresses [46,57]. A reduction in cellulose biosynthesis induces the production of the phytohormones JA and ethylene, resulting in ectopic lignin production [57]. coi1 (impaired in JA-Ile perception) seedlings all exhibit enhanced lignin deposition [57]. With the downregulation of cellulose synthase and upregulation of JA-related genes, it can be considered that cellulose reduction compensates for the lignin content, and to maintain cell proliferation, the increased JA promotes the irregular distribution of lignin under BD.
Similar to the results above, the chitinase-like gene was reported to play a critical role in the biosynthesis of the cell wall, which is regulated by ethylene [58]. Zhong et al. [58] showed that the chitinase-like mutant manifested a reduction of cellulose and the deposition of lignin. Interestingly, we found that half (3 out of 6) of the chitinase genes were upregulated and may take part in other metabolic pathways. The specific functions of these chitinases still merit additional study.
The findings of this study indicated that plant hormones, such as GA, BR, ethylene, auxin, and cytokinin, play important roles in the modulation of lignin biosynthesis under BD. However, their roles in the regulation of lignin synthesis under BD still need to be studied at a deeper level.

Conclusions
In this study, RNA-Seq of the roots of "Nanlin 895" was performed and analyzed. Through the high-throughput sequencing technique, we identified some genes involved in phytohormone (including auxin, BR, GA, ETH, JA, and CTKs) metabolism and signaling, cell elongation (XTH), cell division (CYCD), lignin biosynthesis (including PAL, COMT, F5H, CCR, CAD, LAC and PER), cellulose (CesAs, hexokinase), and flavonol synthesis (CHS, CHI, F3'5'H), as well as transcription factors such as WRKY, TALE, and MYB. These results provide an important theoretical basis and reference data in plant for further research on the mechanism of lignin accumulation under BD.
Supplementary Materials: The supplementary materials for this article are available at http://www.mdpi.com/ 2218-273X/9/4/156/s1. Figure S1. Classification of the raw reads of two libraries. Table S1. Primers used in qRT-PCR. Table S2. Major characteristics of the two libraries. Table S3. Significantly differentially expressed genes. Table S4. GO term analysis of significantly differentially expressed genes. Table S5. Functional categories of significantly differentially expressed genes.