Transcriptomic Analysis of Seasonal Gene Expression and Regulation during Xylem Development in “Shanxin” Hybrid Poplar ( Populus davidiana × Populus bolleana )

: Xylem development is a key process for wood formation in woody plants. To study the molecular regulatory mechanisms related to xylem development in hybrid poplar P. davidiana × P. bolleana , transcriptome analyses were conducted on developing xylem at six different growth stages within a single growing season. Xylem development and differentially expressed genes in the six time points were selected for a regulatory analysis. Xylem development was observed in stem sections at different growth stages, which showed that xylem development extended from the middle of April to early August and included cell expansion and secondary cell wall biosynthesis. An RNA-seq analysis of six samples with three replicates was performed. After transcriptome assembly and annotation, the differentially expressed genes (DEGs) were identiﬁed, and a Gene Ontology (GO) enrichment analysis, Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis and expression analysis of the DEGs were performed on each sample. On average, we obtained >20 million clean reads per sample, which were assembled into 84,733 nonredundant transcripts, of which there were 17,603 unigenes with lengths >1 kb. There were 14,890 genes that were differentially expressed among the six stages. The upregulated DEGs were enriched in GO terms related to cell wall biosynthesis between S1 vs. S2 or S3 vs. S4 and, in GO terms, related to phytohormones in the S1 vs. S2 or S4 vs. S5 comparisons. The downregulated DEGs were enriched in GO terms related to cell wall biosynthesis between S4 vs. S5 or S5 vs. S6 and, in GO terms, related to hormones between S1 vs. S2 or S2 vs. S3. The KEGG pathways in the DEGs related to “phenylpropanoid biosynthesis”, “plant hormone signal transduction” and “starch and sucrose metabolism” were signiﬁcantly enriched among the different stages. The DEGs related to cell expansion, polysaccharide metabolism and synthesis, lignin synthesis, transcription factors and hormones were identiﬁed. The identiﬁcation of genes involved in the regulation of xylem development will increase our understanding of the molecular regulation of wood formation in trees and, also, offers potential targets for genetic manipulation to improve the properties of wood.


Introduction
Due to the short growth period, rapid growth, strong adaptability and strong shading ability, species of Populus are widely planted all over the world. Poplars have a small genome and strong asexual reproduction ability, making it a model plant for genetic improvement of forest tree breeding. P. davidiana × P. bolleana is an interspecific hybrid of two species in Populus sect Populus. P. davidiana, the female parent, was collected from the Gaofeng Forest Farm in Nenjiang County, and the male parent Xinjiang poplar (P. bolleana) was from Urumqi [1]. The hybrid was named "Shanxin" poplar. After 20 years of field testing, the hybrid's performance was considered to be good, and the traits were found to be stable. "Shanxin" poplar has the characteristics of rapid growth, a straight and full trunk, white and delicate material, smooth green bark, graceful posture, a narrow crown and resistance to cold and drought. It is an excellent fast-growing poplar species for afforestation and urban and rural greening in arid and cold regions in the Northern hemisphere [2]. Since "Shanxin" poplar cuttings do not root easily, the survival rate is low. Therefore, research on this hybrid mainly focuses on asexual reproduction, such as cutting propagation and tissue culture and rapid propagation methods. It flowers but does not set seed and is an ideal recipient for genetic transformation to obtain safe genetically modified varieties. The use of transgenic technology can improve the ability of cuttings to root, as well as resistance to salt and alkali and resistance to pests and diseases in a short time. At present, although it has good wood properties, there are few studies focusing on wood formation in "Shanxin" poplar. Understanding the molecular biology of wood formation is important for the potential molecular improvement of "Shanxin" poplar and other species of woody plants.
Broadly speaking, the secondary xylem is the wood. The growth and development of plants and their successful adaptation to various environments depend to a large extent on the conduction of water, nutrients and other important element [3]. Xylem is a special vascular tissue that acts as a pipeline for water and mineral transport that also provides mechanical support for vertical growth. In woody plants, wood is mainly composed of secondary xylem elements [4]. Over the past two decades, several key factors, including hormones, signal transduction, transcriptional regulation and post-transcriptional regulators, have been shown to be the main factors controlling the formation of xylem [5]. The transcript profiling of wood developmental stages during this developmental sequence has revealed the intricate transcriptional regulation of wood formation [6,7]. The formation of xylem is a complex developmental process that requires cell specification and secondary cell wall deposition that is controlled at the transcriptional level. It is often affected by external factors, such as photoperiod, nutrient availability, water content and temperature, and by internal stimuli such as plant hormones and photosynthetic assimilation and the interactions between these factors [8]. Thousands of genes are actively involved in the formation of xylem and have been identified in both herbaceous and woody plants. This will expand our understanding of the regulation, metabolism and biosynthesis processes of xylem formation.
The development of xylem is also subject to seasonal regulation [9,10]. Studies of the cytological changes that occur in Chinese fir cambium have shown that the number of vascular cambium cells is positively correlated with the width of the cambium during the entire seasonal cycle. The number of cambium cells and the radial diameter of cambium cells are closely related to xylem formation [11]. The increase in the perimeter of trees is caused by the activity of the vascular cambium, which produces secondary xylem and phloem [12,13]. In the Northern temperate regions, the vascular cambium of woody plants has a seasonal activity and sleep cycle. The initiation, duration and cessation of vascular cambium activity are usually decisive factors in determining the quantity and quality of wood [11,14].
In our study, we used "Shanxin" poplar as the research subject to conduct anatomical observations of xylem development at different developmental stages during the growing season. Our objectives were to study the expression profile of genes involved in molecular regulation of the seasonal development of xylem, to increase our understanding of the molecular regulatory network that controls the formation of wood properties in "Shanxin" poplar and to provide data and materials for the molecular breeding and improvement of poplar.

Experimental Design
In order to analyze the seasonal expression and regulation of genes related to xylem development in "Shanxin" poplar, we selected six time points during xylem development to analyze the changing seasonal patterns. The samples were taken on 22 April (S1), 12 May (S2), 3 June (S3), 25 June (S4), 15 July (S5) and 4 August (S6). The "Shanxin" poplar is a wild type outside. All samples are sent to Biomarker (Beijing, China) for RNA extraction at the same time.

