Populus trees propagated from the same clone (Populus deltoides × P. euramericana cv. “Nanlin895”) were grown in the same plantation located at Siyang, Suqian, Jiangsu, China (33° 47′ N, 118° 22′ E). Wood-forming tissues were sampled from 1-m trunk above the breast height (1.3 m from ground) from 2-year-old (formation of JW) and 8-year-old trees (formation of MW) at fast growing time (May 2017). After bark was removed, wood-forming tissue (developing secondary xylem) was collected directly into liquid nitrogen and stored at −80°C freezer for later analysis (Song et al., 2011). Three trees as biological replicates were sampled (Supplementary Figure 1), respectively.

After developing xylem was collected, the tree trunk was used for wood analysis. JW and MW were sampled as illustrated in Supplementary Figure 1. Wood tissue was sectioned into 20 μm in thickness and stained with 0.5% phloroglucinol in 12% HCl. Cross sections were observed under a microscope (Olympus, BX53). The number of fibers and vessels and their cross area were counted using Image J. Meanwhile, the wood cells were separated after treatment using acetic acid/hydrogen peroxide (1:1, v/v) solution at 80°C for 6 h. The separated wood cells were then stained with safranine (1% in water), and the length of fiber cells and vessels was measured under a microscope (Olympus, BX53) using Image J.

Air-dried wood sample was ground into powder and filtered through 60-mesh sieve. According to our previous established protocol (Yu et al., 2014), alcohol-insoluble residues (AIR) were firstly obtained by extracting the wood powder with 70% ethanol, chloroform/methanol (1:1, v/v), and acetone. Amylase and pullulanase in 0.1 M sodium acetate buffer (pH 5.0) were used to treat the extracted AIR overnight. For analysis of the sugar in hemicelluloses, AIR was treated with 2 M trifluoroacetic acid (TFA) at 121°C for 90 min. The supernatant was evaporated and incubated in 20 mg/ml fresh sodium borohydride solution at 40° for 90 min. The product was then neutralized with acetic acid and mixed with 1-methylimidazole and acetic anhydride for acetylation. After extraction with dichloromethane, the product was mixed with ethyl acetate for GC-MS (6890N GC system and 5975 Mass detector, Agilent Technologies, equipped with a SP-2380 capillary column, Supelco, Sigma-Aldrich) analysis. Meanwhile, standard sugars were used to calibrate sugar content determined in samples. The insoluble precipitate from the AIR treated with TFA was collected for crystalline cellulose content determination. The Updegraff reagent (acetic acid:nitric acid:water, 8:1:2, v/v) was added to the precipitate and incubated at 100°C for 30 min. After washing with H₂O and acetone, the precipitate was incubated with 72% sulfuric acid at room temperature for 1 h. The content of crystalline cellulose was determined by anthrone assay (Foster et al., 2010b). For lignin measurement, AIR was incubated with freshly prepared acetyl bromide (25%, acetyl bromide in acetic acid) at 50°C for 3 h. After cooling, the AIR was mixed with 2 M NaOH, 0.5 M fresh hydroxylamine hydrochloride, and acetic acid. Lignin content was determined using a microplate reader (Varioskan Flash, Thermo) (Foster et al., 2010a).

Total RNA was extracted from wood-forming tissues using a mirVana miRNA Isolation Kit (Ambion-1561) following the manufacturer's instruction. After being treated by RNase-free DNase I (Sigma, 4716728001), the quality of total RNA was assessed on NanoDrop spectrophotometer (NanoDrop 2000, Thermo Scientific) and on agarose gel electrophoresis. For RNA sequencing (RNA-seq), cDNA library was generated from 5 μg of total RNA with TruSeq Stranded mRNA LTSample Prep Kit (Illumina, RS-122-2101) and Agencourt AMPure XP (BECKMAN COULTER, A63881). cDNA library was qualified through length distribution of fragments using Agilent 2100 (Bioanalyzer). The 150-bp paired-end sequencing was performed using platform of Illumina HiSeq X10. About five million reads per samples were generated.

