Transcriptome Analysis of Cambium Tissue of Paulownia Collected during Winter and Spring

: Paulownia ( Paulownia elongata ) is a fast-growing, multipurpose deciduous hardwood species that grows in a wide range of temperatures from –30 ◦ C to 45 ◦ C. Seasonal cues inﬂuence the secondary growth of tree stems, including cambial activity, wood chemistry, and transition to latewood formation. In this study, a de novo transcriptome approach was conducted to identify the transcripts expressed in vascular cambial tissue from senescent winter and actively growing spring seasons. An Illumina paired-end sequenced cambial transcriptome generated 297,049,842 clean reads, which ﬁnally yielded 61,639 annotated unigenes. Based on non-redundant protein database analyses, Paulownia cambial unigenes shared the highest homology (64.8%) with Erythranthe guttata . KEGG annotation of 35,471 unigenes identiﬁed pathways enriched in metabolic activities. Transcriptome-wide DEG analysis showed that 2688 and 7411 genes were upregulated and downregulated, respectively, in spring tissues compared to winter. Interestingly, several transcripts encoding heat shock proteins were upregulated in the spring season. RT-qPCR expression results of ﬁfteen wood-forming candidate genes involved in hemicellulose, cellulose, lignin, auxin, and cytokinin pathways showed that the hemicellulose genes ( CSLC4, FUT1, AXY4, GATL1 , and IRX19 ) were signiﬁcantly upregulated in spring season tissues when compared to winter tissues. In contrast, lignin pathway genes CCR1 and CAD1 were upregulated in winter cambium. Finally, a transcriptome-wide marker analysis identiﬁed 11,338 Simple Sequence Repeat (SSRs). The AG/CT dinucleotide repeat predominately represented all SSRs. Altogether, the cambial transcriptomic analysis reported here highlights the molecular events of wood formation during winter and spring. The identiﬁcation of candidate genes involved in the cambial growth provides a roadmap of wood formation in Paulownia and other trees for the seasonal growth variation.