Plant Materials
The 5-to 6-year-old test trees were selected from the breeding base of Improved Variety Bases in Fulaerji District, Qiqihar City, Heilongjiang, China, 47 • 15 N, 123 • 45 E. The average maximum temperature at 6 time points was almost 15 • C (S1), 19 • C (S2), 25 • C (S3), 28 • C (S4), 30 • C (S5) and 25 • C (S6). At S1 and S8, the bark was difficult to separate from the xylem. Three well-developed trees with consistent growth were used at each time point as a group, with three biological replicates. Samples were taken every three weeks. A 20-cm stem section was cut at a height of 130 cm from the ground, and the developing xylem was scraped out using a surgical blade. After flash freezing in liquid nitrogen, the samples were stored at −80 • C prior to RNA extraction. At the same time, small blocks of developing xylem were excised from the bottom of each stem and fixed in formalin-alcohol-acetic acid for anatomical observation.

Anatomical Observation of Developing Xylem
To make slides for microscopic observation, stem sections~3 cm were cut and were fixed in FAA (50% ethanol:formaldehyde:acetic acid = 90:5:5). The samples were evacuated for 10 min using a vacuum pump, and the valve was then opened, and the samples were allowed to stand for 10 min. This was repeated several times until there were no more bubbles produced under vacuum, and the material sank to the bottom of the solution. After 1 to 2 weeks at room temperature, the fixed samples were cut into 40-mm slices with a slide-away microtome (Leica 1007, Wetzlar, Germany) and then stained with Safranin and Fast Green for observation under a stereo light microscope (Olympus SZX7, Tokyo, Japan).

Library Construction and Transcriptome Sequencing
Total RNA was extracted using CTAB method according to the paper of Liu et al. [15]. The quality and concentration of total RNA was evaluated with a NanoDrop 2000 UVvisible Spectrophotometer (Thermo Scientific, Waltham, MA, USA). Equal amounts of total RNA from the treated samples were used for transcriptome sequencing.
The purity, concentration and integrity of RNA samples were tested using agarose gel electrophoresis to ensure the use of qualified samples for transcriptome sequencing.
A total of 1 µg of total RNA per sample was used as input material for RNA-seq library preparation. Sequencing libraries were generated using the NEBNext ® Ultra™ RNA Library Prep Kit for Illumina ® (New England Biolabs, Ipswich, MA, USA) following the manufacturer's recommendations, and index codes were added to attribute sequences to each sample. Briefly, mRNA was purified from total RNA using poly-T oligo-attached magnetic beads. Fragmentation was carried out using divalent cations under elevated temperature in NEBNext First Strand Synthesis Reaction Buffer (5X). First-strand cDNA was synthesized by using random hexamer primers and M-MuLV Reverse Transcriptase. Second-strand cDNA synthesis was subsequently performed using DNA Polymerase I and RNase H. The remaining terminal overhangs were converted into blunt ends via exonuclease/polymerase activities. After adenylation of the 3 ends of the DNA fragments, the NEBNext Adaptor with a hairpin loop structure was ligated to prepare for hybridization. In order to select cDNA fragments of~240 bp in length, the library fragments were purified using the AMPure XP system (Beckman Coulter, Beverly, MA, USA). Uracil residues were removed with 3 µL of the USER (Uracil-Specific Excision Reagent) Enzyme (NEB, Ipswich, MA, USA) on the size-selected, adaptor-ligated cDNA at 37 • C for 15 min, followed by 5 min at 95 • C prior to PCR. Amplification was then performed using Phusion High-Fidelity DNA polymerase, Universal PCR primers and the Index (X) Primer. The PCR products were bead-purified (AMPure XP system, Beckman Coulter, Beverly, MA, USA), and the library quality was assessed on the Agilent Bioanalyzer 2100 system.
Clustering of the index-coded samples was performed on a cBot Cluster Generation System using a TruSeq PE Cluster Kit v3-cBot-HS (Illumina), according to the manufacturer's instructions. After cluster generation, the library preparations were sequenced on an Illumina Hiseq 2000 platform, and paired-end reads were generated.

Transcriptome Assembly and Annotation
Transcriptome assembly was accomplished based on the left.fq and right.fq using Trinity [16] with min_kmer_cov set to 2 by default, and the default settings were used for all other parameters. Gene function was annotated based on the following databases: NR (NCBI nonredundant protein sequences), Pfam (Protein family), KOG/COG/eggNOG (Clusters of Orthologous Groups of proteins), Swiss-Prot (a manually annotated and reviewed protein sequence database), KEGG (Kyoto Encyclopedia of Genes and Genomes) and GO (Gene Ontology).

Identification of DEGs and GO Enrichment Analysis
Differential expression analysis of the two conditions/groups was performed using the DESeq R package (1.10.1). DESeq provides statistical routines for determining differential expression in digital gene expression data using a model based on the negative binomial distribution [17] for controlling the false discovery rate. Genes with an adjusted p-value < 0.05 found by DESeq were determined to be differentially expressed. Gene Ontology (GO) enrichment analysis of the DEGs was implemented by the topGO R package based on the Kolmogorov-Smirnov test.

KEGG Pathway Enrichment Analysis
The Kyoto Encyclopedia of Genes and Genomes [18] is a database resource for understanding high-level functions and utilities in biological systems, such as the cell, the organism and the ecosystem, from molecular-level information, especially large-scale molecular datasets generated by genome sequencing and other high-throughput experimental technologies. We used KOBAS [19] software to test the statistical enrichment of DEGs in KEGG pathways.

Histological Observations of Xylem Development
We performed observations on the stem sections of "Shanxin" poplar at six different time points. Figure 1 shows the sections of "Shanxin" poplar stained with Safranin O-Fast green at the different stages. At S1, the xylem has not yet developed. At S2, the developing xylem is dyed green, which suggests that the xylem has begun to develop, and the xylem cells do not yet have lignified secondary walls, suggesting that they are still making their primary wall and/or expanding. At S3, the developing xylem cells are dividing rapidly, and the xylem is expanding, and most of the cells begin to synthesize secondary cell walls. These processes continue until S5, when the rate of xylem expansion is slightly reduced. There were no obvious microscopic differences between S5 and S6, and the formation of latewood was observed ( Figure 1).

