Skip to main content


Transcriptome sequencing and profiling of expressed genes in cambial zone and differentiating xylem of Japanese cedar (Cryptomeria japonica)

The Erratum to this article has been published in BMC Genomics 2016 17:803



Forest trees have ecological and economic importance, and Japanese cedar has highly valued wood attributes. Thus, studies of molecular aspects of wood formation offer practical information that may be used for screening and forward genetics approaches to improving wood quality.


After identifying expressed sequence tags in Japanese cedar tissue undergoing xylogenesis, we designed a custom cDNA microarray to compare expression of highly regulated genes throughout a growing season. This led to identification of candidate genes involved both in wood formation and later cessation of growth and dormancy. Based on homology to orthologous protein groups, the genes were assigned to functional classes. A high proportion of sequences fell into functional classes related to posttranscriptional modification and signal transduction, while transcription factors and genes involved in the metabolism of sugars, cell-wall synthesis and lignification, and cold hardiness were among other classes of genes identified as having a potential role in xylem formation and seasonal wood formation.


We obtained 55,051 unique sequences by next-generation sequencing of a cDNA library prepared from cambial meristem and derivative cells. Previous studies on conifers have identified unique sequences expressed in developing xylem, but this is the first comprehensive study utilizing a collection of expressed sequence tags for expression studies related to xylem formation in Japanese cedar, which belongs to a different lineage than the Pinaceae. Our characterization of these sequences should allow comparative studies of genome evolution and functional genetics of wood species.


Wood represents the main source for terrestrial biomass production and is a major renewable resource for the timber, paper, and bioenergy industries [1]. Genomics approaches have been applied to explore the molecular basis of growth and development in a few forest tree species with economic relevance. Transcript profiling in trees has also specifically focused on wood formation (xylogenesis) because of the ecological significance of forest trees and the economic importance of wood [24]. Wood formation begins from the cambium and generates wood as the end product of secondary vascular system development, which proceeds from cell division to expansion, secondary wall formation, lignification, and finally programmed cell death [5, 6]. Notably, identification of accumulated expressed sequence tags (ESTs) and their expression pattern during wood formation has been achieved in target species for breeding, such as Pinus, Populus and Picea[1, 2, 614].

Japanese cedar (Cryptomeria japonica) is an allogamous coniferous species that relies on wind-mediated pollen and seed dispersal, and it is one of the most important forestry tree species in Japan. The Japanese cedar tree has excellent attributes (straight bole, rapid growth, ease of processing, and pleasant color and scent), and it has been used for house construction, to build wooden ships, barrels, and musical instruments, and for many products intended for daily use for hundreds of years [15]. More than 3,700 Japanese cedar trees have been planted throughout Japan, covering an area of 4.5 million ha and accounting for 44% of Japan’s artificial forests. Seventeen million seedlings are supplied as planting material for forestation every year, making this species very important for Japanese forestry today, as it has been since ancient times [16].

Next-generation sequencing can be a more efficient approach for obtaining functional genomic information. This type of sequencing can result in high transcriptome coverage depth and facilitates the de novo assembly of transcriptomes from species where full genomes do not exist [17, 18]. In addition, by simultaneously measuring the abundance of transcripts for thousands of genes with accumulated sequence information, microarray analysis promises a comprehensive understanding of regulatory gene functions and the growth and development of plants [19]. To understand the molecular mechanisms involved in wood formation and key targets for genetic manipulation and selection of superior wood quality, these techniques will be powerful and efficient tools [20]. The only molecular studies of wood formation in Japanese cedar have identified large numbers of genes that are expressed in male strobili [21]. However, very limited genomics and functional genomics resources related to wood formation are publicly available for Japanese cedar.

The first objective of this paper was to produce an extensive collection of sequenced ESTs found in xylem and cDNA clones to support manufacture of cDNA microarrays and gene discovery efforts in Japanese cedar. The next goal was to elucidate a comprehensive expression profile in the growing season using these microarrays. For this purpose, we identified 55,051 unique sequences by next-generation Roche 454 sequencing using a non-normalized cDNA library from the cambial meristem and its derivatives from Japanese cedar. To gain further insight into seasonal expression patterns, a custom cDNA microarray was designed from the cDNA library obtained and from EST data (inner bark data on ForestGen; [22] and was used to investigate differential gene expression in Japanese cedar during wood formation.

Results and discussion

Microscopic observation of differentiating xylem

Based on anatomical observation of the cambial zone and the differentiating xylem, the tissue underwent seasonal cycles in activity of xylem formation, including cell division, secondary wall formation and lignification, through the growing season (Figure 1A,B). The cambial cells were not active in samples taken on 24 March. An average of only 4.8 cambial cells was found in each radial file, significantly fewer than found in other samples collected in April (p < 0.01), June (p < 0.01) and August (p < 0.05). The expanding tracheids occupied the most space in differentiating xylem in samples taken on 27 April. The number of these tracheids was larger on this date than in samples collected during the growing season. Thus, formation of derived tracheids was most active in samples collected on 27 April. A few secondary wall-forming tracheids were found in some radial files. This indicated that secondary wall formation might have just been reactivated around the day of sampling. The largest number of tracheids at the stage of secondary wall formation and lignification was observed in the 22 June samples. The number was significantly larger than in the other samples collected (p < 0.01). Therefore, the peak activity in xylem formation, including cell differentiation and secondary wall formation, was found in the 22 June samples.

Figure 1

Cross-sections of cambial zone and differentiating xylem of Cryptomeria japonica trees. A. Cross-sections viewed under an ordinary light microscope (24 March) and pairs of ordinary (left) and polarizing (right) light microscope images for the same field (27 April, 22 June, 24 August and 7 October). B. Number of cells in cambial zone and differentiating xylem. Cells in differentiating xylem were categorized into expanding cells, thickening cells and lignifying cells in accordance with ordinary and polarizing light microscope observations.

The number of expanding tracheids in each radial file had significantly decreased from an average of 8.8 cells in samples taken on 22 June to an average of 1.7 cells in samples taken on 24 August (p < 0.01). This indicated that cell division activity in the cambial zone was lower than at earlier stages. Thus, the major activities in differentiating xylem that could be observed microscopically were secondary wall formation and lignification in the samples collected in August and October.

EST sequencing and de novo assembly

Sequencing of cDNA libraries generated a total 308,542 raw reads, with an average length of 405.29 bp. The size distribution of raw reads is shown in Figure 2A, and a summary of sequencing and assembly results is presented in Table 1. After trimming the adaptors and primer sequences, 9,764 sequences were removed due to short length, low complexity, or overall low quality scores. This cleaning and trimming step resulted in 298,778 high-quality reads, corresponding to 96.8% of the original raw sequence. A total of 241,696 high-quality reads was assembled into 11,022 contiguous sequences (contigs over 500 bp), and 40,435 reads were identified as singletons (i.e., reads not assembled into contigs). The size of contigs ranged from 100 to 9,656 bp, with an average length of 1,014 bp. The distribution of contig size is shown in Figure 2B. Contiguous sequences were further assembled into 14,616 isotigs. Isotigs are putative transcripts constructed using the overlapping contig reads provided as input to Newbler cDNA assembler. The size distribution of isotigs ranged from 33 to over 9,656 bp, with an average length of 1,069 bp (Figure 2C). More than 99% of the isotigs were over 100 bp and 50% of the assembled bases were incorporated into isotigs longer than 1,261 bp (N50 = 1,261 bp). The coverage depth for isotigs ranged from 1 to 14, with an average of 1.7 contigs assembled into each isotig (Figure 2D). The isotigs and singletons together resulted in 55,051 unique sequences (Additional file 1: Figure S1).

Figure 2

Summary of EST sequencing and de novo assembly. A) Size distribution and cumulative percentage of raw sequence reads. B) Size distribution and cumulative percentage of contigs. C) Size distribution and cumulative percentage of isotigs. D) Distribution of isotig sequence coverage.

Table 1 C. japonica transcriptome sequencing and assembly summary

Sequence comparison with other species

All unique sequences were searched against the sequences in the National Center for Biotechnology Information (NCBI) non-redundant protein database and The Arabidopsis Information Resource (TAIR) using a BLASTx algorithm E-value of 1e-5 (Figure 3A). A total of 12,606 isotigs (86.2% of all isotigs) and 14,688 singleton sequences (36.3% of all singletons) had significant BLAST matches at NCBI, and 11,958 isotigs (81.8%) and 13,027 singletons (32.2%) had significant BLAST matches at TAIR. When we compared our unique sequences with EST sequences in the Japanese cedar database (ForestGen) and libraries including xylem and cambium tissue from The Gene Index and ForestGen_Xylem (inner bark and sapwood) using a tBLASTx algorithm with an E-value of 1e-5, we found that 25,641 (11,278 isotigs and 14,363 singletons) had significant BLAST matches at ForestGen and 23,524 (11,457 isotigs and 12,067 singletons) to transcripts from Pinus, 24,550 (11,804 isotigs and 12,746 singletons) from Picea, and 15,945 (9,074 isotigs and 6,871 singletons) from Populus and that 13,906 (7,338 isotigs and 6,568 singletons) sequences from ForestGen_Xylem included similar ESTs. The largest overlap was found for the ForestGen database based on lower E-values. A larger overlap was found for other coniferous species than for broadleaf species. Comparison with BLAST results against the ForestGen and ForestGen_Xylem databases indicated that the unique sequences collected in this study were also covered by ESTs previously collected from other organs of Japanese cedar (Figure 3A). When comparing EST sequences in other libraries involving xylem and cambium, though the smallest overlap was found for ForestGen_Xylem, most of the unique sequences overlapped those of three well-known species undergoing xylogenesis (Figure 3B). These results clarified that the previously accumulated Japanese cedar xylogenesis related ESTs were incomplete, whereas our data mostly coincided with ESTs collected from xylem and cambium in other species.

Figure 3

Sequence similarities. A) Number of transcript sequences from C. japonica cambium region similar to sequences in the NCBI, TAIR, ForestGen, ForestGen_Xylem, pine, spruce, and poplar databases according to BLASTx and tBLASTx E-value cutoff values. B) Venn diagram showing the overlap between our collected unique sequences and sequences in four other databases according to a tBLASTx search (E-value <1e-5).

Therefore, our data are expected to be a useful resource for ESTs related to xylem or cambium development in Japanese cedar.

Identifying protein families represented in sequences by Pfam

The unique sequences were investigated for conserved domains using the Pfam database [23] to predict their function. In 55,051 unique sequences, we found that 19,887 (36.1%) of the encoded proteins were similar to members of 4,764 Pfam protein families (E-value < 1e-10). Overall, products of 18,915 (34.4%) of the transcripts from Japanese cedar cambium tissue were similar to members of 4,420 Pfam families when domains of unknown function (DUFs) (317 families) and uncharacterized protein families (UPFs) (27 families) were excluded. The 20 most abundant protein families in cambium tissue of Japanese cedar are shown in Table 2. The frequency of occurrence of members of these families corresponded with previous reports on Japanese cedar male strobili and white spruce [21, 24].