Genomic DNA was extracted from wood-forming tissues using QIAamp DNA Mini kit (Cat.51306, Qiagen). DNA quantification and integrity were determined by a Nanodrop spectrophotometer (Thermo Fisher Scientific, Inc., Wilmington, DE) and 1% agarose electrophoresis, respectively. Before bisulfite treatment, lambda DNA was added to the purified DNA, which was used as an internal reference to calculate the conversion rate. The mixed DNA was then bisulfite treated using a Zymo Research EZ DNA methylaiton-Glod Kit (Zymo, D05005). Bisulfite sequencing (BS-seq) libraries were constructed by TruSeq® DNA Methylation Kit (Illumina, EGMK91396) following the manufacturer's instruction. After libraries were qualified, sequencing was performed on the Illumina HiSeq X Ten platform and 150 bp paired-end reads were generated.

Raw reads of sequencing were processed using NGS QC Toolkit to remove low-quality reads (Patel and Jain, 2012). The cleaned reads were mapped to Populus trichocarpa's genome (http://phytozome.jgi.doe.gov/) using hisat2 with default parameters (Kim et al., 2015). Gene expression level was measured as fragments per kilobase per million reads (FPKM) using cufflinks (Trapnell et al., 2010; Roberts et al., 2011). Read counts for each gene in each sample were obtained using htseq-count and standardized by rlog (Anders et al., 2015). Principle component analysis (PCA) was performed by plotPCA of DEseq2 R package with default parameters. Differential expression genes (DEGs) were identified using the DESeq R package by estimation of Size Factors and nbinomTest. Analysis of DEGs with gene ontology (GO) enrichment and the Kyoto Encyclopedia of Genes and Genomes (KEGG) (Kanehisa et al., 2008) pathway enrichment was performed using R based on the hypergeometric distribution.

The raw reads of BS-seq were cleaned using Fastp (Chen et al., 2018) by removing adapters, ploy-N, and low-quality reads. The remaining high-quality clean reads were mapped to the Populus trichocarpa's genome (http://phytozome.jgi.doe.gov/) using Bismark software with default parameters (Krueger and Andrews, 2011). Methylcytosine (mC) sites were identified using MethylKit (Akalin et al., 2012). With default parameters, MethylKit was applied for PCA analysis. Differentially methylated regions (DMRs) were identified using MethylKit software with a Q value (p-value corrected by FDR method) threshold of 0.05 and an absolute delta cutoff of 10% between the two groups. Analysis of DMGs with GO enrichment and KEGG pathway enrichment was performed according to the same method used for DEGs analysis.

The first-strand cDNA was synthesized from 2 μg of total RNA using a cDNA Synthesis SuperMix (TransGen Biotech, AT311-03). Using cDNA as template, quantitative real-time PCR (qRT-PCR) was performed using Perfectstart™ Green qPCR SuperMix (TransGen, AQ601) and a Quantstudio™ 3 Real-Time PCR Detection System (Thermo). The primers used for selected genes are listed in Supplementary Table 13, and TUB9 was used as an internal control to normalize gene expression.

To examine the properties of the JW and MW produced in Populus, plantation-grown trees that were propagated from a single clone were sampled. Three trees at 2 and 8 years old were collected with trunk at breast height, respectively (Supplementary Figure 1). In wood anatomical section, difference in the ratio of fiber cell/vessel, the length and size of fibers and vessels was observed between JW and MW (Figures 1A–D). MW contained higher ratio of fiber cell/vessel, lower density of vessel cell in wood section, and longer and larger fiber cell and vessel than those in JW (Figures 1E–J). Chemical analysis indicated that MW contained higher content of crystalline cellulose and lower content of lignin compared with JW (Table 1). Sugar composition in hemicelluloses also showed difference between JW and MW. MW contained higher xylose, mannose, glucose, and arabinose but lower galactose compared with JW (Table 1). These results indicated that JW and MW in Populus displayed different cellular structures and chemical compositions.