RNA Sequencing and de Novo Transcriptome Assembly
The cDNA library was constructed using the high-quality RNA extracted from the xylem of "Shanxin" poplar at the six different times. In order to characterize the "Shanxin" poplar transcriptome, the cDNA library was sequenced using the Illumina HiSeq2000 platform. After removing the low-quality reads, adaptors and ambiguous sequences, we generated a total of 121.93 Gb of clean sequencing reads with an average GC content of 44.86%. We then used Trinity [16] to assemble these high-quality reads resulting in 304,694 contigs with an average length of 1633.28 bp and an N50 length of 1764 bp (Table 1). We further assembled 84,733 unigenes with an average length of 804.05 bp. Among these unigenes, 51,512 (60.79%) were >300 bp in length, 31,892 (37.63 %) were >500 bp, 17,603 (20.77 %) were >1000 bp and 9535 (11.25%) were >2000 bp in length ( Figure 2).

RNA Sequencing and De Novo Transcriptome Assembly
The cDNA library was constructed using the high-quality RNA extracted from the xylem of "Shanxin" poplar at the six different times. In order to characterize the "Shanxin" poplar transcriptome, the cDNA library was sequenced using the Illumina HiSeq2000 platform. After removing the low-quality reads, adaptors and ambiguous sequences, we generated a total of 121.93 Gb of clean sequencing reads with an average GC content of 44.86%. We then used Trinity [16] to assemble these high-quality reads resulting in 304,694 contigs with an average length of 1633.28 bp and an N50 length of 1764 bp (Table 1). We further assembled 84,733 unigenes with an average length of 804.05 bp. Among these unigenes, 51,512 (60.79%) were >300 bp in length, 31,892 (37.63 %) were >500 bp, 17,603 (20.77 %) were >1000 bp and 9535 (11.25%) were >2000 bp in length ( Figure 2).

Functional Annotation of the Unigenes
By using BLAST [20]with an E-value parameter of no greater than 1e −5 and HMMER [21] with an E-value no greater than 1e −10 , the assembled unigene sequences were used as search queries against the public databases, including the NCBI nonredundant protein database (Nr) [22]; the Clusters of Orthologous Groups (COGs) [23] database; the euKaryotic Orthologous Groups (KOG) [24] database; the Swiss-Prot [25] protein database; the Gene Ontology (GO) [26] database; the evolutionary genealogy of genes: Non-supervised Orthologous Groups (eggNOG) [27] database, the Protein family (Pfam) [28] database and the Kyoto Encyclopedia of Genes and Genomes (KEGG) [29] database. Altogether, 43,525 of the unigenes were matched to a sequence in at least one of the above-mentioned databases. There were also matches to 40,513 and 37,441 unigenes in the Nr and eggNOG databases, respectively. A small number of unigenes were also used to search another six databases, as shown in Table 2. The GO analysis is widely used to assign recognized gene function to uncharacterized sequences. Based on information from the Nr annotation and the WEGO program, we assigned 40,513 unigenes to 50 different functional groups (Figure 2), and these groups were separated into the three major GO terms "biological process", "cellular component" and "molecular function" using BLAST2GO. In the "biological process" category, the unigenes related to "metabolic process" and "cellular process" accounted for a large proportion, followed by "single-organism process" and "biological regulation". In the "cellular component" category, the GO terms "membrane part", "organelle part", "cell part" and

Functional Annotation of the Unigenes
By using BLAST [20] with an E-value parameter of no greater than 1 × 10 −5 and HMMER [21] with an E-value no greater than 1 × 10 −10 , the assembled unigene sequences were used as search queries against the public databases, including the NCBI nonredundant protein database (Nr) [22]; the Clusters of Orthologous Groups (COGs) [23] database; the euKaryotic Orthologous Groups (KOG) [24] database; the Swiss-Prot [25] protein database; the Gene Ontology (GO) [26] database; the evolutionary genealogy of genes: Non-supervised Orthologous Groups (eggNOG) [27] database, the Protein family (Pfam) [28] database and the Kyoto Encyclopedia of Genes and Genomes (KEGG) [29] database. Altogether, 43,525 of the unigenes were matched to a sequence in at least one of the above-mentioned databases. There were also matches to 40,513 and 37,441 unigenes in the Nr and eggNOG databases, respectively. A small number of unigenes were also used to search another six databases, as shown in Table 2. The GO analysis is widely used to assign recognized gene function to uncharacterized sequences. Based on information from the Nr annotation and the WEGO program, we assigned 40,513 unigenes to 50 different functional groups (Figure 2), and these groups were separated into the three major GO terms "biological process", "cellular component" and "molecular function" using BLAST2GO. In the "biological process" category, the unigenes related to "metabolic process" and "cellular process" accounted for a large proportion, followed by "single-organism process" and "biological regulation". In the "cellular component" category, the GO terms "membrane part", "organelle part", "cell part" and "cell" predominated, while "catalytic activity" and "binding" were the two most represented GO terms in the "molecular function" category.
We used the Cluster of Orthologous Groups (COG) database, which was constructed based on the systematic evolutionary relationships of bacteria, algae and eukaryotes to perform a phylogenetic classification of the "Shanxin" poplar gene sequences (Supplementary Table S1). For example, in the S1 vs. S2 comparison, "Carbohydrate transport and metabolism" represented the largest group, followed by "Signal transduction mechanism". In addition, only a few unigenes grouped in the "Energy production and conversion" and "Intracellular trafficking, secretion, and vesicular transport" categories. Compared with the eggNOG database, which contains the functional description and functional classification of orthologous proteins, we can see the differences between different time points more clearly (Figures 3 and 4). "cell" predominated, while "catalytic activity" and "binding" were the two most represented GO terms in the "molecular function" category. We used the Cluster of Orthologous Groups (COG) database, which was constructed based on the systematic evolutionary relationships of bacteria, algae and eukaryotes to perform a phylogenetic classification of the "Shanxin" poplar gene sequences (Supplementary Table S1). For example, in the S1 vs. S2 comparison, "Carbohydrate transport and metabolism" represented the largest group, followed by "Signal transduction mechanism". In addition, only a few unigenes grouped in the "Energy production and conversion" and "Intracellular trafficking, secretion, and vesicular transport" categories. Compared with the eggNOG database, which contains the functional description and functional classification of orthologous proteins, we can see the differences between different time points more clearly (Figures 3 and 4) Figure 3. Cluster of Orthologous Groups (COG) functional annotation classification of the differentially expressed genes (DEGs) isolated in S1 vs. S2 from the "Shanxin" poplar xylem.  "cell" predominated, while "catalytic activity" and "binding" were the two most represented GO terms in the "molecular function" category. We used the Cluster of Orthologous Groups (COG) database, which was constructed based on the systematic evolutionary relationships of bacteria, algae and eukaryotes to perform a phylogenetic classification of the "Shanxin" poplar gene sequences (Supplementary Table S1). For example, in the S1 vs. S2 comparison, "Carbohydrate transport and metabolism" represented the largest group, followed by "Signal transduction mechanism". In addition, only a few unigenes grouped in the "Energy production and conversion" and "Intracellular trafficking, secretion, and vesicular transport" categories. Compared with the eggNOG database, which contains the functional description and functional classification of orthologous proteins, we can see the differences between different time points more clearly (Figures 3 and 4) Figure 3. Cluster of Orthologous Groups (COG) functional annotation classification of the differentially expressed genes (DEGs) isolated in S1 vs. S2 from the "Shanxin" poplar xylem.