Table 2 Occurrence of the 20 most common Pfam domains in the predicted proteins of unique transcripts from cambium and differentiating xylem of C. japonica

Identifying proteins according to clusters of orthologous groups (COGs) from seven eukaryotic genomes represented in sequences

The unique sequences were searched against the COG database [25] using the BLASTx program. The sequences that showed significant similarity (an E-value < 1e-5) with those in the database were annotated and assigned to designated functional classes. Overall, 22,738 sequences (41.3%) were annotated to known sequences with designated functional classifications, and 3,816 were similar to known genes of unknown function and unassigned sequences in the database (Figure 4). The most frequent functional categories for our data were “Posttranslational modification, protein turnover, chaperones (category symbol O)” and “general function prediction only (R),” which agreed with previous reports [21, 2628]. On the other hand, “signal transduction mechanisms (T)” was the next largest category in the annotated designations of functional classification, including the function unknown (and unassigned) category. This feature differed from previous reports, suggesting that seasonal expression of genes specific to cambial region tissues occurs.

Figure 4

Functional classification and relative levels of ESTs derived from cambium region of C. japonica. Values are shown as percentage of unique transcripts in the pool.

Identification of transcription factors

Transcription factors are proteins that function in controlling the expression of target genes quantitatively, temporally, and spatially [29]. The unique sequences we identified were annotated against the PlnTFDB [30], a recently developed database of transcription factor families for 22 plant species, using the BLASTx program. BLASTx searches revealed 3,085 unique sequences of Japanese cedar with matches against Arabidopsis thaliana and 2,735 unique sequences with matches against Populus trichocarpa with E-values < 1e-5. Of the 82 transcription factor families, these sequences were annotated to 79 in the Arabidopsis genome and 77 in the poplar genome (Tables 3 and 4). The most abundant transcription factor family annotated to A. thaliana was WRKY (WRKY DNA-binding domain), with 218 unique sequences, whereas for P. trichocarpa, it was C3H (Cys4-His-Cys3 zinc finger), with 158 unique sequences. In particular, unique sequences for transcription factors associated with xylogenesis (such as the MYB, NAC and HB transcription factors) were abundant. Therefore, in future studies, it will be necessary to specify which family members are associated with xylogenesis in this species. These transcription factor features were similar to those reported for radiata pine and white spruce [3, 24].

Table 3 Identification of transcripts encoding putative transcription factors in cambium and differentiating xylem in C. japonica against Populus trichocarpa
Table 4 Identification of transcripts encoding putative transcription factors in cambium and differentiating xylem in C. japonica against Arabidopsis thaliana

Comprehensive gene expression changes during xylem formation

In Japanese cedar, physiological and anatomical alterations during cambial activity have been well investigated by anatomical observations [3135], but little is known about the molecular changes that occur. In order to investigate expression of genes in Japanese cedar cambial tissue during xylem formation more comprehensively, we profiled transcripts at five time points during the growing season using microarray analysis. As a result, we identified 10,380 targets that were differentially regulated during xylem formation (p < 0.05, q < 0.2). The differentially expressed genes clustered into 14 different patterns based on their kinetics of gene expression (Figure 5, Additional file 2). These 14 patterns were divided into two expression patterns, associated with upregulation and downregulation, during xylem formation. Overall, we identified 4,019 targets that showed differential expression during the spring reactivation and the peak activity of xylem formation and 6,361 targets that showed differential expression during decreasing cell division and cessation of growth. Cluster A5 was the most abundantly transcribed during spring reactivation and the peak activity of xylem formation. Similarly, the group consisting of abundant clusters B2, 6, and 7 was notably observed during decreasing cell division/cessation of growth. The sequences of all targets in each cluster were also annotated against the COG database (Table 5) using the same settings as the cDNA library. Xylogenesis genes with well-known functions, such as in “Carbohydrate transport and metabolism (G),” “Cell wall/membrane/envelope biogenesis (M),” and “Cytoskeleton (Z),” were abundant categories in cluster A5, which indicated upregulation of these genes during this period. In the major group observed during decreasing cell division/cessation of growth, expressed genes related to tolerance of various conditions and to adjustment of cellular processes, such as “RNA processing and modification (A),” “Signal transduction mechanisms (T),” and “Defense mechanisms (V),” were abundant.

Figure 5

Co-regulation patterns of differentially accumulated transcripts in xylem formation. A total of 10,380 transcripts differentially accumulated in xylem formation were clustered into 14 groups using the Pearson correlation on the Subio platform. The graphs show the average expression profile of each cluster; changes are on a log2 scale. The gene expression pattern is shown as A) upregulation and B) downregulation during xylem formation. The description and expression profile of the individual targets are summarized in Additional file 2.

Table 5 Number of differentially expressed genes according to their cluster and functional COG classification

Cell-cycle related genes

Druart et al [36] reported the expression of 68 homologs of aspen trees based on 80 core cell-cycle genes that were investigated in Arabidopsis[37]. The expression patterns did not correspond to the increasing number of dividing cambial cells during the early phase of cambial cell-cycle activation, leading to the hypothesis of posttranscriptional control of expression after cessation of growth [36]. Similarly, we investigated the expression of Arabidopsis core cell-cycle gene homologs in Japanese cedar. We observed upregulation of 16 of 25 genes from March to April, which suggested that the activation of cell division and induction of cell-cycle genes are correlated in early activity during xylem formation (Figure 6, Additional file 3). Our findings agreed well with data from Arabidopsis[37] and poplar, for which the expression of CDKB and CYCB is regulated seasonally following a rise in temperature [38]. Our data suggested that cambial reactivation occurred in Japanese cedar between 24 March and 27 April based on anatomical observation. To prove the hypothesis of posttranscriptional control of expression after cessation of growth, we therefore harvested samples during this period.

Figure 6

Expression of cell-cycle genes in cambial region during xylem formation. The individual targets are summarized in Additional file 3. All expression data are presented on a log2 scale.

Xylogenesis genes related to phenylpropanoid metabolism

During the development of xylem tissue, primary cell wall biosynthesis, secondary wall deposition, and lignification are important fundamental processes, because of the need for maintaining biological mechanisms conferring adaptability to various environments, compressive strength and defense against pathogens. These processes are also important determinants of wood properties.

The identification and expression profiling of gene family members that are responsible for developmental lignification have been reported for P. trichocarpa, Picea abies, and Pinus teada[3941]. In our study, the expression of the most of these gene family members (Phenylalanine ammonia-lyase (PAL); 4-coumarate: CoA ligase (4CL); Cinnamate-4-hydroxylase (C4H); Hydroxycinnamoyl: CoA shikimate/quinate hydroxycinnamoyl transferase (HCT); p-coumarate-3-hydroxylase (C3H); Caffeoyl-CoA O-methyltransferase; Cinnamyl alcohol dehydrogenase (CCoAOMT); Cinnamoyl-CoA reductase (CCR); and Cinnamyl alcohol dehydrogenase (CAD)) was induced from March to April, and then expression gradually decreased from the peak activity of xylem formation through August (Figure 7A, Additional file 3). The expression of these genes corresponded to our anatomical observation that the number of cambial cells rapidly increased from March to June (Figure 1A,B). These observations indicate that these genes are the main transcripts in developing xylem of Japanese cedar.

Figure 7

Expression of phenylpropanoid metabolism-related genes in cambial region during xylem formation. A) Upregulated genes involved in lignin biosynthesis. B) PAL4 (isotig 10873) and 4CL3 (isotig 11289). The individual targets are summarized in Additional file 3. All expression data are presented on a log2 scale.

Interestingly, enzymes in the early part of the monolignol pathway, acting between PAL and 4CL, are also involved in the biosynthesis of other phenylpropanoids, like flavonoids, coumarins, and stilbene [39]. Lignans, which are monolignol-derived dimers and oligomers involved in such processes as defense reactions, are synthesized through the same pathway [39]. The expression of PAL4 (isotig 10873) and 4CL3 (isotig 11289) was upregulated during dormancy and following cessation of growth (Figure 7B, Additional file 3), which indicates that these enzymes could play roles in defense, such as responses to infection, wounding, drought stress and temperature change.

Lignins result from the oxidative polymerization of p-hydroxycinnamyl alcohols, which can be mediated by both laccase and peroxidase [42]. For the 19 peroxidase superfamily proteins that we examined, the levels of 3 transcripts (isotigs 09523, 13814, and Cj.19051_1) increased during peak xylem formation (Figure 8A, Additional file 3), which corresponded to anatomical observations. The peroxidases that were induced during this period are the strongest candidates for involvement in lignin polymerization. The expression of 7 peroxidase superfamily proteins was upregulated during dormancy and correlated with cessation of growth (Figure 8B, Additional file 3). In P. abies and Pinus sylvestris, high peroxidase activities have also been measured outside the growth period during late autumn, winter, and early spring [43]. Some peroxidase genes generally respond to external stimuli such as wounding, UV-irradiation, bending stress and pathogen infection [39, 44]. These previous reports suggest functions of these genes during the inactive period for the cambium.

Figure 8

Expression of peroxidase superfamily in cambial region during xylem formation. A) Genes upregulated in activity period. B) Genes upregulated during dormancy and cessation of growth. The individual targets are summarized in Additional file 3. All expression data are presented on a log2 scale.

Recently, it was reported that laccase genes Lac4 and Lac17 contribute to constitutive lignification in an A. thaliana mutant [45]. Expression of laccase was 4–5 times higher than peroxidase in developing xylem and young vertical xylem in P. abies[39]. The induction of LAC17 (isotigs 04632, 06065, 08775, and 08985) and LAC4 (isotig 04110) increased rapidly at the peak activity of xylem formation in comparison to other lignification-related genes (Figure 9, Additional file 3). These most highly expressed laccases are also candidates for involvement in lignin polymerization.

Figure 9

Expression of laccase genes in cambial region during xylem formation. The individual targets are summarized in Additional file 3. All expression data are presented on a log2 scale.

Xylogenesis genes related to carbohydrate, cellulose, and hemicellulose metabolism