To dissect the gene expression involved in Populus wood formation, transcripts were profiled in the wood-forming tissues undergoing formation of JW and MW via high-throughput RNA-seq. Assessment of the RNA-seq data and biological repeats validated the high quality of the sequence data generated from the wood-forming tissues (Supplementary Table 1, The raw data in Sequence Read Archive (SRA), ID: PRJNA705066). A total of 300.5 million raw reads were obtained from six samples, and about 284.8 million high-quality reads (more than 94% of raw reads) were obtained after filtering and removal of low-quality reads. More than 86% of the high-quality reads per sample were mapped to the reference genome, corresponding to expression of ~20,000 genes out of 41,335 predicted genes in each sample. Meanwhile, PCA indicated that the transcript profiles showed a clear separation between JW and MW (Supplementary Figure 2A), suggesting the different transcription activities in formation of JW and MW.

A group of 3,992 genes were identified (FPKM ≥3, fold change >2, and p value FDR < 0.05) for their differential expression in JW and MW (Supplementary Figure 2B; Supplementary Table 2). The DEGs included 2,110 higher expression in JW and 1,882 higher expression in MW. As indicated by GO and KEGG enrichment analysis, the DEGs were primarily associated with plant hormones signaling and response, cell wall formation and modification, microtubule, cell organization and biogenesis, transcription, and other biological processes (Figure 2A; Supplementary Tables 3, 4).

Among the detected DEGs were included ample hormone-related genes. Particularly, genes of auxin transportation, brassinosteroids (BR) biosynthesis, and signaling were identified for their remarkable difference of expression in JW and MW (Figure 2B). The homologs of AUX1, LAX3, PILS2, and ABCB19 which are related to auxin transport (Enders and Strader, 2015) had high expression level in JW and MW (FPKM >100). Intriguingly, the homolog (Potri.016G113600) of AUX1 which facilitates auxin influx was expressed in JW higher than in MW, while the homolog (Potri.017G081100) of ABCB19 which facilitates the efflux of auxin was expressed in MW higher than in JW. Furthermore, four homologs of PIN1, PIN3, and PIN6, encoded the auxin transporters that mediate that auxin efflux (Liu et al., 2014; Enders and Strader, 2015), were identified in DEGs, and their expression in MW was much higher than in JW (Figure 2B). In addition, the PIN-LIKES (PILS), which are thought to be located on the endoplasmic reticulum (ER), may transport auxin from the cytoplasm into the ER (Enders and Strader, 2015). One PILS2 homolog (Potri.004G093200) was expressed in MW 13 times higher than in JW (Figure 2B). The differential expression of the auxin-related genes was verified by qRT-PCR determination (Figure 2C). These data indicate that the genes involved in IAA transport were differentially expressed in MW and JW.

Genes involved in BR biosynthesis and signaling were readily noticed among DEGs. The homolog (Potri.008G084800) of DIM/DWF1 which is a key gene for BR biosynthesis in Arabidopsis (Klahre et al., 1998; Youn et al., 2018) was expressed in MW 5.6 times higher than in JW (Figure 2B). Additionally, several genes involved in BR signaling were differentially expressed between JW and MW (Figure 2B). The homolog of BRI1-associated receptor kinase (BAK1) (Potri.001G206700), brassinosteroid signaling positive regulator (BZR1) (Potri.005G126400), and BR signaling kinase 1 (BSK1) (Potri.002G011800) were downregulated in MW compared with JW. qRT-PCR determination confirmed the differential expression (Figure 2C). On the other hand, although the genes involved in other hormones signaling such as gibberellins (GAs) were detected to be differentially expressed between JW and MW (Figure 2B), their expression profiles were unable to be verified. Together, the results suggest that auxin and BR may be involved in the regulation of JW and MW formation.