Identification of Differentially Expressed Genes (DEGs)
In order to fully understand the seasonal response of xylem development in the "Shanxin" poplar and to identify related genes, we constructed six cDNA libraries from samples taken at six time points from 22 April to 4 August (samples S1-S6). We then used DESeq2 to analyze the differential gene expression between the sample groups to obtain the DEG sets. Based on the standard criteria used to identify DEGs, including a low false discovery rate (FDR, <10 −3 ) and |log 2 FPKM| > 1.0 [30], a difference in the number of DEGs in comparison of two samples was obtained. Taking S1 as the control, a total of 4273 unigenes were considered to be differentially expressed between S1 and S2 (2091 upregulated and 2182 downregulated), 8282 unigenes between S1 and S3 (4014 upregulated and 4268 downregulated), 8587 unigenes between S1 and S4 (4266 upregulated and 4321 downregulated), 6852 unigenes between S1 and S5 (3833 upregulated and 3019 downregulated) and 5875 unigenes between S1 and S6 (3308 upregulated and 2567 downregulated). This shows clearly that, as the seasons change, the trend in the number of differential genes shows an initial increase followed by a decrease. The greatest number of DEGs (8587) was enriched in the S1 vs. S4 group, while the fewest (907) were enriched in the S3 vs. S4 group (Table 3). For the S4 vs. S5 comparison, there were more upregulated DEGs than downregulated DEGs. There were more downregulated DEGs than upregulated DEGs in the S5 vs. S6 comparison, whereas the other comparisons gave a roughly similar number of upregulated and downregulated DEGs.

GO Enrichment Analysis of the DEGs
To explore the biological processes involved in xylem development during the seasonal change, pairwise comparisons of the samples were performed, including S1 vs. S2, S2 vs. S3, S3 vs. S4, S4 vs. S5 and S5 vs. S6 (Supplementary Table S2). Significantly enriched GO terms were then determined using a hypergeometric test to identify the functional clusters and biochemical pathways (corrected p-value). GO terms related to the different processes of xylem development were identified, such as terms related to hormones, cell wall polysaccharide biosynthesis and metabolism, lignin biosynthesis and metabolism and cell wall or xylem. Most of the GO terms related to cell wall formation and thickness in the upregulated genes were enriched in S1 vs. S2 and S2 vs. S3, while the GO terms related to hormones were enriched in S1 vs. S2 and S4 vs. S5 (p < 0.05; Figure 5). This suggests that these biological processes may play important roles in the corresponding stages of xylem development. In the downregulated genes, GO terms related to cell wall formation and thickness were enriched in S4 vs. S5 and S5 vs. S6, while GO terms related to hormones were enriched in S1 vs. S2 and S2 vs. S3 (p < 0.05; Figure 6).

Pathway Enrichment Analysis of the DEGs
Using the KEGG database, the pathways with significant changes (Q value ≤ 0.01) in xylem development among the different seasonal stages were determined ( Table 4). The analysis showed that the "phenylpropanoid biosynthesis" pathway was significantly enriched in the DEGs in the S1 vs. S2, S1 vs. S5 and S1 vs. S6 comparisons; "plant hormone signal transduction" was enriched in the DEGs in the S1 vs. S2 and S1 vs. S3 comparisons; "starch and sucrose metabolism" was enriched in DEGs in all comparisons and "amino sugar and nucleotide sugar metabolism" was enriched in the DEGs in the S1 vs. S3 and S1 vs. S4 comparisons. In addition, "endocytosis", "plant-pathogen interaction", "phagosome" and "diterpenoid biosynthesis" were also significantly enriched in some comparisons. Some DEGs were found to be enriched in multiple KEGG pathways, such as c27167.graph_c0 (peroxidase 12-like (Populus euphratica)), c31851.graph_c0 (pathogenesis-related family protein (Populus trichocarpa)) and c20733.graph_c0 (probable fructokinase-1 (Populus euphratica)), which also shows the importance of these genes in xylem development of Shanxin Poplar.

Expression of Genes Related to Cell Expansion and Cell Wall Biosynthesis
Xylem development occurs via complex processes in vascular cambium cells through cell division, cell differentiation and expansion, secondary wall formation and programmed cell death [31]. Expansin and XTH genes are the two most important genes involved in cell expansion and cell wall synthesis [32]. We identified 10 differentially expressed expansin genes in the six different time points, and six of these genes were found to be highly expressed in stages S2 and S3 and downregulated in S5 and S6. Twelve XTH genes were identified, eight of which showed a significant upregulated expression after S2, with a declining expression at S5 (Figure 7). According to the comprehensive analysis of the gene expression and the cross-sectional section of the stem, it is speculated that these EXP and XTH genes may be involved in the regulation of cell expansion and cell wall synthesis in the early stages of xylem development.

Expression of Genes Related to Polysaccharide Biosynthesis and Metabolism
The cell wall includes the primary cell wall and secondary cell wall. Unlike the primary cell walls, which are initiated during cell division and develop along with the expansion of the cells, secondary cell walls are constructed after the cells have stopped growing [33].
The biosynthesis of cellulose in plant cell walls is mainly catalyzed by cellulose synthase (CesA) [34]. In the sequencing results, we identified 21 genes encoding cellulose synthase (CesA), 19 of which were upregulated in S2 and S3, including four that were significantly upregulated ( Figure 8). In S4 and S5, the expression levels were flat, while in S6, the expression decreased slightly. At the same time, we identified six sucrose synthase-related genes, four of which were significantly upregulated in S2 and S3, and their expression levels declined continuously in S4-S6 (Figure 9).