Cellulose is generally the main component of the plant cell wall, and is synthesized at the plasma membrane by a large multimeric cellulose synthase (CesA) complex [46]. The patterns of expression for almost all expressed Cellulose synthase and Cellulose synthase-like (Csl) genes that were upregulated during peak activity of xylem formation are shown (Figure 10 and Additional file 3). Of these genes, CesA1 (isotig 04782), Ces4/IRX5 (isotigs 14272, 08498), CesA6 (isotig 14123), CesA7/IRX3 (isotigs 04866, 09868), and CesA8/IRX1 (isotigs 02784, 14052) were rapidly induced from March to April, and then their expression gradually decreased from the period of peak xylem formation until cessation of growth (Figure 10, Additional file 3), which corresponded with anatomical observations. Interestingly, Ces4/IRX5, CesA7/IRX3 and CesA8/IRX1 are required for cell-wall biosynthesis in vascular tissue of Arabidopsis and rice [47, 48]. Similarly, orthologs identical to these three genes are involved in secondary cell-wall biosynthesis in developing xylem of wood species such as Populus and Pinus[14, 4953]. These findings suggest that the functional roles of these orthologs are conserved in cell-wall synthesis of vascular tissue in herbaceous and woody dicotyldeons, monocotyledons and gymnosperms [14]. Additionally, a membrane-bound endoglucanase, KORRIGAN1 (KOR1), and a glycosylphosphatidylinositol-anchored protein, COBRA (COB), have been implicated in cellulose biosynthesis in Arabidopsis[54, 55]. Both the sequence of the orthologs and their functional roles are reportedly conserved in Populus and Picea[56, 57]. The Japanese cedar homologs of KOR1 and COB were upregulated during peak xylem formation, suggesting a conserved functional role (Figure 11A,B, Additional file 3).

Figure 10

Expression of cellulose synthase and cellulose synthase-like genes in cambial region during xylem formation. The individual targets are summarized in Additional file 3. All expression data are presented on a log2 scale.

Figure 11

Expression of KORRIGAN1 and COBRA in cambial region during xylem formation. A) KORRIGAN1 (KOR1). B) COBRA (COB). The individual targets are summarized in Additional file 3. All expression data are presented on a log2 scale.

Hemicelluloses, including glucomannans and xyloglucans, are major components of the plant secondary cell wall. Most genes related to the glucomannan/galactomannan pathway (such as GDP-D-mannose 4,6-dehydratase1, mannose-1-phosphate guanylyltransferase, phosphomannomutase, mannosyltransferase family protein, galactosyltransferase family protein) and to synthesis of xylan (such as UDP-glucuronic acid decarboxylase: UDP-xylose synthase, β-(D)-xylosidase) and xyloglucan (such as xyloglucan endotransglycosylase, xyloglucan endotransglycosylase/hydrolase) were induced during xylem formation, which suggests that the encoded proteins play an active role in secondary wall formation (Figure 12, Additional file 3). Recently, a number of genes encoding putative glycosyltransferases required for xylan synthesis or deposition have been identified in Arabidopsis using knockout mutants [5863]. The expression of most of these genes (IRX7/FRA8, IRX9, IRX10-like, IRX14, and IRX15) increased during xylem formation, indicating conserved functional roles of these orthologs (Figure 13A, Additional file 3). Xyloglucan is incorporated and modified in the cell-wall network by xyloglucan endotransglycosylases and hydrolases (XTHs, also known as XET/hydrolases and XEHs) [64, 65]. We observed 10 genes involved in cell-wall biosynthesis to be upregulated at peak xylem formation; however, the expression of 12 genes was downregulated (Figure 13B, Additional file 3). Some XTH genes are induced in dormant cambium and cold-stressed organs [65]. In Japanese cedar, these 12 genes may be candidates for this functional role.

Figure 12

Expression of hemicellulose-related gene family in cambial region during xylem formation. The individual targets are summarized in Additional file 3. All expression data are presented on a log2 scale.

Figure 13

Expression of glycosyltransferase and xyloglucan endotransglycosylase in cambial region during xylem formation. A) Glycosyltransferase. B) Xyloglucan endotransglycosylase. The individual targets are summarized in Additional file 3. All expression data are presented on a log2 scale.

Sucrose synthase (Sus) catalyzes formation of UDP-glucose, the immediate substrate for cellulose biosynthesis. Members of the Sus gene family in many plant species are divergent in function and differentially expressed during plant development [66]. In Pinus and Populus, some Sus genes showed an expression pattern identical to that of Ces genes in developing xylem [5, 14]. Similarly, the expression of Sus (isotig 12351) was upregulated from April to June, which suggested it as a robust candidate gene for involvement in xylem formation (Figure 14A, Additional file 3). Interestingly, all other Sus genes were rapidly downregulated from March to June, and then expression gradually increased through October. Sus gene activity is considered to be associated with environmental stresses, such as cold, drought and O2 deficiency [36, 67]. We observed upregulation of these genes in March and October, months showing markedly low temperature, which suggested synthesis of cryoprotectants and responses to cold stress (Figure 14A, Additional file 3). Additionally, because reactivation of the cambium in the spring occurs before any significant photosynthesis activity, the induction of the Sus gene and various invertases indicates that sucrose catabolism generates hexose that can be metabolized via glycolysis during this period [36]. The reaction catalyzed by sucrose phosphate synthase (SPS) plays an important regulatory role in controlling Sus genes in plants [68]. In hybrid poplar, an AtSPS transgenic hybrid has altered phenology, such as timing of leaf senescence and bud break, compared to wild type [69]. Some invertase and SPS genes also showed an expression pattern identical to that of the Sus genes (Figure 14B,C, Additional file 3). In Japanese cedar, these expression profiles could indicate that these genes are involved in providing an alternative source of energy and carbon skeletons in the early period of cambial reactivation in spring.

Figure 14

Expression of sucrose synthase, invertase, and sucrose phosphate synthase in cambial region during xylem formation. A) Sucrose synthase. B) Invertase. C) Sucrose phosphate synthase. The individual targets are summarized in Additional file 3. All expression data are presented on a log2 scale.

Transcription factors

Several transcription factor family members, such as NAC, MYB, zinc finger proteins, and proteins with a Lim domain or a homeodomain are thought to help regulate secondary cell-wall biosynthesis [7076]. In particular, some transcription factors of the NAC and MYB subfamilies are master switches in the transcriptional network for secondary cell-wall biosynthesis [75]. In Arabidopsis, some MYB genes (AtMYB20, 43, 46, 52, 63, 83, 85, 99, 103, and 118) are significantly upregulated in expression just when xylem vessel elements actively form [77]. In conifers, Picea (PgMYB2, 4, 8) and Pinus (PtMYB1, 4) MYB genes are also involved as transcriptional regulators in lignin metabolism and/or wood formation in stem and root [7880]. Most of these MYB genes are clustered in a phylogenetic tree of MYBs from spruce, Pinus, Arabidopsis, and the nearest sequences from other species [78]. We found that 34 MYB genes were upregulated during the peak activity of xylem formation (Figure 15A, Additional file 3). AtMyb20 (isotig 05701) and AtMyb43 (Cj.5920_1) were expressed preferentially; however, AtMyb103 (isotig 03851) was downregulated during this period. These findings suggest that the functional roles of MYB20 and MYB43 orthologs are conserved in cell-wall synthesis of vascular tissue in Japanese cedar and other species. Recent molecular and genetic studies have revealed that a subgroup of Arabidopsis NAC domain transcription factors (SND1, NST1, VND6, 7) are master switches regulating a cascade of downstream transcription factors, leading to activation of secondary wall biosynthesis [8186]. Although among these NAC domains, only a VND6 homolog (Cj17576_1) was included on our array; its expression was moderately decreased, as with 9 NAC family homologs (anac2, 8, 28, and 45) during peak xylem formation (Figure 15B, Additional file 3). The homologs of other cell-wall biosynthesis-related transcription factors (LIM, HB, b-ZIP, WRKY) were induced during this time (Figure 15C, Additional file 3), implying that these genes could be important in regulating downstream genes.

Figure 15

Expression of cell wall-related transcription factors in cambial region during xylem formation. A) MYB. B) NAC. C) Other cell wall-related transcription factors (LIM, HB, b-ZIP, WRKY). The individual targets are summarized in Additional file 3. All expression data are presented on a log2 scale.

Hormonal regulation of the activity-dormancy cycle

Auxin has been implicated as a key signal regulating cambial cell proliferation and cambial meristem identity [12]. In the cambial region, the amount of IAA varies seasonally, and rapid induction in cambial activity occurs in spring to early summer [8789]. The IAA distribution shows a radial gradient and is most concentrated in the cambial region [90, 91]. The positive correlation observed between the regions with high IAA levels and the number of cells in the same region suggests that the gradient in endogenous IAA level controls the number of cambial cells [92]. The expression of some auxin signaling and transport component genes (Aux1, IAA16, 27, Auxin efflux carrier (PIN1, 2), Auxin response factor (ARF1, 2, 4), and SAUR-like auxin-responsive protein) was upregulated in April, and then gradually decreased through October (Figure 16A, Additional file 3). The concentration of IAA in P. sylvestris is high at the start of cambial reactivation, declines when the number of differentiating tracheids begins to increase, and then rises as the number of cells decreases [88]. In Japanese cedar, our results indicate that these genes are regulated early in xylem formation. Auxin signaling is mediated through the ubiquitin-proteasome pathway, in which AUX/IAA proteins are degraded through SCFTIR1 complexes, composed of cullin, SKP1, F-box protein, and RBX1 protein [93, 94]. We observed inverse expression of these genes relative to the expression of some auxin signaling and transport component genes, suggesting that auxin signaling and transport component genes are repressed during cessation of growth.

Figure 16

Expression of hormonal regulation-related genes in cambial region during xylem formation. A) Auxin signaling and transport component. B) GA biosynthesis-related and signaling genes. C) ABA biosynthesis-related and signaling genes. The individual targets are summarized in Additional file 3. All expression data are presented on a log2 scale.

Gibberellins (GAs) act synergistically with auxin in stimulating cambial growth [95]. The analysis of transgenic aspen indicated that GAs are required both in xylogenesis, which is likely mediated by GA signaling in the cambium, and in fiber elongation in the developing xylem [96]. In angiosperm trees, application of GA results in the formation of wood fibers with enhanced thickness of the inner layers of cell walls [97, 98]. We found that a homolog of GA3-oxidase (GA3ox, Cj.17342_1), implicated in the last step of GA biosynthesis, and the receptor gene GID1 (Cj.5192_1, isotig13598) were moderately upregulated at peak xylem formation (Figure 16B, Additional file 3). The genes encoding GA biosynthetic enzymes GA20-oxidase (GA20ox) and GA3ox are particularly important for control of bioactive GA levels [99]. GA signaling operates as a derepressible system that is moderated by DELLA-domain proteins, which are transcriptional regulators that repress GA responses [99]. Like DELLA-domain proteins, the homologs of RGA (Cj 552_1, 1674_1) were expressed inversely to these genes (Figure 16B, Additional file 3). These findings suggest that genes involved in GA signaling have an important role in xylem formation.