Introduction
Paulownia (Paulownia elongata) is an extremely fast-growing woody plant growing up to 20 feet in one year when young. Some Paulownia spp., when in plantation, can be harvested for saw timber in as little as five years. The genus Paulownia consists of nine species of deciduous, fast-growing, multi-purpose, hardwood trees [1] that have long been shown to be extremely adaptive to wide environmental variations in both edaphic and climatic factors, as well as being capable of growing on marginal lands [2,3]. Species of Paulownia are native to Asia and are widely cultivated in China, Laos, Vietnam, Japan, and Korea. It has now been introduced and cultivated in Australia, Europe, and North and Central America. Ten-year-old Paulownia trees, in natural conditions, can attain 30-40 cm in diameter at breast height (DBH) and a timber volume of 0.3-0.5 m 3 [1]. Craftworkers in Japan and other countries have used this valuable wood to create intricate carvings, surfboards, musical instruments, toys, and furniture. Paulownia wood has a high ignition point of 420-430 • C compared to other hardwoods, which generally range from 220-225 • C, thus making Paulownia wood fire retardant [4,5]. Paulownia plants bear abundant flowers that are highly nectariferous and yield premium honey [6], adding to the rural economy. By adding Paulownia wood flour (25-40%) to plastics, an attractive, equally strong, environmentally agreeable, and economically important biocomposite can be produced [7][8][9] to serve many industries. Additionally, the lightweight and strong nature of the wood makes it widely applicable in the music industry for preparing the finer components of instruments. Biochar produced from Paulownia is also a desirable organic soil amendment that allows the growth of beneficial microbes in the porous holes of the biochar [10]. Recently, researchers found its potential use as an animal feed resource [11].
Wood synthesis provides one of the most important sinks for atmospheric carbon dioxide [12]. Wood formation is a result of the regulated accumulation of secondary xylem cells (fibers, vessels, and rays in dicots) differentiated from the vascular cambium that involves wall thickening. This wall thickening is accompanied by the biosynthesis of wall components, lignin, cellulose, and hemicelluloses and is terminated by programmed cell death [13,14]. In order to survive multiple growing seasons, perennial plant species have adapted a dormancy regulation system that allows active growth during the desirable time of year and vegetative dormancy when climatic conditions are unfavorable for growth [15]. The Paulownia tree, due to its fast-growing nature, is capable of producing~45 kg/tree in the first growing year and~90 kg/tree at the end of the second year [16]. Being a perennial tree, Paulownia harvest is not limited to a small seasonal window but can be conducted year-round with proper management practices. Another beneficial property of Paulownia is coppicing, which is defined as the production of multiple sprouts from a stump after the removal of the tree or shrub. Harvest cycles of 2-3 years could be implemented to establish a fast-growing bioenergy crop. Since Paulownia is a short-rotation fast-growing perennial tree and serves as a good candidate to produce lignocellulosic biofuel, which can eliminate dependence on fossil fuel.
Transcriptome analyses from various tree species indicate that the putative role of gene families belonging to receptor kinases, transcription factors, and secondary wall biosynthesis are highly expressed in wood-forming cells [17][18][19][20][21]. At a molecular level, studies related to Paulownia gene expression profiling published in the last few years have primarily focused on drought tolerance and host-pathogen interactions [22][23][24][25]. Cambial development (the initiation and activity of the vascular cambium) leads to an accumulation of wood (secondary xylem tissue). Seasonal cues play a significant role in determining cambial growth as perennial plants growing in temperate and high-latitude regions show termination of cell division in the meristems [26] and reversal of growth arrest during long days [27]. Time-coursed RNA sequence studies identified downregulated phytohormonerelated genes (IAA, ARF, and SAURs) and upregulated circadian genes (PIF3 and PRR5; [28,29]. A transcriptome-wide profiling study in P. elongata identified a subset of candidate genes that contribute to the production of wood [30] by investigating the differential expression of transcripts of the vascular cambium. Furthermore, there is evidence for the microR-NAs controlling gene expression in Paulownia tomentosa cambial tissues in response to seasonal changes [31]. Transcriptomic analyses have been carried out to profile gene expression regulations for biotic and abiotic stresses and growth responses. However, to the best of our knowledge, no study has described how the gene expression profile changes in woody tissues under seasonal variations. In this study, we sequenced and analyzed the transcriptomes of cambium tissues collected during the winter and spring seasons to assess the impact of two seasons on biomass. A transcriptome-wide analysis identified 61,639 annotated unigenes, and 2688 and 7411 transcripts were up-and downregulated, respectively, in the spring season. Interestingly, among the selected wood-forming genes, hemicellulose-specific genes were upregulated in spring. Finally, 11,338 simple sequence repeats (SSRs) were identified from the transcriptome data. The identification of genes and pathways involved in cambial growth will be useful to further investigate the regulation of wood formation in Paulownia and other trees.

Tissue Sampling and RNA Isolation
Samples were selected from trees growing at the Fort Valley State University (FVSU) Paulownia Bioenergy Farm located at 32 • 31 15.04 N and 83 • 52 12.95 W. Paulownia trees bear flowers first, and after 3-4 weeks of flowering, leaves appear. Samples in replicates were collected from two seasonal points, each representing a different physiological state. The first sample (winter wood; hereafter referred to as WW) was collected before flowering on 1 March 2015 (min temp 33.1 • F and max temp 54 • F) and represented the senescent winter wood ( Figure 1A). The second sample (spring wood; hereafter referred to as SW) was collected on 20 May 2015 (min temp 64 • F, max temp 91 • F) during spring and was representative of the actively growing spring wood ( Figure 1B). Samples were harvested from twigs located at approximately 1.0-1.5 m from the ground with an average diameter of 2.5 cm. Since Paulownia has an opposite branching pattern, the WW and SW samples were taken from the mirrored axis of nodes. The samples were labeled, wrapped in aluminum foil, flash-frozen in liquid nitrogen, and subsequently stored in a −80 • C freezer until further use. Biological replicates were labeled as WW1, WW2, and WW3 for WW, and SW1, SW2, and SW3 for SW, respectively. impact of two seasons on biomass. A transcriptome-wide analysis identified 61,639 annotated unigenes, and 2688 and 7411 transcripts were up-and downregulated, respectively, in the spring season. Interestingly, among the selected wood-forming genes, hemicellulose-specific genes were upregulated in spring. Finally, 11,338 simple sequence repeats (SSRs) were identified from the transcriptome data. The identification of genes and pathways involved in cambial growth will be useful to further investigate the regulation of wood formation in Paulownia and other trees.

Tissue Sampling and RNA Isolation
Samples were selected from trees growing at the Fort Valley State University (FVSU) Paulownia Bioenergy Farm located at 32° 31′15.04″ N and 83° 52′12.95″ W. Paulownia trees bear flowers first, and after 3-4 weeks of flowering, leaves appear. Samples in replicates were collected from two seasonal points, each representing a different physiological state. The first sample (winter wood; hereafter referred to as WW) was collected before flowering on March 1, 2015 (min temp 33.1 °F and max temp 54 °F) and represented the senescent winter wood ( Figure 1A). The second sample (spring wood; hereafter referred to as SW) was collected on 20 May 2015 (min temp 64 °F, max temp 91 °F) during spring and was representative of the actively growing spring wood ( Figure 1B). Samples were harvested from twigs located at approximately 1.0-1.5 m from the ground with an average diameter of 2.5 cm. Since Paulownia has an opposite branching pattern, the WW and SW samples were taken from the mirrored axis of nodes. The samples were labeled, wrapped in aluminum foil, flash-frozen in liquid nitrogen, and subsequently stored in a −80 °C freezer until further use. Biological replicates were labeled as WW1, WW2, and WW3 for WW, and SW1, SW2, and SW3 for SW, respectively. For high-quality, intact total RNA extraction, vascular cambium tissues were harvested from the frozen samples by first slicing a shallow, longitudinal cut into the outer bark with a sterile scalpel ( Figure 1C). The bark was then removed using sterile forceps in a large, single piece. The frozen green vascular cambium was then gently scraped from the wood below into small strips using a sterile scalpel. A 100 mg vascular cambium tissue sample was powdered into microvials containing zirconia beads (BioSpec, Bartlesville, OK, USA) and 550 µL of TRIzol reagent (Invitrogen, Carlsbad, CA, USA) in a MagNA Lyser (Roche, Basel, Switzerland) as described earlier [32]. RNA was purified using Direct-Zol™ RNA mini-prep kit (Zymo Research, Irvine, CA, USA), and removal of genomic DNA contamination was carried out using DNase I treatment. RNA quality and quantity were analyzed using NanoDrop 1100 (NanoDrop, Wilmington, DE, USA) and Agilent 2100 Bioanalyzer (Agilent, Santa Clara, CA, USA). For high-quality, intact total RNA extraction, vascular cambium tissues were harvested from the frozen samples by first slicing a shallow, longitudinal cut into the outer bark with a sterile scalpel ( Figure 1C). The bark was then removed using sterile forceps in a large, single piece. The frozen green vascular cambium was then gently scraped from the wood below into small strips using a sterile scalpel. A 100 mg vascular cambium tissue sample was powdered into microvials containing zirconia beads (BioSpec, Bartlesville, OK, USA) and 550 µL of TRIzol reagent (Invitrogen, Carlsbad, CA, USA) in a MagNA Lyser (Roche, Basel, Switzerland) as described earlier [32]. RNA was purified using Direct-Zol™ RNA mini-prep kit (Zymo Research, Irvine, CA, USA), and removal of genomic DNA contamination was carried out using DNase I treatment. RNA quality and quantity were analyzed using NanoDrop 1100 (NanoDrop, Wilmington, DE, USA) and Agilent 2100 Bioanalyzer (Agilent, Santa Clara, CA, USA).

cDNA Synthesis and RNA Sequencing
RNA samples of each biological replicate from both treatments (a total of 6 samples) were sequenced at BGI International (http://bgi-international.com/). Briefly, magnetic beads coated with Oligo (dT) were used to isolate mRNA from the total RNA, and size selected mRNA were adapter-ligated and sequenced on the Illumina Hiseq 2000 platform following the manufacturer's protocol.

Assembly and Annotation
Raw reads generated from pair-end sequencing were filtered for adapters and lowquality sequences by BGI proprietary software "filter_fq". Clean reads were assembled into transcripts using Trinity (http://trinityrnaseq.sourceforge.net/) [33]. Three modules in Trinity (Inchworm, Chrysalis, and Butterfly) were applied sequentially to process raw RNAseq reads to assemble them into contigs and full-length transcripts (unigenes). Butterfly also determined alternatively spliced isoforms of genes Unigenes were aligned with five databases, namely, KEGG (Kyoto Encyclopedia of Genes and Genomes), COG (Clusters of Orthologous Groups), NT (NCBI nucleotide database), NR (NCBI non-redundant protein database), and Swiss-Prot (Protein sequence database). The KEGG database [34] was used to perform a systematic analysis of metabolic pathways and the function of gene products within a cell. The COG database (http://www.ncbi.nlm.nih.gov/COG/) was used to classify orthologous gene products into clusters. COG clusters predicted the possible function of the transcripts. Both NCBI's NR database and the Swiss-Prot (http://www. uniprot.org/uniprot/) annotated protein databases and added additional information about the possible function of the transcripts.

Gene Ontology and Coding Sequences
Gene ontology (GO) was employed to standardize functional gene classification, including molecular function, cellular components, and biological processes. The Blast2GO program (http://www.blast2go.com/b2ghome) was used to retrieve GO functional classification for all transcripts [35]. In order to determine the CDS for the transcripts, unigenes were first aligned to the protein databases, listed in order of priority of NR, Swiss-Prot, KEGG, and COG by using a local blastx (ftp://ftp.ncbi.nlm.nih.gov/blast/executables/ release/LATEST/), with a significant cut-off value of e <0.00001 , of the unigene sequences. Unigenes with alignments to higher priority databases, for example, the NR database, were not aligned to lower priority databases. The highest-ranking proteins in the blastx results were used to decide the coding region sequences of unigenes. Results of the blastx alignment used a standard codon table to translate the nucleotide query sequence into a translated amino acid sequence. Unigenes that could not be aligned to any database were further scanned by ESTScan [36], producing nucleotide sequence (5 →3 ) direction and amino sequence of the predicted coding region.

Gene Expression Analysis
In order to determine the expression pattern of the unigenes, clean reads were first mapped to unigenes using the program Bowtie2 (v. 2.2.5) [37]. SAM files generated through Bowtie were used with the RSEM (RNA-Seq by Expectation-Maximization) software package (http://deweylab.github.io/RSEM/; v1.2.12) in R (v1.03; http://www.r-project.org/) to measure the expression level of unigenes. RSEM software was used to estimate gene expression levels from RNA-seq data [38] in FPKM (fragments per kilobase of transcript per million mapped reads) format, which was subsequently used to perform differential gene expression analysis in this study. To detect differentially expressed genes (DEGs), the program NOIseq was utilized [39]. In this study, WW and SW unigenes with a fold change of ≥2 and a probability ≥0.8 were considered to be significantly differentially expressed. Principal component analysis (PCA) was performed using the R statistical package.

Validation of Wood-Forming Candidate Genes with RT-qPCR
Based on existing research information, fifteen wood-forming candidate genes corresponding to the biosynthesis of cellulose, hemicellulose, or lignin [18,[40][41][42][43], were selected for validation of expression level. Using the local blast utility (ftp://ftp.ncbi.nlm.nih.gov/ blast/executables/release/LATEST/), a database of all Paulownia unigenes was created. The mRNA sequences acquired for each of the selected genes were then aligned to the database of all unigenes using the local blastx utility. The unigenes showing maximum homology for each of the genes were selected for a two-step RT-qPCR primer design. The software Primer Express v3.01 (Applied Biosystems, Foster City, CA) was used to design primers for the unigenes corresponding to the selected wood forming genes. cDNA was synthesized from Paulownia vascular cambium total RNA using SuperScript™ II Reverse Transcriptase (Invitrogen, Waltham, Massachusetts, USA) with the suggested protocol. FastStart SYBR Green Master Mix (Roche, Grenzacherstrasse, Basel, Switzerland) reagent was used in combination with primers and cDNA. RT-qPCR of three biological replications with no-template control (NTC) involved StepOnePlus Real-Time PCR System (Applied Biosystems, Foster City, CA, USA) and FastStart SYBR Green (Roche). The expression of selected genes was normalized with reference gene 18S rRNA. Finally, the relative gene expression was measured using the 2 -∆∆Ct method [44].
All raw files of sequencing data have been submitted to the NCBI Sequenced Read Archive (SRA) database with accession numbers SRX9298494-SRX9298499.

RNA-Seq and Transcriptome Assembly of Paulownia Cambial Tissue
To obtain the candidate genes associated with cambium development of the empress tree (common name for all Paulownia species) during seasonal growth, transcriptome sequencing analysis for tissues representing winter and spring seasons ( Figure 1A,B) was carried out by collecting cambial tissues from tree twigs ( Figure 1C). A total of 305,882,370 (~29 Gb) raw reads were generated. Removal of adapter sequences, lowquality reads, and ambiguous sequences resulted in 297,049,842 clean reads (Q20 > 97.73%) with an average length of 100 nucleotides. WW generated more raw reads when compared to spring wood samples ( Table 1). The de novo assembling of clean reads produced 129,428 and 104,388 contigs for WW and SW cambial tissues, respectively. Clustering and assembly of these contigs resulted in 64,142 and 45,671 unigenes for WW and SW with an average length of 960 and 842 nucleotides, respectively. Among the unigenes, all unigenes were sub-classified according to nucleotide length ( Figure S1). A total of 40,814 unigenes were greater than 1 Kbs, and 58,654 unigenes were greater than 500 nucleotides in length. The mean length of the contigs (~340 bp) was shorter than that of the unigenes (>1000 bp) in most previously reported studies. The paired-end reads resulted in longer unigenes (average of~900 bp) than those reported in previous transcriptome studies on trees [45,46]. The mean length of unigenes (900 nucleotides) was less than those reported in previous studies on tetraploid Paulownia australis and drought exposed P. tomentosa [23,47]. Most of our assembled unigenes showed homology to nucleotide sequences deposited in six public nucleotide databases. The unmatched unigenes are most likely to represent Paulownia-specific genes especially related to the winter and spring seasons.

Functional Annotation of Paulownia Cambial Transcriptome
A total of 61,639 transcripts were annotated by performing a BLAST search of the sequences across six databases as described above. The BLAST search using 61,639 transcripts showed 72.47% similarity to non-redundant (Nr) protein database (Tables S1 and S2).
Of the annotated sequences in the Nr database, 39.3% of the mapped unigenes had very significant homology to known sequences (e-value 10-100), 35.1% showed significant homology (10-100; e-value 10-30), and 25.6% showed weak homology (e-value 10-30 to 10-5) (Figure 2A). We also performed the sequence conservation analysis of Paulownia transcripts with proteomes of all sequenced plant species. As depicted in Figure 2B, the E-value distribution analysis of transcripts showed that 47.0% of unigenes had a similarity of more than 80% with plant species. We used a BLAST search to study the relationship of Paulownia with other plant species to identify proteins and pathways that would be unique to Paulownia. The sequence conservation analysis of transcripts showed homology to nucleotide sequences from Erythranthe guttata (64.8%), followed by Vitis vinifera (6.3%), Solanum tuberosum (5.0%), Solanum lycopersicum (3.0%), Theobroma cacao (2.7%), and others ( Figure 2C). Erythranthe guttata, a yellow bee-pollinated annual or perennial plant, is a model organism for biological studies. Paulownia transcripts shared strong homology with Erythranthe species, and this could be due to close phylogenetic relationships between these two species [48]. However, transcripts from a drought-related transcriptomic study of Paulownia [28,47] showed homology to Vitis vinifera (45-48%), which could be due to the selection of seasonal abiotic tissue for cambial transcriptome study.
Further analysis of 43,780 transcripts using the COG database classified these transcripts into 25 different protein families with potential functions in transcription, replication and recombination, posttranslational modification, signal transduction, and others.
Interestingly, only a few transcripts exhibited their potential role in the extracellular structure and nuclear structure (17 and 4 transcripts, respectively). Importantly, many unigenes have been assigned to a wide range of COG classifications (Figure 3), indicating that a wide diversity of transcripts were involved in wood formation as in Chinese fir [31].
The Gene Ontology (GO) classification classified 42,588 out of 62,639 transcripts into ontologies related to molecular functions, cellular components, and biological processes ( Figure 4). We identified a significantly higher number of transcripts involved in metabolic processes (19,418) and related to cellular processes (18,047) when compared to others, such as rhythmic processes. The most represented category for cellular components was cells (GO: 0005623; 18,186 genes), followed by organelles (GO:0043226; 14,053 genes). Genes and pathways putatively responsible for dormant winter and active spring growth in Paulownia were identified in this study. In Populus, PtrHB7, a class III HD-Zip gene, is known to play a critical role in the regulation of vascular cambium differentiation [49] and homeobox gene ARBORKNOX1 regulates the shoot apical meristem and the vascular cambium [50]. In our study, Unigene2201, Unigene3374, and Unigene4121 which were downregulated in their expressions, belong to GO:0005488 (molecular function: binding) and are homologs of the KNOX gene. KNOX family gene KNAT7 negatively regulates secondary wall formation in Arabidopsis and Populus [51]. Since KNAT7 is a negative regulator of secondary wall biosynthesis, these Paulownia homologs might positively regulate cambium growth during the active spring season. Further analysis of 43,780 transcripts using the COG database classified these transcripts into 25 different protein families with potential functions in transcription, replication and recombination, posttranslational modification, signal transduction, and others.
Interestingly, only a few transcripts exhibited their potential role in the extracellular structure and nuclear structure (17 and 4 transcripts, respectively). Importantly, many unigenes have been assigned to a wide range of COG classifications (Figure 3), indicating that a wide diversity of transcripts were involved in wood formation as in Chinese fir [31].  (Tables  S3 and S4). With the help of the KEGG database, we could further analyze the metabolic pathways and functions of gene products, which can help in studying the complex biological behaviors of genes. The majority of representative unigenes were annotated to metabolic pathways, biosynthesis of secondary metabolites, plant-pathogen interaction, plant hormone signaling, spliceosome, and phenylpropanoid biosynthesis using the KEGG database, which led us to conclude that most of the genes identified in this study are involved in cambial differentiation and wood formation.  The Gene Ontology (GO) classification classified 42,588 out of 62,639 transcripts into ontologies related to molecular functions, cellular components, and biological processes ( Figure 4). We identified a significantly higher number of transcripts involved in metabolic processes (19,418) and related to cellular processes (18,047) when compared to others, such as rhythmic processes. The most represented category for cellular components was cells (GO: 0005623; 18186 genes), followed by organelles (GO:0043226; 14053 genes). Genes and pathways putatively responsible for dormant winter and active spring growth in Paulownia were identified in this study. In Populus, PtrHB7, a class III HD-Zip gene, is known to play a critical role in the regulation of vascular cambium differentiation [49] and homeobox gene ARBORKNOX1 regulates the shoot apical meristem and the vascular cambium [50]. In our study, Unigene2201, Unigene3374, and Unigene4121 which were downregulated in their expressions, belong to GO:0005488 (molecular function: binding) and are homologs of the KNOX gene. KNOX family gene KNAT7 negatively regulates secondary wall formation in Arabidopsis and Populus [51]. Since KNAT7 is a negative regulator of secondary wall biosynthesis, these Paulownia homologs might positively regulate cambium growth during the active spring season.   (Tables S3 and S4). With the help of the KEGG database, we could further analyze the metabolic pathways and functions of gene products, which can help in studying the complex biological behaviors of genes. The majority of representative unigenes were an-

Transcriptional Profiling of Cambial Tissues in Winter and Spring
A total of 10,099 (12.33%) transcripts were found to be significantly differentially expressed between two tissue samples. Of these differentially expressed genes (DEGs), 2688 (26.61%) transcripts were found to be upregulated (>1.6 fold) in the spring season, whereas 7411 (73.39%) were downregulated (<−1.6 fold) when compared to the winter season ( Figure S2). Hierarchical clustering of the DEGs identified in winter and spring conditions led to the detailed overall structure of clustering. This indirectly indicated that more genes were upregulated, active, and required during the senescent winter season to keep tissues dormant. Out of 2688 genes, the top 20 genes with log 2 Fold change >8.00 are summarized in Table 2. This included APC/C cyclosome complex, phosphoenolpyruvate carboxykinase, different classes of heat shock proteins, actin depolymerization factor, anaphase-promoting complex subunit (>12-fold expression), etc. Similarly, many key genes, including synthases such as galactinol synthase (<−12-fold expression), rosmarinate synthase, and valencene synthase, kinases such as receptor-like protein kinase, serine/threonine-protein kinase, CBL-interacting protein kinase, and hormone-regulated genes such as auxin efflux carrier family protein and ethylene-responsive transcription factor were downregulated ( Table 3). The cell cycle is one of the most important biological processes in the cambial zone and plays a central role in regulating the growth and development of organisms, including plants. The anaphase-promoting complex/cyclosome (APC/C; homolog in our study Unigene8688), a well-known ubiquitin ligase, acts to accomplish basic cell-cycle control. The APC/C must be turned off at the end of the G1 phase to allow the S phase cyclins to accumulate and cells to begin DNA replication [52]. This is very key during the spring season for cell multiplication and growth. The Cyclin U2 (Unigene22553), one of the major cyclins involved in cell cycle control, like cyclins A and B on maximum gene expression in the poplar cambium zone [53], was upregulated in Paulownia. The high abundance of cyclin transcripts in active cambium during the spring season also reflected a positive correlation between cambium cell division and key cell cycle gene expression.  As shown in Figure 5, GO analysis indicated that most of the differentially expressed transcripts for biological processes were involved in the metabolic process (1016), cellular process (890), and response to stimulus (350). Similarly, GO cellular component analysis revealed that a large number of transcripts encoded functions related to cell (824), cell part (789), membrane (436), and organelle (399) synthesis. Meanwhile, GO molecular function analysis showed that the DEGs predominantly contributed to catalytic activity (836), followed by transporter activity and structural molecule activity. KEGG enrichment analysis of DEG showed that these genes were involved in various pathways in the Paulownia plant during seasonal changes ( Table 4). Most of the DEGs were enriched in metabolic pathways (ko01100; 1387), biosynthesis of secondary metabolites (ko01110; 827), and plant hormone signal transduction (ko04075; 320). It was also found that starch and sucrose metabolism (ko00500; 172; Figure S3) and phenylpropanoid biosynthesis (ko00940; 106; Figure S4), which correspond to the production of several key wood-forming genes, were within the top 25 most DEG enriched KEGG pathways (Table 4). Nineteen (Unigene11539, Unigene12164, Unigene12788, Unigene16018, Unigene17615, Unigene18048, Unigene18594, Unigene18926, Unigene22808, Unigene24462, Unigene24837, Unigene3634, Unigene4753, Unigene6221, Unigene891, and Unigene9670 (K00430), Unigene16856 (K11188), Unigene22599 and Uni-gene25305 (K03782)) out of 497 unigenes (total unigenes) involved in lignin synthesis in the phenylpropanoid biosynthesis pathway (Ko00940) were identified and differentially regulated ( Figure S4). Lignin plays a vital role in keeping the structural integrity of the cell wall and protecting plants from pathogens [54], as well as a main component of wood. Of these 19, different types of peroxidases (Unigene11539, Unigene18926, Unigene4753 and Unigene6221) were upregulated during winter. Recently, a notable remodeling of the transcriptome was reported in Norway Spruce, where monolignol biosynthesis genes showed high expression during the period of secondary cell wall formation as well as the second peak in midwinter. Interestingly, this midwinter peak expression did not trigger lignin deposition [55]. These genes could be preparing for the biosynthesis and distribution of guaiacyl (G), p-hydroxyl phenol (H), and syringyl (S) lignin in developing biomass as soon as the onset of spring. rsity 2021, 13, x FOR PEER REVIEW 12 of 22 Figure 5. GO function analysis of the differentially expressed genes. GO function analysis results for the differentially expressed genes in cambial tissues due to winter and spring seasons into biological processes, cellular components, and molecular functions.

Expression of Lignocellulosic Pathway Genes and Their Validation
Wood, the secondary xylem, is produced from the activity of vascular cambium that is composed of two meristematic initials, fusiform and ray initials [56], with the sequential developmental process including differentiation of vascular cambium cells into secondary xylem mother cells, cell expansion, and massive deposition of secondary walls (where a number of genes involved in vascular tissue differentiation and secondary wall biosynthesis are) [57]. When the wood compression starts, the expression of a number of genes involved in the synthesis of lignocellulosic components (cellulose, hemicellulose, and lignin) and lignans was upregulated in maritime pine [58]. In addition, the onset of wood formation undergoes three periods: winter shrinkage, spring rehydration (32-47 days), and summer transpiration in the stem [59].
In order to explore the roles of cell wall-and hormone-related genes for the seasonal cues, fifteen candidate genes were identified from previous studies (  [71], YUC1 [72,73]), and cytokinin (IPT1 [29] synthesis/pathways). RT-qPCR was employed to study the expression of these wood formation genes in Paulownia during the winter and spring seasons ( Figure 6). Cellulose is synthesized in plant cell walls by large membrane-bound protein complexes proposed to contain several copies of the catalytic subunit of the cellulose synthase, CesA. Here, we found CesA1 and CesA6 were upregulated during spring while CesA3 was moderately downregulated during the winter season. In hybrid aspen, expression analyses of the CesA family showed a specific location in normal wood undergoing xylogenesis, while PttCesA2 seems to be activated on the opposite side of a tension wood [60]. However, in Arabidopsis, the expression levels of the three primary cell wall genes (AtCesA2, AtCesA5, AtCesA6) was increased, but not AtCesA3, AtCesA9 or secondary cell wall AtCesA7 [74]. Our results, along with these studies, indicated that the expression of major primary wall CesA genes accelerate the primary wall CesA complex.
Several proteins encoded by the cellulose synthase-like (CSL) gene family, including CSLA proteins, which synthesize β-(1→4)-linked mannans, and CSLC proteins, which are thought to synthesize the β-(1→4)-linked glucan backbone of xyloglucan are known to be involved in the synthesis of cell-wall polysaccharides [61]. Higher expression of CLSC4 in Paulownia during the spring season indicated that it might involve cellulose synthesis. The fucosyltransferase (FUT1) is an enzyme that transfers an L-fucose sugar from a GDP-fucose (guanosine diphosphate-fucose) donor substrate to an acceptor substrate. The Arabidopsis mur1 (AtFUT1) mutant study [63] exhibited a dwarf growth habit and decreased wall strength indicating the indispensable role of FUT1 function in wood formation. Another key gene family of O-acetyl substituents seems to be very important for various plant tissues and during plant development [75], suggesting an important functional role in the plant. Mutants lacking AXY4 transcript resulted in a complete lack of O-acetyl substituents on xyloglucan in several tissues, except seeds [64]. Biosynthesis of xylan in woody plants is a major pathway for plant biomass. Populus genes PdGATL1.1 and PdGATL1.2, the closest orthologs to the Arabidopsis PARVUS/GATL1 gene, have been shown to be important for xylan synthesis but may also have a role(s) in the synthesis of other wall polymers [65]. The expression of GATL1 homolog was six-fold increased ( Figure 6) in the spring season tissue of Paulownia, implying more xylan biosynthesis. Collapsed xylem phenotypes of Arabidopsis [76] and Physcomitrella patens [66] mutants (irx10) identified mutants deficient in cellulose deposition in the secondary cell wall due to lack of synthesis of the glucuronoxylan. Acetyl transferases are involved in cellulose biosynthesis in plants. In Arabidopsis, the ESKIMO1 (ESK1) gene participated in several functions, and esk1 mutants indicated that ESK1 is necessary for the synthesis of functional xylem vessels in forming secondary cell wall components [67]. Our RNA seq data indicated that all the genes associated with cellulose and hemicellulose synthesis were expressed during the spring season to complete wood formation. (irx10) identified mutants deficient in cellulose deposition in the secondary cell wall due to lack of synthesis of the glucuronoxylan. Acetyl transferases are involved in cellulose biosynthesis in plants. In Arabidopsis, the ESKIMO1 (ESK1) gene participated in several functions, and esk1 mutants indicated that ESK1 is necessary for the synthesis of functional xylem vessels in forming secondary cell wall components [67]. Our RNA seq data indicated that all the genes associated with cellulose and hemicellulose synthesis were expressed during the spring season to complete wood formation. . Expression quantity of the calibrator sample (winter tissue) was set to 1. Data are the mean ± SD. Student's t-test was used to compare significant changes in spring tissues compared to winter tissues. *, p < 0.1; **, p < 0.01; ***, p < 0.001; ns, no significance. . Expression quantity of the calibrator sample (winter tissue) was set to 1. Data are the mean ± SD. Student's t-test was used to compare significant changes in spring tissues compared to winter tissues. *, p < 0.1; **, p < 0.01; ***, p < 0.001; ns, no significance.

Analysis of Hormone-Specific Genes and Their Validation
Plant growth regulator auxins play a key role in regulating wood formation through their action on cambial activity and xylem development [77]. Auxin has been shown to be required for maintaining the cambium in a meristematic state as depleting the auxin in cambium leads to differentiation of cambial cells into axial parenchyma cells [78]. Cytokinins, on the other hand, have a well-established function in cell division during growth and development, and they are called central regulators of cambial activity [79]. The interaction between auxin and cytokinin seems to be essential for the induction of phenylalanine ammonialyase activity in support of lignification [80]. TAA1, which performs the first two reactions in the auxin pathway, is a Trp aminotransferase that converts Trp to IPA in the IPA auxin biosynthesis branch in Arabidopsis [73]. Higher-order mutants in TAA1 showed auxin-related multiple phenotypes. Later, it was identified that the TAA1 gene was essential for hormone crosstalk with ethylene for plant development [71]. Later, new putative functions of IAA production via IPyA and transport were identified [81].
Another group of auxin biosynthesis gene family, YUCCA flavin monooxygenases, control the formation of floral organs and vascular tissues in Arabidopsis [72]. When the TAA family of aminotransferases converts tryptophan to indole-3-pyruvate (IPA), YUCCA (YUC) family participates in converting IPA to IAA [73]. In addition, the authors found that YUC and TAA work genetically in the same pathway and that YUC is downstream of TAA. From our transcriptome and gene expression studies, we observed TAA1 was strongly expressed, but YUC1 expression was negligible during the spring season. Our study identified several transcripts involved in auxin biosynthesis through the tryptophan pathway for cell enlargement and plant growth ( Figure S5).
In Arabidopsis, cambial activity responded to small changes in cytokinin levels indicating that cytokinins are central regulators of cambium activity [79]. Isopentenyltransferase, the rate-limiting step of cytokinin biosynthesis, is an important enzyme playing key roles in meristem maintenance and organ development. Arabidopsis quadruple mutants lacking AtIPT1, AtIPT3, AtIPT5, and AtIPT7 were unable to form cambium and showed reduced thickening of the root and stem, though single mutant atipt3 showed moderately decreased levels of cytokinins without any other recognizable morphological changes. Similarly, increased cytokinin biosynthesis stimulates the cambial cell division rate and increases the production of trunk biomass in transgenic Populus trees [29]. Surprisingly, IPT1 expression was high in winter and moderately reduced in spring. There could be many other members in the IPT family that complements the function of cambial development. Auxin and cytokinin display distinct distribution profiles across the cambium, and elevated cytokinin content leads to an increased cambial auxin concentration [29]. Together, it is very interesting to see the interaction of lignocellulosic pathways genes along with major hormone-regulated genes and their cross-talks to maintain the balance of cambial activities for quality wood formation with alternative seasonal changes ( Figure S6).

Analysis of Simple Sequence Repeats (SSRs)
SSR markers are very useful for multiple applications in plant genetics because of their co-dominance, high level of polymorphism, multi-allelic variance, and abundance, and cross-species transferability [82,83]. In the present study, SSR was identified utilizing the transcriptome of Paulownia cambial tissues because EST-SSR markers have a relatively higher transferability than genomic SSRs [84]. Recent studies showed that abundant EST-SSRs from RNA-seq have agronomic potential and constitute a scientific basis for future studies on the identification, classification, molecular verification, and innovation of germplasms in hawthorn and Lei bamboo [85,86].
We identified 11,338 SSRs from the annotated 61,639 unigenes. We detected 3036 mononucleotides, 5492 dinucleotides, 2493 trinucleotides, 204 tetranucleotides, 194 pentanucleotides, and 344 hexanucleotides motifs (Figure 7). Among the dinucleotide and trinucleotide SSRs, AG/CT repeats represented 2997 SSRs, and AAG/CTT repeats represented 582 SSRs. In mononucleotide, dinucleotide, trinucleotide, quadnucleotide, pentanucleotide, and hexanucleotide repeat categories, the occurrences of repeats were twelve, six, five, five, four, and four, respectively (Table S4). Finally, 6773 oligonucleotide pairs were generated for these identified SSR markers. SSRs and SNPs are the most useful and robust molecular markers for genetics and plant breeding applications [87]. This study provided a set of SSR markers that could be used, for example, in diversity analysis of Paulownia species. In addition, Paulownia tree breeding programs will benefit from the availability of these SSR markers identified from our RNAseq data. Mononucleotide SSRs would be excluded because of the frequent homopolymer errors found in sequencing data and fewer polymorphisms; dinucleotides (46.6%) and trinucleotides (21.2%) contributed most in Paulownia. This is consistent with the EST-SSRs distributions reported in other plant species [88,89]. In plants, SNPs are predominantly beneficial in the construction of high-resolution genetic maps, positional cloning, marker-assisted selection (MAS) of important genes, genome-wide linkage disequilibrium associate analysis, and species origin, relationship, and evolutionary studies [90].
would be excluded because of the frequent homopolymer errors found in sequencing data and fewer polymorphisms; dinucleotides (46.6%) and trinucleotides (21.2%) contributed most in Paulownia. This is consistent with the EST-SSRs distributions reported in other plant species [88,89]. In plants, SNPs are predominantly beneficial in the construction of high-resolution genetic maps, positional cloning, marker-assisted selection (MAS) of important genes, genome-wide linkage disequilibrium associate analysis, and species origin, relationship, and evolutionary studies [90].

Conclusions
Paulownia is a fast-growing, multipurpose timber tree suitable for use as a dedicated lignocellulosic bioenergy crop. To understand the genes involved in the formation of woody biomass related to seasonal cues, a de novo transcriptome study was conducted on vascular cambium tissue from senescent winter vascular cambium tissue and actively growing spring vascular cambium tissue. To the best of our knowledge, this is the first transcriptome-based study on P. elongata, as well as the first transcriptome study performed on Paulownia vascular cambium tissue focusing on seasonal difference. A set of

Conclusions
Paulownia is a fast-growing, multipurpose timber tree suitable for use as a dedicated lignocellulosic bioenergy crop. To understand the genes involved in the formation of woody biomass related to seasonal cues, a de novo transcriptome study was conducted on vascular cambium tissue from senescent winter vascular cambium tissue and actively growing spring vascular cambium tissue. To the best of our knowledge, this is the first transcriptome-based study on P. elongata, as well as the first transcriptome study performed on Paulownia vascular cambium tissue focusing on seasonal difference. A set of transcripts was specifically expressed in the same set of tissues at two different time points. The transcript abundance data confirms the differential pattern of expression of cellulosic, hemicellulosic, lignin biosynthesis specific, and hormone pathway-specific genes. By analyzing the transcriptome from two different temporal treatments (winter and spring), representing two distinct physiological states of the plant, DEGs were identified from both treatments. Cell division is one of the key processes taking place in the cambial zone, and the majority of the cell cycle genes were upregulated during the active stage. The onset of cambial activity began between the end of March and the beginning of April as the increased vacuolization of meristematic cells, and the mitotic activity occurs. However, our current study showed that more genes were downregulated in the spring season. Overall, data from this study will serve as a benchmark for further analysis of molecular mechanisms of wood formation in Paulownia and other trees.