Expression of Genes Related to Cell Expansion and Cell Wall Biosynthesis
Xylem development occurs via complex processes in vascular cambium cells through cell division, cell differentiation and expansion, secondary wall formation and programmed cell death [31]. Expansin and XTH genes are the two most important genes involved in cell expansion and cell wall synthesis [32]. We identified 10 differentially expressed expansin genes in the six different time points, and six of these genes were found to be highly expressed in stages S2 and S3 and downregulated in S5 and S6. Twelve XTH genes were identified, eight of which showed a significant upregulated expression after S2, with a declining expression at S5 (Figure 7). According to the comprehensive analysis of the gene expression and the cross-sectional section of the stem, it is speculated that these EXP and XTH genes may be involved in the regulation of cell expansion and cell wall synthesis in the early stages of xylem development.

Expression of Genes Related to Polysaccharide Biosynthesis and Metabolism
The cell wall includes the primary cell wall and secondary cell wall. Unlike the primary cell walls, which are initiated during cell division and develop along with the expansion of the cells, secondary cell walls are constructed after the cells have stopped growing [33]. The biosynthesis of cellulose in plant cell walls is mainly catalyzed by cellulose synthase (CesA) [34]. In the sequencing results, we identified 21 genes encoding cellulose synthase (CesA), 19 of which were upregulated in S2 and S3, including four that were significantly upregulated (Figure 8). In S4 and S5, the expression levels were flat, while in S6, the expression decreased slightly. At the same time, we identified six sucrose synthaserelated genes, four of which were significantly upregulated in S2 and S3, and their expression levels declined continuously in S4-S6 ( Figure 9).    The sugar metabolism, including pentose phosphate and glycolysis, is ubiquitous in most organisms and plays an important role in the development of xylem in plants [34]. As shown in Supplementary Table S3, there are a total of 66 DEGs in S1 vs. S2, 149 in S1 vs. S3, 149 in S1 vs. S4, 110 in S1 vs. S5 and 107 in S1 vs. S6 related to polysaccharide biosynthesis and metabolism. It can be seen that the number of DEGs is the highest in S3, The sugar metabolism, including pentose phosphate and glycolysis, is ubiquitous in most organisms and plays an important role in the development of xylem in plants [34]. As shown in Supplementary Table S3, there are a total of 66 DEGs in S1 vs. S2, 149 in S1 vs. S3, 149 in S1 vs. S4, 110 in S1 vs. S5 and 107 in S1 vs. S6 related to polysaccharide biosynthesis and metabolism. It can be seen that the number of DEGs is the highest in S3, in which 73 are upregulated and 76 are downregulated. Among these genes, there are more DEGs involved in regulating the "starch and sucrose metabolism" pathways than other pathways. This finding indicates that the "starch and sucrose metabolism" pathway in the xylem of "Shanxin" poplar is more important during seasonal changes and that June is an important period in xylem development.
Changes in energy metabolism alter SCW accumulation and composition. Transgenic tomato plants with altered enzyme activities associated with the tricarboxylic acid (TCA) cycle exhibited a substantial reduction in cellulose and lignin during xylem vessel development in the roots [35]. The starch and sucrose metabolism pathway participates in the synthesis of ADP-glucose and GDP-glucose by providing a-D-glucose.
It is worth noting that expression of the c48648.graph_c0 (pectin esterase inhibitor 33 OS), c51200.graph_c0 (vacuolar invertase family protein) and c46943.graph_c0 (Probable fructokinase-5 OS) genes is almost undetectable in S1, and that the expression levels peak in S3. We speculate that these genes play positive roles in regulating xylem development.

Identification of Phytohormone-Related Genes in the Xylem during the Seasonal Change
By analyzing the DEGs and comparing NR database annotations, we can get genes related to phytohormones from six different periods (Supplementary Table S5). Different types of hormones have different expression patterns, and most of them have obvious changes in S3 and S4. Compared with S1, there were 37 DEGs related to plant hormones

Identification of Phytohormone-Related Genes in the Xylem during the Seasonal Change
By analyzing the DEGs and comparing NR database annotations, we can get genes related to phytohormones from six different periods (Supplementary Table S5). Different types of hormones have different expression patterns, and most of them have obvious

Identification of Phytohormone-Related Genes in the Xylem during the Seasonal Change
By analyzing the DEGs and comparing NR database annotations, we can get genes related to phytohormones from six different periods (Supplementary Table S5). Different types of hormones have different expression patterns, and most of them have obvious changes in S3 and S4. Compared with S1, there were 37 DEGs related to plant hormones found in S2, of which 14 were upregulated and 23 were downregulated. In S3, 50 related DEGs were found, and 26 were upregulated and 24 were downregulated. There were 49 upregulated and 26 downregulated DEGs in S4. S5 had 48 DEGs, 25 of which were upregulated and 23 were downregulated. In S6, there were 46 DEGs, 25 of which were upregulated and 21 were downregulated ( Figure 12). During this period, the gene c53437.graph_c0 (auxin-responsive family protein (Populus trichocarpa)) had the largest increase in expression. The only DEGs related to abscisic acid were c43055.graph_c0 (PREDICTED: abscisic acid receptor PYL4 (Populus euphratica)), c50201.graph_c0 (PREDICTED: abscisic acid receptor PYL8-like isoform X1 (Populus euphratica)), c49741.graph_c0 (PREDICTED: abscisic acid receptor PYL8-like (Populus euphratica)) and c16442.graph_c0 (abscisic acid and environmental stress-inducible protein (Jatropha curcas)). The remaining DEGs are all related to auxin, gibberellin, cytokinin and ethylene. This indirectly explained the importance of related genes in the xylem during seasonal changes.