Abscisic acid (ABA) content increases during abiotic stress, and especially protects plant water status. In poplar cambium, ABA levels are increased by short days and by short days with low temperature in late autumn and during cambial reactivation in early spring [36]. Genes related to ABA biosynthesis and signaling, such as ABA4, NCED, CYP707A, PP2C (HAI1, 2), SnRK2.6 (OST1), ABRE (ABI1, 5, ABF1), and PYL (1,4,10), were upregulated in March and October (Figure 16C, Additional file 3). Most of these genes were rapidly downregulated from March to April, suggesting that their downregulation is coincident with release from cold hardiness and the improvement in water deficit on cambial reactivation. The Arabidopsis CYP707A gene family (CYP707A1, 2), involved in ABA catabolism, controls seed dormancy [100]. Therefore, our observations suggest that ABA is degraded during cambial reactivation in Japanese cedar. In the apex of hybrid aspen, some 9-cis-epoxycarotenoid genes (such as NCED), which are involved in ABA biosynthesis, are induced after 5 weeks of short-day treatment, which also induces growth cessation [101]. As seen in our data (Figure 16C, Additional file 3), these genes (Cj13501_1, 2567_1, 8387_1) were upregulated in accordance with changes in day length. Other ABA biosynthesis-related and signaling genes were also upregulated from August to October, indicating they may be induced in response to several abiotic stresses (such as cold and drought) that also lead to cessation of growth (Figure 16C, Additional file 3).

Development of cold hardiness in activity-dormancy cycle

On 24 March, before cambial reactivation, cold hardiness was maintained in cambial cells (Figure 1A,B). On the other hand, the number of expanding cells and cells depositing secondary walls, as well as temperature and day length, rapidly decreased from 24 August to 7 October (Figure 1A,B, Additional file 4: Figure S2), suggesting acquisition of cold hardiness on growth cessation. The transcriptional regulators and modulating genes involved in the acquisition of cold hardiness of Japanese cedar have not been identified. In Arabidopsis, the ICE1 (inducer of CBF expression) and CBF (C-repeat binding factor) family transcription factors are respectively upstream and downstream regulators of the cold-responsive transcriptome and freezing tolerance [102, 103]. The role of the CBF family as transcriptional activators in cold acclimation of Arabidopsis has been maintained in Populus; in particular, CFB1 and CFB3 show significant induction in Populus stems [104]. The only ICE1 homologs presented on our array were induced in conditions consistent with maintaining and acquiring cold hardiness (isotigs 05865, 14021) (Figure 17, Additional file 3), so these genes are candidates for this functional role. Druart et al. [36] listed and clustered the expression pattern of genes involved in cold hardiness from three data sets: 1. poplar genes induced by low temperature and atCBF overexpression; 2. a poplar homolog of an Arabidopsis gene induced by low temperature; and 3. poplar homologs of genes involved in the development of cold hardiness in three other tree species [36]. We clustered the expression patterns of these genes based on their timing of induction. Clustering of the 138 homologs of these genes yielded three main groups, two associated with autumn transition and early spring (clusters 1, 2), and another with spring reactivation (cluster 3) (Figure 18, Additional file 3). Their expression rapidly decreased from March to April, which suggested repression and release of cold hardiness in cluster 1. In cluster 2, expression (for example, of Cold regulated gene and dehydrin family protein) was up- or downregulated in accordance with changes in day length. In acquisition of cold hardiness, these genes were moderately upregulated prior to reduction in temperature (Additional file 4: Figure S2). This finding implies that a signal other than low temperature (such as short days) must trigger the induction of these genes in autumn under natural conditions [36, 105]. In cluster 3, some homologs were superinduced from March to April, which corresponded to previous findings [36]. The exact role of this superinduction is unclear; however, it might reflect a need to protect the very sensitive dividing cambial cells from sudden drops in temperature during early spring [36].

Figure 17

Expression of ICE1 in cambial region during xylem formation. The individual targets are summarized in Additional file 3. All expression data are presented on a log2 scale.

Figure 18

Expression of low temperature-induced genes in cambial region during xylem formation. Low temperature-induced and cold hardiness-related genes were clustered using the Pearson correlation on the Subio platform into three main patterns of expression during the autumn transition and early spring (clusters 1, 2) and during spring reactivation (cluster 3). The individual targets are summarized in Additional file 3. All expression data are presented on a log2 scale.

Because of limited photosynthesis in autumn, plants must derive the energy and carbon required for the acquisition of cold hardiness from some other source [36]. Conversion of starch to sugar is a key metabolic process associated with the entry into dormancy, as starch-derived sugars serve several purposes, for example as cryoprotectants as well as a source of energy [12]. Transcriptional induction of a key enzyme of the starch breakdown pathway occurs in poplar cambium in autumn and during dormancy [12, 36]. Our data showed that most homologs involved in the process, such as β-amylase and phosphoglucan-water dikinase, were induced in March and October (Figure 19, Additional file 3). These profiles suggested them as candidate genes for regulating the availability of an alternative energy and carbon source during limitations on photosynthetic activity.

Figure 19

Expression of starch-breakdown related genes in cambial region during xylem formation. The individual targets are summarized in Additional file 3. All expression data are presented on a log2 scale.

Cytoskeleton-related genes

The cytoskeleton regulates cellular polarity, morphology, and movement through its involvement in cellular events such as vesicular traffic, organellar movement, abiotic and biotic stress sensing, signal transduction, and cell wall biosynthesis [106]. Plant cell morphogenesis relies on the organization and function of two polymer arrays separated by the plasma membrane: the cortical microtubule cytoskeleton and cellulose microfibrils in the cell wall [107]. In the large S2 layer of secondary fiber cell walls, the orientation of microfibril deposition, which is directed by cortical microtubules, is an important trait determining wood quality and wood stiffness or elasticity, and is referred to in trees as the microfibril angle [108]. The expression of genes for α- and β-tubulins, which comprise dynamic arrays of cortical microtubules, appear to play a role in determining these characteristics during xylem development in Populus and Eucalyptus[108, 109]. Our data showed that most tubulin gene family members are highly expressed during peak xylem formation (Figure 20A, Additional file 3). Similarly, most other cytoskeleton-related genes such as Actin and genes encoding actin-related or -interacting proteins (actin binding protein, actin related protein, actin depolymerization protein, villin, fibrin profiling, capping protein), microtubule-associated protein (MAP), microtubule-motor family protein, microtubule end binding protein, and kinesin were also induced (Figure 20B,C,D, Additional file 3).

Figure 20

Expression of cytoskeleton related genes in cambial region during xylem formation. A) α/β-Tubulin. B) Actin, actin-related and actin interacting proteins. C) Kinesin gene family. D) Microtubule-associated protein gene family. The individual targets are summarized in Additional file 3. All expression data are presented on a log2 scale.

Actin forms microfilament structures by self-polymerization and interactions with numerous actin-binding proteins. In our data, the four homologs of atACT7 (isotigs 05994, 09744, 11932, and 14133) were upregulated during peak xylem formation, along with actin-related/interacting protein and a gene encoding a kinesin family protein (Figure 20B,C, Additional file 3). The homologs clustered closely in a group with a Pinus homolog in a phylogenetic tree of actin from Arabidopsis and the nearest sequences from other species [110]. atACT7 is preferentially expressed in younger, rapidly developing tissue, such as during germination and root growth in Arabidopsis[110, 111]. These findings correspond with our findings in developing xylem.

The organization and dynamics of microtubules are regulated by MAPs[108]. Our study found a gene encoding a MAP (MAP65-1: isotig 05735, 09873) that was more strongly transcribed than other MAP genes (Figure 20D, Additional file 3). AtMAP65-1 is able to promote tubulin polymerization, enhance microtubule nucleation, and decrease the critical concentration for tubulin polymerization [112]; this role agrees reasonably well with what would be expected from our expression pattern and anatomical observations.

Validation of microarray expression of 12 selected genes by qRT-PCR

Microarray expression data of 12 differentially transcribed genes selected in this study were validated by qRT-PCR using the same RNA samples used for the microarray experiments. Transcript accumulation measured by qRT-PCR was fairly consistent with the microarray results for all 12 validated genes (Figure 21), particularly in the ranking of magnitude of expression, indicating that the microarray experiments in this study were sufficiently reliable for the identification of genes that may influence xylem formation in the cambium of Japanese cedar.

Figure 21

Validation of microarray expression of 12 selected genes by qRT-PCR. A total of 12 genes were selected in the validation using qRT-PCR: (A) Phenylalanine ammonia-lyase (PAL: isotig 10873), (B) Cinnamate-4-hydroxylase (C4H: isotig 09462), (C) 4-Coumarate:CoA ligase (4CL: isotig 04988), (D) Hydroxycinnamoyl-CoA shikimate/quinate hydroxycinnamoyl transferase (HCT: isotig 10178), (E) p-Coumarate-3-hydroxylase (C3H: isotig 02271), (F) Caffeoyl-CoA O-methyltransferase; Cinnamyl alcohol dehydrogenase (CCoAOMT: isotig 10196), (G) Cinnamoyl-CoA reductase (CCR: isotig 05079), (H) Cinnamyl alcohol dehydrogenase (CAD: isotig 02638), (I) Sucrose synthase (Sus: isotig 12351), (J) Cellulose synthase (Ces: isotig 08498), (K) α-Tubulin (Tub α: isotig 02753), (L) β−Tubulin (Tub β: isotig 13384). In qRT-PCR, Ubiquitin (UBQ) was used as a reference gene, and the data were calibrated relative to transcript levels in the 24 March sample. Error bars show standard deviation for three replicates.


In this study, we obtained 55,051 unique sequences by sequencing a non-normalized cDNA library from the cambial meristem and derivative cells of Japanese cedar. A custom cDNA microarray was designed based on this library and EST data to investigate seasonal gene expression in Japanese cedar. This is the first comprehensive study of an extensive collection of EST sequences and expression studies related to xylem formation in Japanese cedar. Because Japanese cedar belongs to a different lineage than the Pinaceae, comparison of data could lead to significant findings on genome evolution in coniferous species. Our data may also be a useful resource for forward genetics and functional genetics studies in wood species.


Plant material

Tissue from the cambium region (including phloem and the differentiating xylem) was taken from four 15-year-old trees of Cryptomeria japonica plus-trees, clones Chousui8, Iiyama9, Nisihkawa10 and Tano1, in Hitachi, Ibaraki Prefecture for molecular analysis. The daily minimum and maximum temperatures were also recorded during the study (Additional file 4: Figure S2). The harvested tissues were immediately frozen in liquid nitrogen in the field, and then stored in the laboratory at −80°C for later RNA extraction. A square block (approximately 1 cm2) was collected for microscopy and fixed in FAA (formalin: acetic acid: 50% alcohol, 5:5:90) in the field. To evaluate how gene expression and morphological development in the cambial region changed over a single growing season, tissues from this region were collected from different trees at the same time (around 10 AM) on 15 different dates from 2010 to 2011: on 9 March, 9 April, 10 May, 1 June, 24 June, 16 July, 16 August, 19 September, 29 September, and 29 October in 2010 for construction of cDNA libraries and on 24 March, 27 April, 22 June, 22 August, and 7 October in 2011 for cDNA microarrays and for anatomical observation.

Anatomical observation of the cambial zone and the differentiating xylem