The identified 3,992 DEGs included 305 transcription factor (TF) genes, 7.6% of all DEGs, which is higher than 6% of TF genes in Populus genome (PlantTFDB, http://planttfdb.cbi.pku.edu.cn) (Jin et al., 2014) (Figure 3A; Supplementary Table 5), implying that expression of transcription factor genes is altered in higher proportion. Among the TF genes, 105 TF genes were upregulated in MW, and the rest were downregulated in MW. GO annotation analysis indicated that the TFs in the DEGs were primarily associated with hormone biosynthesis and responses, plant growth and development, cell fate, cell wall formation, and response to abiotic and biotic stimuli (Figure 3B).

Several TF genes that are involved in auxin signaling were differentially expressed in JW and MW (Supplementary Table 5). The DEGs contained nine auxin response factor (ARF) genes with four upregulated and five downregulated in JW. In addition, three homologs of SHI-RELATED SEQUENCEs (SRSs) that play a role in activation of auxin biosynthesis (Eklund et al., 2010) were upregulated in MW. One homolog of AINTEGUMENTA-LIKE 6 (AIL6) that is involved in regulation of auxin biosynthesis (Pinon et al., 2013) displayed higher expression in MW. The differential expression determined by RNA-seq and qRT-PCR analyses was well-correlated (Figure 3C; Supplementary Figure 3).

The DEGs included a number of TFs related to regulation of cell wall formation, such as the homologs of SND1, NST1 VND1, and VND4 (Supplementary Table 5), which were key TFs for regulation of secondary cell wall biosynthesis (Zhong et al., 2010; Hussey et al., 2013; Kumar et al., 2016). PtrMYB26, PtrMYB90, and PtrMYB152 which were reported to directly regulate lignin biosynthesis in Populus (Zhong et al., 2011; Wang et al., 2014; Li et al., 2015) were expressed higher in JW than in MW. This is in agreement with the higher lignin content in JW. Furthermore, expression of PtrMYB170 and PtrMYB121, which may play a role in resource acquisition and allocation for xylem development (Romano et al., 2012), were upregulated in MW. One homolog (Potri.006G241700) of MYB3R1, which played a role in activating expression of the genes in cell cycle (Haga et al., 2011), showed higher expression in MW. PtrERF118 (Potri.018G028000), which was identified as a TF-regulating xylem cell expansion in Populus (Vahala et al., 2013; Seyfferth et al., 2018), was upregulated in MW. These results suggest that the transcriptional networks in relation to auxin biosynthesis and signaling, secondary cell wall formation, and xylem cell differentiation were differentially regulated in the formation of JW and MW.

The DEGs included a large number of genes responsible for cell wall biosynthesis (Figure 4; Supplementary Table 6). It is worth noting that among the DEGs, a large number of genes related to turgor maintenance and cell expansion were identified, of which the majority were upregulated in MW (Figure 4A; Supplementary Table 6). Cell turgor pressure is closely related to cell expansion, it can induce irreversible cell expansion (Genard et al., 2001). Five TIPs and three PIPs which are related to maintenance of turgor pressure were upregulated in MW. On the other hand, nine EXPAs and four XTHs which were related to cell wall loosening (Mcqueenmason et al., 1992; Van Sandt et al., 2007; Nishikubo et al., 2011) were identified upregulated in MW (Figure 4A). In addition, the homologs of FUC1, PAEs, PMEs, and PMRs, which play a role in modifying cell wall for cell wall loosening (Gou et al., 2012; Kato et al., 2018), were also upregulated in MW. Meanwhile, the homolog of HA11 (a PM H+-ATPase) were upregulated in MW. The PM H+-ATPase can reduce the pH in the apoplast space to activate expansins and other cell wall loosening proteins as well as promote the absorption of water to provide turgor pressure for cell expansion (Spartz et al., 2014). These changes of the gene expression were confirmed by qRT-PCR (Figure 4B). These results support the observation that MW is formed with larger and longer fibers and vessel elements.

Expression of the genes for monolignol biosynthesis was consistently lower in MW, including five PAL homologs, C4H1 and C4H2, 4CL1 and 4CL5, HCT1, C3H3, CSE2, CCoAOMT1 and CCoAOMT2, CCR2, F5H1 and F5H2, COMT1, and CAD (Figure 4A). This may reflect less lignin biosynthesis in MW. After biosynthesis, monolignols are polymerized by laccases or peroxidases. Among detected, 11 laccase (LAC) genes including homologs of LAC2, LAC10, LAC11, and LAC17 were downregulated in MW. Interestingly, among the 10 detected peroxidase genes, two peroxidases were expressed higher in JW and the other eight were expressed higher in MW. It is worthy of further investigating whether the different members of LAC or PRX act in different phases of wood formation. Cellulose is synthesized by cellulose synthase complex (CSC) of synthases (CesAs) (Song et al., 2010; McFarlane et al., 2014). It is interesting to notice that the homologs of CesA4, CesA7, and CesA8, which form CSCs for cell wall thickening (Song et al., 2010; Watanabe et al., 2015; Xi et al., 2017), were downregulated in MW (Figures 4A,B). As cellulose synthesis is affected by CesA modifications at protein level (Polko and Kieber, 2019), it is unclear whether modification of CesAs are involved in regulating cellulose synthesis in JW and MW.

A number of genes for biosynthesis of hemicelluloses were differentially expressed in JW and MW (Figure 4A). UDP-glucuronic acid decarboxylase (UXS) catalyzes UDP-glucuronic acid (UDP-GlcA) to biosynthesis of UDP-Xyl (Kuang et al., 2016), which is a donor for biosynthesis of xylan, a major secondary cell wall hemicellulose. Cellulose synthase-like D (CSLD) involved xylan synthesis (Bernal et al., 2007). Cellulose synthase-like A (CSLA) encodes mannan synthase (Liepman et al., 2005; Suzuki et al., 2006; Verhertbruggen et al., 2011). Homologs of UXS, CSLC, CSLD, and CLSA showed higher expression in MW. The results support a higher content of xylan and mannan deposited in MW than in JW (Table 1).

The differential gene expression in different growth phases prompted us to examine the whole genome bisulfite sequencing (WGBS). The bisulfite sequencing showed that 87.6–91.9% of the reads was qualified for methylation assay against the Populus genome (http://phytozome.jgi.doe.gov/) (Supplementary Table 7, The raw data in Sequence Read Archive (SRA), ID: PRJNA705570). Overall, the methylation level was different within cytosine methylation contexts (CG, CHG, and CHH). The context of CG had higher methylation level, while CHG and CHH had lower (Supplementary Table 8). The DNA methylation context patterns displayed a similarity with those previously observed in Populus (Vining et al., 2012; Su et al., 2018). PCA showed that JW and MW had distinct DNA methylation (Figure 5A). Comparison of the DNA methylation in JW and MW revealed 12,176 differentially methylated regions (DMRs) (with methylation difference ≥10, Q-value < 0.05). Majority of DMRs were in the contexts of CG sites (10,303) and CHG sites (1,663) (Supplementary Figure 4A; Supplementary Table 9), Among them, 10,237 DMRs were located in gene body and/or flanking regions (±2 kb), named differentially methylated genes (DMGs). In MW, 5,414 DMGs showed higher methylation while JW contained 4,849 DMGs with higher methylation (Figure 5B; Supplementary Table 10), suggesting that different DNA methylations occurred in the formation of JW and MW. Analysis of the correlation between DMGs and DEGs indicated that DMRs in gene promoter region were more likely to affect gene expression (Supplementary Figure 4B). About 20% DEGs (802) displayed different methylation (Supplementary Table 11). These DEGs were closely related to plant hormone signaling and response, cell wall formation and modification, metabolic process, transcription and translation, etc. (Supplementary Figure 4C; Supplementary Table 12). For example, the homologs of ARFs, BAK1, BSK1, and BZR1, which are involved in auxin and BR signaling, were differential methylated in their different gene regions in JW and MW (Figure 5C; Supplementary Table 11). Furthermore, several genes related to cell wall formation such as XTH30, PAEs, WND1B, CESA4, CESA7, and CESA8 (Figure 5C; Supplementary Table 11) showed differential methylation in JW and MW. In addition, several DMRs in intergenic region were neighbored to the homologs of PILS2, AUX1, PIN7, WND2A, MYBs, and PAL1, which are involved in auxin distribution and cell wall biosynthesis (Supplementary Table 11). In summary, the results revealed that DNA methylation displayed a clear difference in the formation of JW and MW, which may play a role in regulating gene expression in different growth phases, particularly for the genes involved in hormone signal transduction, cell division, and cell wall biosynthesis in wood formation.