Identification of Transcription Factor (TF) Genes during the Seasonal Change
The secondary cell wall of the xylem is a special structure produced after the differentiation of the procambium. Its thickening process is regulated by a complex multilevel transcription network composed of a large number of transcription factors. These transcription factors mainly include a series of NAC (NAM, ATAF1/2 and CUC2) transcription factors and MYB transcription factors [44]. We found that the relative abundance of mRNAs specific for 44 TF genes showed a high degree of dynamic change during the change of seasons, including two NACs (c56971.graph_c0 and c56580.graph_c0) and nine MYBs (c50601.graph_c1, c27343.graph_c0, c47250.graph_c0, c51583.graph_c1, c48456.graph_c0, c55770.graph_c0, c51063.graph_c1, c55858.graph_c0 and c58751.graph_c0). In addition, there are genes for other transcription factors, such as WRKYs, bZIPs, GATAs and AP2/ERFs. Based on their expression patterns, the differentially expressed TF genes can be divided into four subgroups ( Figure 13) [45]. In subgroup I, compared with the control group S1, the transcription level peaked at S4 and then continued to decline at S5 and S6. The TF genes in subgroup II generally showed a steady decreasing trend over the five time points; the expression of only a few individual genes (c52963.graph_c0, c52341.graph_c0, c35420.graph_c0 and c52941.graph_c0) increased again at S6. Genes in subgroup III had the highest transcript abundance. It is interesting that the TF gene expression peaked in the S1 group. The expression of TF genes in subgroup IV reached a peak at S6. In subgroup IV, there were genes for three WRKY, two

Identification of Transcription Factor (TF) Genes during the Seasonal Change
The secondary cell wall of the xylem is a special structure produced after the differentiation of the procambium. Its thickening process is regulated by a complex multilevel transcription network composed of a large number of transcription factors. These transcription factors mainly include a series of NAC (NAM, ATAF1/2 and CUC2) transcription factors and MYB transcription factors [44]. We found that the relative abundance of mRNAs specific for 44 TF genes showed a high degree of dynamic change during the change of seasons, including two NACs (c56971.graph_c0 and c56580.graph_c0) and nine MYBs (c50601.graph_c1, c27343.graph_c0, c47250.graph_c0, c51583.graph_c1, c48456.graph_c0, c55770.graph_c0, c51063.graph_c1, c55858.graph_c0 and c58751.graph_c0). In addition, there are genes for other transcription factors, such as WRKYs, bZIPs, GATAs and AP2/ERFs. Based on their expression patterns, the differentially expressed TF genes can be divided into four subgroups ( Figure 13) [45]. In subgroup I, compared with the control group S1, the transcription level peaked at S4 and then continued to decline at S5 and S6. The TF genes in subgroup II generally showed a steady decreasing trend over the five time points; the expression of only a few individual genes (c52963.graph_c0, c52341.graph_c0, c35420.graph_c0 and c52941.graph_c0) increased again at S6. Genes in subgroup III had the highest transcript abundance. It is interesting that the TF gene expression peaked in the S1 group. The expression of TF genes in subgroup IV reached a peak at S6. In subgroup IV, there were genes for three WRKY, two heat shock transcription factors, one MYB and one MADS box TF (Supplementary Table S6).

Illumina Sequencing and Functional Annotation Can Provide Valuable Transcriptomi sources for the Analysis of Seasonal Gene Expression
Trees allow the harvest of large quantities of secondary xylem at different dev ment stages in a long period of a growth season. Anatomical observation showed th ferences of xylem development among six time points (Figure 1). Analyzing the ex sion profiles of genes in these stages could reveal the molecular mechanisms behin regulation of xylem development and identify the genes related to different stag wood formation in "Shanxin" poplar. In this study, we used Illumina high-throu nucleotide sequencing to analyze the transcriptome of hybrid Populus by RNA-seq mina sequencing, which is widely used for transcriptome assembly, produces rela short reads compared with some other sequencing technologies but produces higher scriptome coverage at a lower cost. We sequenced 18 transcriptome libraries of six points and obtained a total of 121.93 Gb of clean data. These reads were de novo a bled into 84,733 contigs, of which 17,603 were >1 kb in length. This indicates that t sembled gene fragments are of high quality and may contain many full-length cD These results indicate that our "Shanxin" poplar database will be a valuable transcrip resource for gene discovery and functional analysis (Figure 2).

Illumina Sequencing and Functional Annotation Can Provide Valuable Transcriptomic Resources for the Analysis of Seasonal Gene Expression
Trees allow the harvest of large quantities of secondary xylem at different development stages in a long period of a growth season. Anatomical observation showed the differences of xylem development among six time points (Figure 1). Analyzing the expression profiles of genes in these stages could reveal the molecular mechanisms behind the regulation of xylem development and identify the genes related to different stages of wood formation in "Shanxin" poplar. In this study, we used Illumina high-throughput nucleotide sequencing to analyze the transcriptome of hybrid Populus by RNA-seq. Illumina sequencing, which is widely used for transcriptome assembly, produces relatively short reads compared with some other sequencing technologies but produces higher transcriptome coverage at a lower cost. We sequenced 18 transcriptome libraries of six time points and obtained a total of 121.93 Gb of clean data. These reads were de novo assembled into 84,733 contigs, of which 17,603 were >1 kb in length. This indicates that the assembled gene fragments are of high quality and may contain many full-length cDNAs. These results indicate that our "Shanxin" poplar database will be a valuable transcriptome resource for gene discovery and functional analysis (Figure 2).

Gene Expression in Response to Seasonal Changes
There have been studies to show that some known genes are involved in the regulation of xylem development and secondary cell wall biosynthesis, and examples of the induction of gene expression and inhibition of the corresponding xylem development have been reported, such as MYB [46] and NAC [47] transcription factor genes, lignin and cellulose biosynthesis genes. In our study, the change of seasons resulted in the differential expression of genes in the xylem of "Shanxin" poplar. Among them, using the S1 sample as the control group, 2091, 4014, 4266, 3833 and 3308 genes were upregulated in S2-S6, and 2182, 4268, 4321, 3019 and 2567 genes were downregulated, respectively. In S3, the upregulated genes were the most abundant (Table 3). GO terms related to cell wall formation and thickness in the upregulated genes were also enriched in S1 vs. S2 and S2 vs. S3 (Figure 5). At these time points, the developing xylem elements expand rapidly, and the secondary cell walls are deposited, and various related genes show significant differential expressions, which reflect that a large number of related genes are involved in the active development of xylem in late May to early June in Harbin when the temperature rises rapidly in the spring. Subsequently, the number of DEGs slowly decreased and the downregulated genes related to cell wall formation and thickness were enriched in S4 vs. S5 and S5 vs. S6, indicating that the xylem developmental process gradually slowed, which agrees with the histological observations. These data implied that the xylem development and its gene expression might be correlated with seasonal temperature change. The onset of wood formation (the first cell begins to grow) in Northern hemisphere conifers is primarily driven by the photoperiod and mean annual temperature (MAT) and only secondarily by the spring forcing temperature, winter chilling and moisture availability. The photoperiod interacts with the MAT and plays a dominant role in regulating the onset of secondary meristem (xylem) growth. Understanding the unique relationships between exogenous factors and wood formation could aid in predicting how forest ecosystems respond and adapt to climate warming while improving the prediction accuracy of the earth system model for carbon, water and energy cycles and assess the carbon sequestration potential of forests [48].