Small blocks were collected from stems corresponding to those used for microarray analysis. Thin sections were prepared from embedded tissue in blocks of LR White Resin (London Resin Co., Basingstoke, UK) and stained with safranin and Alcian blue 8GX. Anatomical observations were carried out under both an ordinary light microscope and a polarizing light microscope. The number of cells in the cambial zone and the number of expanding tracheids, secondary wall forming tracheids and lignified tracheids in each radial file were counted under the microscope. The number of cells at each growth stage was statistically compared by a Student’s t-test between samples.

RNA extraction and pyrosequencing

Total RNA was isolated from tissue of the cambium region and differentiating xylem of plus-trees using an RNeasy Plant Mini kit (Qiagen, Gaithersburg, MD, USA) for Chousui8 samples from ten different dates. The quality of total RNA was assessed by measuring the ratio of absorption at 260 nm and 280 nm via an Agilent Bioanalyzer 2100 (Agilent Technologies, Palo Alto, CA, USA). cDNA synthesis from a mixture of ten RNA samples, nebulization, adaptor ligation, emulsion PCR and sequencing were done at Hokkaido System Science Co., Ltd. (Sapporo, Hokkaido, Japan). Sequencing was performed using a Roche 454 Genome Sequencer platform (Roche/454 Life Sciences, Branford, CT, USA) with FLX or Titanium technology.

Assembly of ESTs from sequences obtained on the 454 platform