Genes for Polysaccharide Biosynthesis and Metabolism Were Helpful for Xylem Development during Seasonal Change
The process of xylem development is accompanied by polysaccharide biosynthesis during the formation of cell walls. The synthesis of glucoxylan and cellulose relies on various sugar metabolism pathways. We found that, in the progression from stages S1 to S6, the DEGs involved in sugar metabolism showed a certain trend in the different periods, and the highest number of DEGs was in S4, according to the differential expression and GO enrichment analyses (Table 3). In higher plants, cellulose is synthesized at the plasma membrane by the cellulose synthase (CesA) enzyme complex [49]. Sucrose synthase (SuSy) is one of the key enzymes in plant sucrose metabolism, and it is ubiquitous in most plant tissues [50]. SuSy participates in many metabolic processes in plants, including the synthesis of starch and cellulose, as well as the distribution of carbon sources [51]. In this study, the expression of CesA-and SuSy-related genes were upregulated at S3 (Figures 8 and 9), indicating that cellulose biosynthesis was an important bioprocess during this period of xylem development. Some of the cellulose synthase genes continue to show an upregulated expression at S4 compared with S1; after which, the expression levels slowly decrease. Among these CESA genes, two (c57251.graph_c0 and c55736.graph_c1) are homolog to PtrCESA4 and PtrCESA7 of Populus trichocarpa wood and c32828.graph_c0 is a homolog to PtrCESA8, which suggested that they were involved in the cellulose biosynthesis of the secondary cell wall during this period (Supplementary Table S7). The study showed that several glycosyltransferases are involved in the synthesis of glucoxylan (GX) and cellulose in poplar trees [52]. We also found that there were 101 different genes related to glycosyltransferase in S4 compared with S1. A probable polygalacturonase gene and a sucrose-phosphate synthase were upregulated compared to S1 and reached a peak at S4. The polygalacturonase gene has been proven to affect the hydraulic properties of the poplar xylem [53]. Sucrose-phosphate synthase expression influences the poplar phenology [54]. The obvious changes in the DEGs in sugar metabolism pathways, cellulose synthesis and sucrose synthase indicate that S2, S3 and S4 are important stages in xylem development, and these genes play important roles in that process.

Expression of Genes Related to Lignin Synthesis Varied in the Different Stages of Xylem Development
The xylem contains a large amount of lignin [55], which is related to the accumulation of lignin during the development of the xylem. The biosynthesis of lignin is a complex process that explains the growth state of xylem at different developmental periods. Gene Ontology enrichment revealed an early and a late upregulation of the genes associated with the phenylpropanoid pathway among the six time points of xylem development of the Shanxin poplar. It has been previously shown that lignin biosynthesis is the main phenylpropanoid product during wood development and that it can occur even in the walls of dead xylem tracheary elements and possibly fibers, thanks to the neighboring, still living fibers and xylem parenchyma cells (e.g., ray cells) [56][57][58]. The molecular processes involved in monolignol biosynthesis are not restricted to the cambial activity timeframe but continue after the end of cambium cell proliferation. In temperate regions, latewood is produced when cambial activity declines with the approach of autumnal dormancy. The temporal pattern analysis showed an accumulation pattern of the genes involved in the phenylpropanoid pathway during latewood formation [59]. The protein identification of the seasonal variations of xylem also identified the proteins involved in lignin biosynthesis in mature xylem [10]. The late upregulation of the genes in the phenylpropanoid pathway might play roles in these tissues and cells. We analyzed the 14 intermediate product-related genes involved and found that most of the DEGs showed an upregulated expression trend throughout the period, and the expression levels peaked at S4 compared to S1 (Supplementary Table S4). This shows that lignin synthesis reflects the changes in the xylem during its development. The enrichment analysis of KEGGs and DEGs in the transcriptome sequencing of the "Shanxin" poplar revealed many genes that regulate lignin synthesis. The expression levels of the genes encoding CAD2 and CAD3 increased at either S5 or S6 ( Figure 10A). Two CCO and one CAD genes were downregulated at S4 compared to S1 ( Figure 10B). The gene c51274.graph_c0 (cinnamoyl CoA reductase), which promotes the formation of feruloyl-CoA and sinapoyl-CoA to coniferyl aldehyde and sinapaldehyde, was upregulated at S4 compared to S1 ( Figure 10C). Two PAL, one 4CL and two CAD genes were upregulated from S2 to S6 compared to S1 ( Figure 10D); for example, the gene c54100.graph_c0 (PAL1) reached its peak at S4, and its expression was significantly enriched (Supplementary Table S4). Combined with the results of the histological observation, these genes play roles in lignin biosynthesis in addition to xylem development. This can also be seen in the KEGG analysis of the lignin synthesis pathway in genes such as c57344.graph_c0 (2-dehydro-3-deoxyphosphoheptonate aldolase family protein) and c53880.graph_c0 (coniferyl aldehyde 5-hydroxylase 2). A coniferyl aldehyde 5hydroxylase in aspen (P. tremuloides) was demonstrated to regulate monolignol biosynthesis in the xylem [60]. In the present study, the expression level of the gene coniferyl aldehyde 5-hydroxylase 2 is very interesting; the RPKM values at S1 were 0.7, 0.6 and 0.89, and, at S2, they were 3.85, 3.18 and 5.62. However, the RPKM values at S3 reached an astonishing 944, 1103.07 and 883.1, and, at S4, they were 1769.84, 1939.91 and 1985.57, reaching their highest levels at this time point (Supplementary Table S4). It is worth noting that the expression of gene c51159.graph_c0 (aldehyde dehydrogenase family 2 member C4) was downregulated, which directly affects the synthesis of coniferyl aldehyde and sinapaldehyde into ferulic acid and sinapic acid [61]. More interesting is that ferulic and sinapic acids are the preproducts of feruloyl-CoA and sinapoyl-CoA. The sequencing results show that these two genes must share some kind of mutual relationship. Our results suggest that these genes are involved in the response of xylem development to the seasonal changes in "Shanxin" polar.