Using GS FLX pyrosequencing software, we selected high-quality sequences (> 99.0% accuracy on single base reads) for further processing and assembly. Trimmed and cleaned sequences were assembled using the cDNA assembly feature of Roche Newbler software v. 2.3 (Roche/454 Life Sciences). To obtain clean ESTs, adapter trimming and poly(A/T) removal were performed by the cutadapt tool [113], then short sequences (< 50 bp) were removed and the remaining sequences evaluated using the BLASTN algorithm against C. japonica microsatellite sequences obtained from NCBI ( [114], and Arabidopsis thaliana retrotransposon sequences obtained from TAIR ( [115]; reads with alignment length of 20 nt or more and percent identity of 90% or more were considered “hit reads” against these sequences. De novo assembly was performed using GS De Novo Assembler v2.3 (provided with the Roche GS FLX sequencer) with default parameters (minimum overlap length of 40, minimum percent identity of 90).

Functional annotation with the BLAST program

The assembled unique sequences putatively encoding proteins were searched against the Arabidopsis protein database in TAIR [115] and the NCBI non-redundant database [114] using the BLASTx algorithm. In addition, our transcripts were also searched against the ForestGEN database [22] using the tBLASTx algorithm. A typical cutoff E-value < 1e-5 was used. To identify known protein families, the unique sequences were also searched for the presence of Pfam domain sequences (release 21.0) using the blastx algorithm (E-value < 1e-10) [23]. Similarities to ESTs from libraries derived form xylem and/or cambium of Pinus, Picea, Poplar and Japanese cedar were determined with the tBLASTx program. We used Pinus Gene Index release 6.0 (PGI_libraries PHJ, PHM, ONA, PJD, ERF, 2NV, CJQ, 11 F, 0TU, PJQ, M7S, M7N, 9UQ, 72B, 0U0, 0TT-2, M7Q, M7R, M7P, M7O, PJT, PJM, ERE, CER, PHN, NIL, ERD, BTR, CCS, CJP, 8FB, PJR, ONB, OI1, CJS, ERB, 9NV, 5BN, 1RR, and 0TV), Spruce Gene Index release 2.0 (Sgi_libraries KH2, H5M, H5L, FKG, KGV, EOT, PHL, EOR, KH1, KH0, FH7, FKM, F7N, LCC, F7O, F7U, LCF, IQE, EOS, LCD, LCN, LCM, IQG, IQD, FH9, F7V, IQF, FKL, EOQ, and LCE), Poplar Gene Index release 3.0 (PplGI_libraries EA1, 9BN, BMF, G26, NIQ, EA2, ASV, NL3, FKA, DRG, F8V, F8D, 1CV, LRS, G22, G21, DRC, LRR, DRF, G27, and DQP) from The Gene Index Project website ( [116], and the ForestGen database (inner bark and sapwood data) as EST databases [22].

PlnTFDB [30], a recently developed database of transcription factor families for 22 plant species, was used to identify putative transcription factors expressed during Japanese cedar wood formation. Blastx searches were performed on matches against A. thaliana and P. trichocarpa in the PlantTFDB with E-values < 1e-5.

The unique sequences were searched locally against a database of clusters of orthologous groups (COGs) from seven eukaryotic genomes [25]. The COGs are comprised of three databases containing orthologous proteins from at least three out of seven eukaryotic species (KOGs), proteins from two species (TWOGs), and lineage-specific expansion groups (LSEs). Sequences with E-values < 1e-5 were considered to have significant homology, and were classified following the KOG functional classification.

The sequences of ESTs have been submitted to the DNA Data Bank of Japan under accession numbers DC882454 through DC883482.

Microarray analysis

We built a custom microarray platform containing 60-mer oligonucleotide probes designed based on 14,612 isotigs (probes to 4 isotigs could be not designed) from all isotigsin proprietary NGS data and 3,470 EST sequences from the “sapwood” and “inner bark” categories (including a full-length cDNA library) in the ForestGen database [22]. A set of 18,082 probes was selected and accommodated in the NimbleGen 4 × 72 K array format (Roche-Nimblegen Inc., Waldkraiburg, Germany), which can examine the expression levels of up to 20,000 genes for four samples at the same time. Therefore, in this format, 18,082 probes were accommodated at least in triplicate in our custom array. For microarray analysis of five sampling dates, we used four biological replicates and three technical replicates for each sample (Additional file 1: Figure S1). Total RNA was extracted with a Plant RNeasy Mini Kit, and DNase was treated in-column with an RNase-Free DNase set (Qiagen). The A260/A280 ratios of RNA samples used for hybridization ranged from 1.7 to 2.0. An Agilent 2100 Bioanalyzer analyzed the integrity of RNA samples. RNA integrity values of samples used for hybridization ranged from 8.1 to 10.0. Double-stranded cDNA was synthesized using a SuperScript double-stranded cDNA synthesis kit (Invitrogen, Carlsbad, CA, USA) with random 6-mers following the manufacturer’s protocol. Cy3 labeling and hybridization were performed by NimbleGen using standard procedures. Labeled and hybridized slides were scanned using a NimbleGen MS 200 microarray scanner to generate paired files.

Because there were three or four spots for each target, the paired files contained redundant signal intensities for all probes. We took medians as representing intensities to avoid the effect of outliers, and loaded them into Subio platform software (Subio Inc., [117]. Intensity values were normalized at the 75th percentile, and then transformed into log2 ratios based on the average of the 60 samples, which were composed of 5 time points with 12 replicates each. The data presented in this study have been deposited in NCBI’s Gene Expression Omnibus and are accessible through GEO Series access number GSE53034.

Of the total 18,082 target genes, 748 with raw signal intensities not exceeding 1,000 in any samples were filtered out. We calculated the averages of log2 ratios at each time point, and excluded an additional 6,273 genes with expression levels hardly varying over time (between −0.5 and 0.5). We tested the 11,061 genes by ANOVA (p < 0.05 and BH-FDR < 0.2) to extract 10,380 genes with expression levels that varied for at least one time point. Hierarchical clustering (unweighted pair group method with arithmetic mean, Pearson correlation) was used to identify groups of co-expressed genes. We extracted clusters from tree nodes (Figure 5). We additionally created trees with gene sets manually selected based on biological knowledge.

Validation of quantitative RT-PCR

Independent verification of microarray results was carried out by qRT-PCR analysis using total RNA from the cambium region tissues used for microarray experiments. Total RNA (500 ng) was reverse-transcribed using the PrimeScript II 1st strand cDNA synthesis kit (Takara Bio, Otsu, Shiga, Japan) with random 6-mers following the manufacturer’s instructions. The resulting first-strand cDNA was diluted 1:5 in water before real-time PCR. Primers were designed using Primer Express software ver. 3.0 (Applied Biosystems, Foster City, CA, USA), with a melting temperature (Tm) between 60 and 65°C, and produced amplicons between 100 and 250 bp. Specific primer pairs were designed for each gene: Phenylalanine ammonia-lyase (PAL) (isotig 10873: forward 5′-GACCCAGGACGGGAAAGAG-3′, reverse 5′-TAGGCTGGAGTTCAAACGGTTT-3′); 4-Coumarate:CoA ligase (4CL) (isotig 04988: forward 5′-CAGTCGTCGCCAACTATGACA-3′, reverse 5′-ACGGCATCTTCCAGGTCCTT-3′) Cinnamate-4-hydroxylase (C4H) (isotig 09462: forward 5′-CGTTGAGAAGCTGCCGTATCT-3′, reverse 5′-CGTCAAGGGAGGCTTCTTCA-3′); Hydroxycinnamoyl-CoA shikimate/quinate hydroxycinnamoyl transferase (HCT) (isotig 10178: forward 5′-GCCCATCCATGATGCAGATT-3′, reverse 5′-GACTGGGCAAAATGAAACCAA-3′); p-Coumarate-3-hydroxylase (C3H) (isotig 02271: forward 5′-TCACATGGACCCCTCCTGAA-3′, reverse 5′-CGGTAGAGATGCTCAGGCAAT-3′); Caffeoyl-CoA O-methyltransferase; Cinnamyl alcohol dehydrogenase (CCoAOMT) (isotig 10196: forward 5′-ACTGCAGAGGCTTCCAAGGA-3′, reverse 5′-TCGCTCTGAAGGAGACTCTTGTG-3′); Cinnamoyl-CoA reductase (CCR) (isotig 05079: forward 5′-CAGGAGCGGGAGGATTTATTG-3′, reverse 5′-CCTCTGGATTGCGAACTGTTC-3′); Cinnamyl alcohol dehydrogenase (CAD)(isotig 02638: forward 5′-GCAGAGGCAGGCAAGAGATG-3′, reverse 5′-AGTCACATGATGCCCAAATGC-3′); Cellulose synthase (Ces) (isotig 08498: forward 5′-CATGGCCTGGGAACAACACT-3′, reverse 5′-ATGCGAGGCAGTTCGTTACC-3′); Sucrose synthase (Sus) (isotig 12351: forward 5′-ACGACTGTTCTTGGCAAACCAT-3′, reverse 5′-ATTGAGCGACCGGAACAAAC-3′); α-Tubulin (Tub α) (isotig 02753: forward 5′-CATCCTTGGGCACAACATCTC-3′, reverse 5′-TGCCTTTGAGCCTTCTTCCAT-3′); βTubulin (Tub β) (isotig 13384: forward 5′-TACACTGGTGAGGGCATGGA-3′, reverse 5′-GCATCCTCATCCGCAGTTG-3′); and the endogenous control Ubiquitin (UBQ) (forward 5′-CGTTAAAGCCAAGATCCAGGACAA-3′, reverse 5′-TCCATCCTCAAGCTGTTTCCCA-3′). For each sample, triplicate quantitative PCR assays were performed using Power SYBR Green PCR master mix (Applied Biosystems) with ROX reference dye according to the manufacturer’s protocol. Amplification was carried out with a StepOnePlus system (Applied Biosystems). After an initial 10-min activation step at 95°C, 40 cycles (95°C for 15 s and 60°C for 1 min) were performed, and a single fluorescent reading was obtained after each cycle immediately following the annealing/elongation step at 60°C. Preliminary quantitative PCR assays were performed to evaluate primer pair efficiency and absence of genomic DNA contamination using a negative control. A melting curve analysis was performed at the end of cycling to ensure amplification of a single product. For relative quantification and comparisons, we used the delta-delta-Ct method with Ubiquitin as the normalization internal control gene.


  1. 1.

    Li X, Wu HX, Southerton SG: Seasonal reorganization of the xylem transcriptome at different tree ages reveals novel insight into wood formation in Pinus radiata . New Phytol 2010, 187:764–776.

  2. 2.

    Pavy N, Boyle B, Nelson C, Paule C, Giguère I, Caron S, Parsons LS, Dallaire N, Bedon F, Bérubé H, Cooke J, Mackay J: Identification of conserved core xylem gene sets: conifer cDNA microarray development, transcript profiling and computational analyses. New Phytol 2008, 180:766–786.

  3. 3.

    Li X, Wu HX, Dillon SK, Southerton SG: Generation and analysis of expressed sequence tags from six developing xylem libraries in Pinus radiata D. Don. BMC Genomics 2009, 10:41.

  4. 4.

    Spicer R, Groover A: Evolution of development of vascular cambia and secondary growth. New Phytol 2010, 186:577–592.

  5. 5.

    Hertzberg M, Aspeborg H, Schrader J, Andersson A, Erlandsson R, Blomqvist K, Bhalerao R, Uhlén M, Teeri TT, Lundeberg J, Sundberg B, Nilsson P, Sandberg G: A transcriptional roadmap to wood formation. Proc Natl Acad Sci USA 2001, 98:14732–14737.

  6. 6.

    Wang M, Qi X, Zhao S, Zhang S, Lu MZ: Dynamic changes in transcripts during regeneration of the secondary vascular system in Populus tomentosa Carr. revealed by cDNA microarrays. BMC Genomics 2009, 10:215.

  7. 7.

    Sterky F, Regan S, Karlsson J, Hertzberg M, Rohde A, Holmberg A, Amini B, Bhalerao R, Larsson M, Villarroel R, Van Montagu M, Sandberg G, Olsson O, Teeri TT, Boerjan W, Gustafsson P, Uhlén M, Sundberg B, Lundeberg J: Gene discovery in the wood-forming tissues of poplar: analysis of 5,692 expressed sequence tags. Proc Natl Acad Sci USA 1998, 95:13330–13335.

  8. 8.

    Allona I, Quinn M, Shoop E, Swope K, St Cyr S, Carlis J, Riedl J, Retzel E, Campbell MM, Sederoff R, Whetten RW: Analysis of xylem formation in pine by cDNA sequencing. Proc Natl Acad Sci USA 1998, 95:9693–9698.

  9. 9.

    Lee C-H, Chan M-H, Wang Y-N, Chu F-H: Gene investigation into the inner bark of Taiwania ( Taiwania cryptomerioides ). Bot Stud 2006, 47:111–118.

  10. 10.

    Barros E, van Staden CA, Lezar S: A microarray-based method for the parallel analysis of genotypes and expression profiles of wood-forming tissues in Eucalyptus grandis . BMC Biotechnol 2009, 9:51.

  11. 11.

    Dharmawardhana P, Brunner AM, Strauss SH: Genome-wide transcriptome analysis of the transition from primary to secondary stem development in Populus trichocarpa . BMC Genomics 2010, 11:150.

  12. 12.

    Schrader J, Moyle R, Bhalerao R, Hertzberg M, Lundeberg J, Nilsson P, Bhalerao RP: Cambial meristem dormancy in trees involves extensive remodeling of the transcriptome. Plant J 2004, 40:173–187.

  13. 13.

    Paiva JA, Garcés M, Alves A, Garnier-Géré P, Rodrigues JC, Lalanne C, Porcon S, Le Provost G, Perez Dda S, Brach J, Frigerio JM, Claverol S, Barré A, Fevereiro P, Plomion C: Molecular and phenotypic profiling from the base to the crown in maritime pine wood-forming tissue. New Phytol 2008, 178:283–301.

  14. 14.

    Nairn CJ, Lennon DM, Wood-Jones A, Nairn AV, Dean JF: Carbohydrate-related genes and cell wall biosynthesis in vascular tissues of loblolly pine ( Pinus taeda ). Tree Physiol 2008, 28:1099–1110.

  15. 15.

    Ohba K: Clonal Forestry with sugi ( Cryptomeria japonica ). In Clonal Forestry II, Conservation and Application. Edited by: Ahja M-R, Libby WJ. Berlin, Heidelberg, New York: Springer; 1993:66–90.

  16. 16.

    Forest Agency, Ministry of Agriculture, Forest and Fisheries, Japan: Annual report on forest and forestry in Japan (in Japanese). In Forest management. Edited by: Forestry agency. Tokyo: National Forestry Extension Association in Japan; 2012:71–72.

  17. 17.

    Bellin D, Ferrarini A, Chimento A, Kaiser O, Levenkova N, Bouffard P, Delledonne M: Combining next-generation pyrosequencing with microarray for large scale expression analysis in non-model species. BMC Genomics 2009, 10:555.

  18. 18.

    Vera JC, Wheat CW, Fescemyer HW, Frilander MJ, Crawford DL, Hanski I, Marden JH: Rapid transcriptome characterization for a nonmodel organism using 454 pyrosequencing. Mol Ecol 2008, 17:1636–1647.

  19. 19.

    Egertsdotter U, van Zyl LM, MacKay J, Peter G, Kirst M, Clark C, Whetten R, Sederoff R: Gene expression during formation of earlywood and latewood in loblolly pine: expression profiles of 350 genes. Plant Biol 2004, 6:654–663.

  20. 20.

    Le Provost G, Paiva J, Pot D, Brach J, Plomion C: Seasonal variation in transcript accumulation in wood-forming tissues of maritime pine ( Pinus pinaster Ait) with emphasis on a cell wall glycine-rich protein. Planta 2003, 217:820–830.

  21. 21.

    Futamura N, Totoki Y, Toyoda A, Igasaki T, Nanjo T, Seki M, Sakaki Y, Mari A, Shinozaki K, Shinohara K: Characterization of expressed sequence tags from a full-length enriched cDNA library of Cryptomeria japonica male strobili. BMC Genomics 2008, 9:383.

  22. 22.

    ForestGen. []

  23. 23.

    Bateman A, Coin L, Durbin R, Finn RD, Hollich V, Griffiths-Jones S, Khanna A, Marshall M, Moxon S, Sonnhammer ELL, Studholme DJ, Yeats C, Eddy SR: The Pfam protein families database. Nucleic Acids Res 2004, 32:D138-D141.

  24. 24.

    Pavy N, Paule C, Parsons L, Crow JA, Morency M-J, Cooke J, Johnson JE, Noumen E, Guillet-Claude C, Butterfield Y, Barber S, Yang G, Liu J, Stott J, Kirkpatrick R, Siddiqui A, Holt R, Marra M, Seguin A, Retzel E, Bousquet J, Mackay J: Generation, annotation, analysis and database integration of 16,500 white spruce EST clusters. BMC Genomics 2005, 6:144.

  25. 25.

    Tatusov RL, Fedorova ND, Jackson JD, Jacobs AR, Kiryutin B, Koonin EV, Krylov DM, Mazumder R, Mekhedov SL, Nikolskaya AN, Rao BS, Smirnov S, Sverdlov AV, Vasudevan S, Wolf YI, Yin JJ, Natale DA: The COG database: an updated version includes eukaryotes. BMC Bioinformatics 2003, 4:41.

  26. 26.

    Futamura N, Ujino-Ihara T, Nishiguchi M, Kanamori H, Yoshimura K, Sakaguchi M, Shinohara K: Analysis of expressed sequence tags from Cryptomeria japonica pollen reveals novel pollen-specific transcripts. Tree Physiol 2006, 26:1517–1528.

  27. 27.

    Nanjo T, Futamura N, Nishiguchi M, Igasaki T, Shinozaki K, Shinohara K: Characterization of full-length enriched expression sequence tag of stress-treated poplar leaves. Plant Cell Physiol 2004, 45:1738–1748.

  28. 28.

    Yoshida K, Futamura N, Nishiguchi M: Collection of expressed genes from the transition zone of Cryptomeria japonica in the dormant season. J Wood Sci 2012, 58:89–103.

  29. 29.

    Yamaguchi M, Demura T: Transcriptional regulation of secondary wall formation controlled by NAC domain proteins. Plant Biotech 2010, 27:237–242.

  30. 30.

    Guo AY, Chen X, Gao G, Zhang H, Zhu QH, Liu XC, Zhong YF, Gu XC, He K, Luo JC: PlantTFDB: a comprehensive plant transcription factor database. Nucleic Acids Res 2008, 36:D966-D969.

  31. 31.

    Itoh T, Hayashi S, Kishima T: Cambial activity and radial growth in SUGI trees (Japanese Cryptomeria). Wood Res 1968, 45:23–35.

  32. 32.

    Begum S, Nakaba S, Oribe Y, Kubo T, Funada R: Cambial sensitivity to rising temperatures by natural condition and artificial heating from late winter to early spring in the evergreen conifer Cryptomeria japonica . Trees 2010, 24:43–52.

  33. 33.

    Begum S, Nakaba S, Oribe Y, Kubo T, Funada R: Changes in the localization and levels of starch and lipids in cambium and phloem during cambial reactivation by artificial heating of main stems of Cryptomeria japonica trees. Ann Bot 2010, 106:885–895.

  34. 34.

    Begum S, Nakaba S, Yamagishi Y, Yamane K, Islam Md A, Oribe Y, Ko JH, Jin HO, Funada R: A rapid decrease in temperature induces latewood formation in artificially reactivated cambium of conifer stems. Ann Bot 2012, 110:875–885.

  35. 35.

    Oribe Y, Kubo T: Effect of heat on cambial reactivation during winter dormancy in evergreen and deciduous conifers. Tree Physiol 1997, 17:81–87.

  36. 36.

    Druart N, Johansson A, Baba K, Schrader J, Sjödin A, Bhalerao RR, Resman L, Trygg J, Moritz T, Bhalerao RP: Environment and hormonal regulation of the activity-dormancy cycle in the cambial meristem involves stage-specific modulation of transcriptional and metabolic network. Plant J 2007, 50:557–573.

  37. 37.

    Menges M, de Jager SM, Gruissem W, Murray JA: Global analysis of the core cell cycle regulators of Arabidopsis identifies novel genes, reveals multiple and highly specific profiles of expression and provides a coherent model for plant cell cycle control. Plant J 2005, 41:546–566.

  38. 38.

    Li WF, Ding Q, Chen JJ, Cui KM, He XQ: Induction of PtoCDKB and PtoCYCB transcription by temperature during cambium reactivation in Populus tomentosa Carr. J Exp Bot 2009, 69:2621–2630.

  39. 39.

    Koutaniemi S, Warinowski T, Kärkönen A, Alatalo E, Fossdal CG, Saranpää P, Laakso T, Fagerstedt KV, Simola LK, Paulin L, Rudd S, Teeri TH: Expression profiling of the lignin biosynthetic pathway in Norway spruce using EST sequencing and real-time RT-PCR. Plant Mol Biol 2007, 65:311–328.

  40. 40.

    Shi R, Sun YH, Li Q, Heber S, Sederoff R, Chiang VL: Towards a systems approach for lignin biosynthesis in Populus trichocarpa : transcript abundance and specificity of the monolignol biosynthetic genes. Plant Cell Physiol 2010, 51:144–163.

  41. 41.

    Anterola AM, Jeon JH, Davin LB, Lewis NG: Transcriptional control of monolignol biosynthesis in Pinus taeda : factors affecting monolignol ratios and carbon allocation in phenylpropanoid metabolism. J Biol Chem 2002, 277:18272–18280.

  42. 42.

    Barcelo R: Xylem parenchyma cells deliver the H 2 O 2 necessary for lignification in differentiating xylem vessels. Planta 2005, 220:747–756.

  43. 43.

    Marjamaa K, Lehtonen M, Lundell T, Toikka M, Saranpää P, Fagerstedt KV: Developmental lignification and seasonal variation in beta-glucosidase and peroxidase activities in xylem of Scots pine, Norway spruce and silver birch. Tree Physiol 2003, 23:977–986.

  44. 44.

    Hiraga S, Yamamoto K, Ito H, Sasaki K, Matsui H, Honma M, Nagamura Y, Sasaki T, Ohashi Y: Diverse expression profiles 21 rice peroxidase genes. FEBS Lett 2000, 471:245–250.

  45. 45.

    Berthet S, Demont-Caulet N, Pollet B, Bidzinski P, Cézard L, Le Bris P, Borrega N, Hervé J, Blondet E, Balzergue S, Lapierre C, Jouanin L: Disruption of LACCASE4 and 17 results in tissue-specific alterations to lignification of Arabidopsis thaliana stems. Plant Cell 2011, 23:1124–1137.

  46. 46.

    Endler A, Persson S: Cellulose synthases and synthesis in Arabidopsis . Mol Plant 2011, 4:199–211.

  47. 47.

    Taylor NG, Howells RM, Huttly AK, Vickers K, Turner SR: Interactions among three distinct CesA proteins essential for cellulose synthesis. Proc Natl Acad Sci USA 2003, 100:1450–1455.

  48. 48.

    Tanaka K, Murata K, Yamazaki M, Onosato K, Miyao A, Hirochika H: Three distinct rice cellulose synthase catalytic subunit genes required for cellulose synthesis in the secondary wall. Plant Physiol 2003, 133:73–83.

  49. 49.

    Wu L, Joshi CP, Chiang VL: A xylem-specific cellulose synthase gene from aspen ( Populus tremuloides ) is responsive to mechanical stress. Plant J 2000, 22:495–502.

  50. 50.

    Samuga A, Joshi CP: A new cellulose synthase gene ( PtrCesA2 ) from aspen xylem is orthologous to Arabidopsis AtCesA7 (irx3) gene associated with secondary cell wall synthesis. Gene 2002, 296:37–44.

  51. 51.

    Joshi CP: Xylem-specific and tension stress-responsive expression of cellulose synthase genes from aspen trees. Appl Biochem Biotechnol 2003, 105–108:17–25.

  52. 52.

    Bhandari S, Fujino T, Thammanagowda S, Zhang D, Xu F, Joshi CP: Xylem-specific and tension stress-responsive coexpression of KORRIGAN endoglucanase and three secondary wall-associated cellulose synthase genes in aspen trees. Planta 2006, 224:828–837.

  53. 53.

    Nairn CJ, Haselkorn T: Three loblolly pine CesA genes expressed in developing xylem are orthologous to secondary cell wall CesA genes of angiosperms. New Phytol 2005, 166:907–915.

  54. 54.

    Nicol F, His I, Jauneau A, Vernhettes S, Canut H, Höfte H: A plasma membrane-bound putative endo-1, 4-beta-D-glucanase is required for normal wall assembly and cell elongation in Arabidopsis . EMBO J 1998, 17:5563–5576.

  55. 55.

    Roudier F, Fernandez AG, Fujita M, Himmelspach R, Borner GH, Schindelman G, Song S, Baskin TI, Dupree P, Wasteneys GO, Benfey PN: COBRA, an Arabidopsis extracellular glycosyl-phosphatidyl inositol-anchored protein, specifically controls highly anisotropic expansion through its involvement in cellulose microfibril orientation. Plant Cell 2005, 17:1749–1763.

  56. 56.

    Takahashi J, Rudsander UJ, Hedenström M, Banasiak A, Harholt J, Amelot N, Immerzeel P, Ryden P, Endo S, Ibatullin FM, Brumer H, del Campillo E, Master ER, Scheller HV, Sundberg B, Teeri TT, Mellerowicz EJ: KORRIGAN1 and its aspen homolog PttCel9A1 decrease cellulose crystallinity in Arabidopsis stems. Plant Cell Physiol 2009, 50:1099–1115.

  57. 57.

    Maloney VJ, Samuels AL, Mansfield SD: The endo-1, 4-beta-glucanase Korrigan exhibits functional conservation between gymnosperms and angiosperms and is required for proper cell wall formation in gymnosperms. New Phytol 2012, 193:1076–1087.

  58. 58.

    Peña MJ, Zhong R, Zhou G-K, Richardson EA, O’Neill MA, Darvill AG, York WS, Ye Z-H: Arabidopsis irregular xylem8 and irregular xylem9 : implications for the complexity of glucuronoxylan biosynthesis. Plant Cell 2007, 19:549–563.

  59. 59.

    Brown DM, Goubet F, Wong VW, Goodacre R, Stephens E, Dupree P, Turner SR: Comparison of five xylan synthesis mutants reveals new insight into the mechanisms of xylan synthesis. Plant J 2007, 52:1154–1168.

  60. 60.

    Brown DM, Zhang Z, Stephens E, Dupree P, Turner SR: Characterization of IRX10 and IRX10-like reveals an essential role in glucuronoxylan biosynthesis in Arabidopsis . Plant J 2009, 57:732–746.

  61. 61.

    Brown D, Wightman R, Zhang Z, Gomez LD, Atanassov I, Bukowski JP, Tryfona T, McQueen-Mason SJ, Dupree P, Turner S: Arabidopsis genes IRREGULAR XYLEM ( IRX15 ) and IRX15L encode DUF579-containing proteins that are essential for normal xylan deposition in the secondary cell wall. Plant J 2011, 66:401–413.

  62. 62.

    York WS, O’Neill MA: Biochemical control of xylan biosynthesis — which end is up? Curr Opin Plant Biol 2008, 11:258–265.

  63. 63.

    Wu A-M, Rihouey C, Seveno M, Hörnblad E, Singh SK, Matsunaga T, Ishii T, Lerouge P, Marchant A: The Arabidopsis IRX10 and IRX10-LIKE glycosyltransferases are critical for glucuronoxylan biosynthesis during secondary cell wall formation. Plant J 2009, 57:718–731.

  64. 64.

    Rose JKC, Braam J, Fry SC, Nishitani K: The XTH family of enzymes involved in xyloglucan endotransglucosylation and endohydrolysis: current perspectives and a new unifying nomenclature. Plant Cell Physiol 2002, 43:1421–1435.

  65. 65.

    Geisler-Lee J, Geisler M, Coutinho PM, Segerman B, Nishikubo N, Takahashi J, Aspeborg H, Djerbi S, Master E, Andersson-Gunnerås S, Sundberg B, Karpinski S, Teeri TT, Kleczkowski LA, Henrissat B, Mellerowicz EJ: Poplar carbohydrate-active enzymes. Gene identification and expression analyses. Plant Physiol 2006, 140:946–962.

  66. 66.

    Chen A, He S, Li F, Li Z, Ding M, Liu Q, Rong J: Analyses of the sucrose synthase gene family in cotton: structure, phylogeny and expression patterns. BMC Plant Biol 2012, 12:85.

  67. 67.

    Déjardin A, Sokolov LN, Kleczkowski LA: Sugar/osmoticum levels modulate differential abscisic acid-independent expression of two stress-responsive sucrose synthase genes in Arabidopsis . Biochem J 1999, 344:503–509.

  68. 68.

    Stitt M, Wilke I, Feil R, Heldt HW: Coarse control of sucrose phosphate synthase in leaves: alterations of the kinetic properties in response to the rate of photosynthesis and the accumulation of sucrose. Planta 1988, 174:217–230.

  69. 69.

    Park J-Y, Canam T, Kang K-Y, Unda F, Mansfield SD: Sucrose phosphate synthase expression influences poplar phenology. Tree Physiol 2009, 29:937–946.

  70. 70.

    Kawaoka A, Kaothien P, Yoshida K, Endo S, Yamada K, Ebinuma H: Functional analysis of tobacco LIM protein Ntlim1 involved in lignin biosynthesis. Plant J 2000, 22:289–301.

  71. 71.

    Oh S, Park S, Han KH: Transcriptional regulation of secondary growth in Arabidopsis thaliana . J Exp Bot 2003,54(393):2709–2722.

  72. 72.

    Rogers LA, Dubos C, Surman C, Willment J, Cullis IF, Mansfield SD, Campbell MM: Comparison of lignin deposition in three ectopic lignification mutants. New Phytol 2005, 168:123–140.

  73. 73.

    Groover A, Robischon M: Developmental mechanisms regulating secondary growth in woody plants. Curr Opin Plant Biol 2006, 9:55–58.

  74. 74.

    Du J, Groover A: Transcriptional regulation of secondary growth and wood formation. Trends Plant Sci 2010, 15:17–27.

  75. 75.

    Zhong R, Lee C, Ye Z-H: Evolutionary conservation of the transcriptional network regulating secondary cell wall biosynthesis. Trends Plant Sci 2010, 15:625–632.

  76. 76.

    Pimrote K, Tian Y, Lu X: Transcriptional regulatory network controlling secondary cell wall biosynthesis and biomass production in vascular plants. Afr J Biotechnol 2012, 11:13928–13937.

  77. 77.

    Nakano Y, Nishikubo N, Goué N, Ohtani M, Yamaguchi Katayama Y, Demura T: MYB transcription factors orchestrating the developmental program of xylem vessels in Arabidopsis roots. Plant Biotechnol 2010, 27:267–272.

  78. 78.

    Bedon F, Grima-Pettenati J, Mackay J: Conifer R2R3-MYB transcription factors: sequence analyses and gene expression in wood-forming tissues of white spruce (Picea glauca). BMC Plant Biol 2007, 7:17.

  79. 79.

    Patzlaff A, McInnis S, Courtenay A, Surman C, Newman LJ, Smith C, Bevan MW, Mansfield S, Whetten RW, Sederoff RR, Campbell MM: Characterisation of a pine MYB that regulates lignification. Plant J 2003, 36:743–754.

  80. 80.

    Patzlaff A, Newman LJ, Dubos C, Whetten RW, Smith C, McInnis S, Bevan MW, Sederoff RR, Campbell MM: Characterisation of PtMYB1 , an R2R3-MYB from pine xylem. Plant Mol Biol 2003, 53:597–608.

  81. 81.

    Kubo M, Udagawa M, Nishikubo N, Horiguchi G, Yamaguchi M, Ito J, Mimura T, Fukuda H, Demura T: Transcription switches for protoxylem and metaxylem vessel formation. Genes Dev 2005, 19:1855–1860.

  82. 82.

    Zhong R, Demura T, Ye ZH: SND1 , a NAC domain transcription factor, is a key regulator of secondary wall synthesis in fibers of Arabidopsis . Plant Cell 2006, 18:3158–3170.

  83. 83.

    Zhong R, Lee C, Zhou J, McCarthy RL, Ye Z-H: A battery of transcription factors involved in the regulation of secondary cell wall biosynthesis in Arabidopsis . Plant Cell 2008, 20:2763–2782.

  84. 84.

    Mitsuda N, Iwase A, Yamamoto H, Yoshida M, Seki M, Shinozaki K, Ohme-Takagi M: NAC transcription factors, NST1 and NST3, are key regulators of the formation of secondary walls in woody tissues of Arabidopsis. Plant Cell 2007, 19:270–280.

  85. 85.

    Mitsuda N, Ohme-Takagi M: NAC transcription factors NST1 and NST3 regulate pod shattering in a partially redundant manner by promoting secondary wall formation after the establishment of tissue identity. Plant J 2008, 56:768–778.

  86. 86.

    Yamaguchi M, Kubo M, Fukuda H, Demura T: Vascular-related NAC-DOMAIN7 is involved in the differentiation of all types of xylem vessels in Arabidopsis roots and shoots. Plant J 2008, 55:652–664.

  87. 87.

    Sandberg G, Ericsson A: Indole-3-acetic acid concentration in the leading shoot and living stem bark of Scots pine: seasonal variation and effects of pruning. Tree Physiol 1987, 3:173–183.

  88. 88.

    Sundberg B, Little CHA, Cui K, Sandberg G: Level of indole-3-acetic acid in the stem of Pinus sylvestris in relation to seasonal variation of cambial activity. Plant Cell Environ 1991, 14:241–246.

  89. 89.

    Funada R, Kubo T, Tabuchi M, Sugiyama T, Fushitani M: Seasonal variations in endogenous indole-3-acetic acid and abscisic acid in the cambial region of Pinus densiflora stems in relation to earlywood-latewood transition and cessation of tracheid production. Holzforschung 2001, 55:128–134.

  90. 90.

    Uggla C, Moritz T, Sandberg G, Sundberg B: Auxin as a positional signal in pattern formation in plants. Proc Natl Acad Sci USA 1996, 93:9282–9286.

  91. 91.

    Tuominen H, Puech L, Fink S, Sundberg B: A radial concentration gradient of indole-3-acetic acid is related to secondary xylem development in Hybrid Aspen. Plant Physiol 1997, 116:577–585.

  92. 92.

    Uggla C, Mellerowicz EJ, Sundberg B: Indole-3-acetic acid controls cambial growth in scots pine by positional signaling. Plant Physiol 1998, 117:113–121.

  93. 93.

    Dharmasiri N, Estelle M: Auxin signaling and regulated protein degradation. Trends Plant Sci 2004, 9:302–308.

  94. 94.

    Parry G, Estelle M: Auxin receptors: a new role for F-box proteins. Curr Opin Plant Biol 2006, 18:152–156.

  95. 95.

    Israelsson M, Sundberg B, Moritz T: Tissue-specific localization of gibberellins and expression of gibberellin-biosynthetic and signaling genes in wood-forming tissues in aspen. Plant J 2005, 44:494–504.

  96. 96.

    Mauriat M, Moritz T: Analyses of GA20ox - and GID1 -over-expressing aspen suggest that gibberellins play two distinct roles in wood formation. Plant J 2009, 58:989–1003.

  97. 97.

    Funada R, Miura T, Shimizu Y, Kinase T, Nakaba S, Kubo T, Sano Y: Gibberellin-induced formation of tension wood in angiosperm trees. Planta 2008, 227:1409–1414.

  98. 98.

    Nugroho WD, Yamagishi Y, Nakaba S, Fukuhara S, Begum S, Marsoem SN, Ko JH, Jin HO, Funada R: Gibberellin is required for the formation of tension wood and stem gravitropism in Acacia mangium seedlings. Ann Bot 2012, 110:887–895.

  99. 99.

    Fleet CM, Sun T-P: A DELLAcate balance: the role of gibberellin in plant morphogenesis. Curr Opin Plant Biol 2005, 8:77–85.

  100. 100.

    Matakiadis T, Alboresi A, Jikumaru Y, Tatematsu K, Pichon O, Renou JP, Kamiya Y, Nambara E, Truong HN: The Arabidopsis abscisic acid catabolic gene CYP707A2 plays a key role in nitrate control of seed dormancy. Plant Physiol 2009, 149:949–960.

  101. 101.

    Karlberg A, Englund M, Petterle A, Molnar G, Sjödin A, Bako L, Bhalerao RP: Analysis of global changes in gene expression during activity-dormancy cycle in hybrid aspen apex. Plant Biotechnol 2010, 27:1–16.

  102. 102.

    Gilmour SJ, Zarka DG, Stockinger EJ, Salazar MP, Houghton JM, Thomashow MF: Low temperature regulation of the Arabidopsis CBF family of AP2 transcriptional activators as an early step in cold-induced COR gene expression. Plant J 1998, 16:433–442.

  103. 103.

    Chinnusamy V, Ohta M, Kanrar S, Lee BH, Hong X, Agarwal M, Zhu JK: ICE1 : a regulator of cold-induced transcriptome and freezing tolerance in Arabidopsis . Genes Dev 2003, 17:1043–1054.

  104. 104.

    Benedict C, Skinner JS, Meng R, Chang Y, Bhalerao R, Huner NP, Finn CE, Chen TH, Hurry V: The CBF1-dependent low temperature signalling pathway, regulon and increase in freeze tolerance are conserved in Populus spp. Plant Cell Environ 2006, 29:1259–1272.

  105. 105.

    Park S, Keathley DE, Han KH: Transcriptional profiles of the annual growth cycle in Populus deltoids . Tree Physiol 2008, 28:321–329.

  106. 106.

    Oda Y, Hasezawa S: Cytoskeletal organization during xylem cell differentiation. J Plant Res 2006, 119:167–177.

  107. 107.

    Gutierrez R, Lindeboom JJ, Paredez AR, Emons AM, Ehrhardt DW: Arabidopsis cortical microtubules position cellulose synthase delivery to the plasma membrane and interact with cellulose synthase trafficking compartments. Nat Cell Biol 2009, 11:797–806.

  108. 108.

    Spokevicius AV, Southerton SG, MacMillan CP, Qiu D, Gan S, Tibbits JF, Moran GF, Bossinger G: beta-Tubulin affects cellulose microfibril orientation in plant secondary fibre cell walls. Plant J 2007, 51:717–726.

  109. 109.

    Oakley RV, Wang YS, Ramakrishna W, Harding SA, Tsai CJ: Differential expansion and expression of alpha- and beta-tubulin gene families in Populus . Plant Physiol 2007, 145:961–973.

  110. 110.

    McDowell JM, Huang S, McKinney EC, An Y-Q, Meagher RB: Structure and evolution of the actin gene family in Arabidopsis thaliana . Genetics 1996, 142:587–602.

  111. 111.

    Gilliland LU, Pawloski LC, Kandasamy MK, Meagher RB: Arabidopsis actin gene ACT7 plays an essential role in germination and root growth. Plant J 2003, 33:319–328.

  112. 112.

    Mao T, Jin L, Li H, Liu B, Yuan M: Two microtubule-associated proteins of the Arabidopsis MAP65 family function differently on microtubules. Plant Physiol 2005, 138:654–662.

  113. 113.

    The Cutadapt Tool. []

  114. 114.

    NCBI Non Redundant Database. []

  115. 115.

    The Arabidopsis Information Resource (TAIR). []

  116. 116.

    The Gene Index Project. []

  117. 117.

    Subio Platform Software. []

Download references


We thank Dr. Mine Nose for providing us with the endogenous control primer for qRT-PCR. This work was supported by grants from the Japan Society for the Promotion of Science Grant-in-Aid for Scientific Research (number 24380098), the Program for Promotion of Basic and Applied Researches for Innovations in Bio-oriented Industry, the Forestry and Forest Products Research Institute (number 201004), and the project “Development of mitigation and adaptation techniques to global warming in the sectors of agriculture, forestry, and fisheries” financed by the Ministry of Agriculture, Forestry and Fisheries of Japan.

Author information

Correspondence to Atsushi Watanabe.

Additional information

Competing interests

The authors declare that they have no competing interests.

Authors’ contributions

KM carried out construction of the cDNA library and microarray experiments, qRT-PCR validation, and manuscript preparation. MT conducted part of the analysis of EST data against databases. TF, TI, KK, KY and YF helped sample materials and performed anatomical observations. TF and AW proposed the research project and guided the research process. All authors have read and approved the final version of the manuscript.

An erratum to this article is available at

Electronic supplementary material

Additional file 1: Figure S1.: Chart of this study. (PDF 138 KB)

Additional file 2:: Differentially expressed targets during the growing season. A set of 10,380 C. japonica genes that were differentially expressed. Gene order is the same as in Figure 5 cluster diagrams. (XLSX 2 MB)

Additional file 3:: Description and expression profile of the individual targets listed for Figures 6,7,8,9,10,11,12,13,14,15,16,17,18,19, and 20. (XLSX 139 KB)

Additional file 4: Figure S2.: Daily maximum and minimum temperatures measured at the sampling site. Sampling days are indicated: 24 March, 27 April, 22 June, 24 August and 7 October, 2011. (PDF 132 KB)

Rights and permissions

Reprints and Permissions

About this article


  • Cryptomeria japonica
  • cDNA library
  • Microarray
  • Cambium