Phytohormone Signaling during the Seasonal Change
The plant hormones auxin and gibberellin play important roles in the development of xylem. Our transcriptome data shows that the expression levels of the signal transduction pathway-related genes related to several plant hormones were significantly different at the different time points (Supplementary Table S8). The GO enrichment analysis of the DEGs showed that terms related to hormones were enriched in the S1 vs. S2, S2 vs. S3 and S4 vs. S5 comparisons, and both of the upregulation and downregulation genes in GO terms relate to hormones were enriched in S1 vs. S2. The KEGG enrichment analysis also showed the presence of the hormone signal pathway, which suggests that hormone signals play important roles at these stages. Auxin is a key regulator in the growth of the vascular cambium and in xylem development [62]. It has long been known that the activation of preprotocells depends on auxin signaling and transport. According to the auxin tube theory, the initial auxin flux from the auxin source to the auxin pool creates its own transport channel. There were 29 DEGs related to auxin among the six time points (Supplementary  Table S5). The RPKM value of the gene c53437.graph_c0 (homolog to the auxin-responsive family protein (Populus trichocarpa)) related to the auxin signal transduction pathway was 0.15 at S1, and it was upregulated to 27.60 at S2. At S4, the expression peaked at 281.1, which indicates that auxin is involved in xylem development. Gibberellin (GA) is a key signaling molecule that induces vessel growth, fiber differentiation and lignification. Gene c45848.graph_c0 (Gibberellin 2-beta-dioxygenase OS) is the most abundant gene in the gibberellin pathway in our libraries. The study showed that the accumulation of Gibberellin 2-beta-dioxygenase maintains the hypocotyl elongation [63]. However, the molecular mechanism of how GA affects the xylem elongation and secondary wall development in woody tree species is still unclear. It has been shown previously that the GA signal in the xylem mediates the occurrence of lignin and promotes fiber elongation [64]. There were 12 DEGs related to GAs among the six time points. The sequencing results showed that the expression of gibberellin-related genes was significantly different in the different periods, and the expression was the most significant at S4. The results indicate that gibberellin may play an important role in xylem development.

Transcription Factors Involved in Xylem Development
The transcription of genes encoding CesA, PAL and other enzymes related to the biosynthesis of cellulose, xylan and lignin is regulated by the transcription factors NAC and MYB [65][66][67]. SND and MYB46 are thought to play a key role in these processes [68]. In this study, many TFs were found to play an important role in xylem development; this included 44 genes predicted to encode TFs such as NACs, MYBs, WRKYs, bZIPs, GATAs and AP2/ERFs that were differentially expressed at the six time points (Figure 13). These transcription factor genes could be clustered into different groups at different stages of xylem development in a year, suggesting that each group of genes played an important role in the corresponding stages ( Figure 13). The largest category belonged to the WRKY and MYB TF families, both of which consisted of eight members. At present, the involvement of WRKY transcription factors in the regulation of lignification and xylem development has only been proposed indirectly. The transcription profile analysis showed that some WRKY TF genes are differentially expressed during protoplast cell wall regeneration, secondary cell wall biosynthesis in the xylem along the axial development gradient or during the formation of tension wood [69]. Together, with these results, we can hypothesize that WRKY TFs play roles in the xylem development, especially in the autumn xylem development process, because some WRKY transcription factor genes were enriched in the S6 stage (Supplementary Table S6). In vivo assays showed that MYB46 transcriptionally activates downstream target genes for transcription factors, three of which (AtC3H14, MYB52 and MYB63) were shown to be able to activate secondary wall biosynthesis genes [70]. The R2R3-MYB family, one of the largest transcription factor families in higher plants, controls a wide variety of plant-specific processes, including, notably, phenylpropanoid metabolism and secondary cell wall formation [71]. MYB TFs can regulate the expression of the genes involved in the synthesis of phenylpropanes, thereby affecting the content of lignin and the accumulation of phenols [72][73][74]. In our xylem RNA-seq libraries, the expression of MYB genes showed different patterns. For example, the expression level of c50601.graph_c1 (homolog to the myb family transcription factor family protein (Populus trichocarpa)) in S2 was 4.85, and it was 19.98, 18.83, 14.17 and 0.67 in stages S3, S4, S5 and S6, respectively. Gene c56971.graph_c0(NAC domain-containing protein 7 OS) has the highest expression level at S3 and is similar to c50601.graph_c1 (Supplementary Table S6). The c56971.graph_c0 homolog to PtVNS/PtrWND in poplar, which induced ectopic secondary wall thickening in poplar leaves [75]. These are basically the same trend found in the expression of some lignin synthesis genes, which suggests that these genes are involved in the xylem development response to seasonal changes. We find that, in addition to the known differences in the expression levels of genes for transcription factors such as WRKY, MYB and NAC, many genes encoding transcription factors similar to LTM and BATA have also been detected, which implies their roles in xylem development in "Shanxin" polar.

Conclusions
In this study, using high-throughput RNA sequencing, we constructed a "Shanxin" poplar xylem transcriptome data set containing 84,733 predicted transcripts. Based on the slices of xylem development at different time points, further studies on six time points taken from April to August showed the dynamic changes in gene expression that occur during a growing season. By analyzing the different genes in the growing season, including cell expansion and cell wall biosynthesis, polysaccharide biosynthesis and metabolism, ligninrelated genes, plant hormone-related genes and transcription factor (TF) genes, as well as the analysis of different enrichment pathways to understand its main biological processes, it is speculated that it may play an important role in the development of xylem. The transcriptome assembly and gene expression profile analysis of xylem development will provide data and materials that can be used in the molecular breeding of "Shanxin" poplar.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/f12040451/s1: Table S1: COG and eggNOG Function classification. Table S2: GO classification enrichment analysis. Table S3: The table of glucose metabolism differential genes. Table S4: The table of lignin synthesis differential genes. Table S5: The table of hormone differential genes. Table S6:  The table of