Effects of day length- and temperature-regulated genes on annual transcriptome dynamics in Japanese cedar (Cryptomeria japonica D. Don), a gymnosperm indeterminate species

Seasonal phenomena in plants are primarily affected by day length and temperature. The shoot transcriptomes of trees grown in the field and a controlled-environment chamber were compared to characterize genes that control annual rhythms and the effects of day length- and temperature-regulated genes in the gymnosperm Japanese cedar (Cryptomeria japonica D. Don), which exhibits seasonally indeterminate growth. Annual transcriptome dynamics were clearly demonstrated by principal component analysis using microarray data obtained under field-grown conditions. Analysis of microarray data from trees grown in a controlled chamber identified 2,314 targets exhibiting significantly different expression patterns under short-day (SD) and long-day conditions, and 2,045 targets exhibited significantly different expression patterns at 15°C (LT; low temperature) versus 25°C. Interestingly, although growth was suppressed under both SD and LT conditions, approximately 80% of the SD- and LT-regulated targets differed, suggesting that each factor plays a unique role in the annual cycle. The top 1,000 up-regulated targets in the growth/dormant period in the field coincided with more than 50% of the SD- and LT-regulated targets, and gene co-expression network analysis of the annual transcriptome indicated a close relationship between the SD- and LT-regulated targets. These results indicate that the respective effects of day length and temperature interact to control annual transcriptome dynamics. Well-known upstream genes of signaling pathways responsive to environmental conditions, such as the core clock (LHY/CjLHYb and CCA1/CjLHYa) and PEBP family (MFT) genes, exhibited unique expression patterns in Japanese cedar compared with previous reports in other species, suggesting that these genes control differences in seasonal regulation mechanisms between species. The results of this study provide new insights into seasonal regulation of transcription in Japanese cedar.

Background Day length and temperature, which exhibit particularly large annual changes, have a marked influence on seasonal phenomena, such as bud flush and formation, growth pattern, and dormancy induction [1][2][3][4][5][6]. A number of previous reports have described how annual cycles are controlled at the transcript level and the effects of day length and temperature in the model tree, Populus [5][6][7][8][9]. Although the molecular mechanisms regulating seasonal cycles in angiosperm trees such as Populus have been extensively studied, these cycles may be controlled by different mechanisms in gymnosperm trees, given that angiosperms and gymnosperms separated evolutionarily 300 million years ago [10]. Indeed, recent studies of gymnosperm trees described unique characteristics of a number of genes that function as important regulators of environmental responses. For example, while the expression of homologs of FT (FLOWERING LOCUS T), which regulates flowering time and seasonal growth, decreased in Populus (PtFT1), it increased in Picea abies (Norway spruce, PaFT4) during the period of growth cessation [11][12][13][14]. The expression of core clock genes, which play roles in adaptation to day length and temperature changes and the regulation of seasonal phenomena [15,16], was shown to be arrested under conditions of continuous light or dark in gymnosperm conifers, in contrast to angiosperms [17,18]. Differences in diurnal rhythms of clock genes between Arabidopsis and Pseudotsuga menziesii (Douglas-fir) have also been reported [19]. Collectively, these data suggest the importance of studying gymnosperms to understand their unique seasonal regulatory mechanisms.
To date, few studies of evergreen coniferous species of gymnosperms have used time-series transcriptome analyses to investigate transcripts of needles, the key perennial organ that senses changes in environmental conditions, in order to elucidate seasonal changes and the effects of environmental factors [19][20][21]. Global changes in gene expression from late summer to early winter (August to December) have been reported in Picea sitchensis (Sitka spruce) using microarrays [20]. In P. menziesii, >80% of all identified transcripts were reportedly responsive to variations in environmental conditions during the summer (May to September) in the field [21], and analysis of annual samples indicated that 58.7% of transcripts exhibited a circannual cycle [19]. Although these studies reported seasonal dynamics and estimated the effect of environmental factors on transcript levels using data collected under field conditions [19,21], there are no published reports of studies verifying the effects of individual environmental factors on seasonal transcriptome dynamics. As multiple environmental factors, including day length and temperature, change synchronously under natural conditions, it is difficult to estimate the effect of a single factor. Comparative studies of trees grown in the field and controlled-environment chambers could enable clearer estimations of the roles of individual environmental factors in regulating annual transcriptome dynamics.
Japanese cedar (Cryptomeria japonica D. Don), a major forestry species in Japan, exhibits various unique characteristics compared with other coniferous species. Japanese cedar grows continuously until environmental or internal factors cause growth to cease (indeterminate species), whereas the amount of annual growth is regulated endogenously and environmental factors have a minor effect on the timing of growth cessation in many other coniferous species (determinant species), such as those of the genera Picea and Pinus [4,16,22]. Different mechanisms may regulate the growth of indeterminate and determinate species in response to environmental factors and control annual growth, particularly during the transition to dormancy. A previous report on Japanese cedar indicated that annual growth is influenced primarily by photoperiod and temperature [3]. The duration of growth (or conversely, the timing of dormancy induction) is important during preparation for harsh winter conditions. Investigating the contribution of photoperiod and temperature to annual transcriptome dynamics may help elucidate the mechanism by which seasonal phenomena affect the growth and environmental adaptation of indeterminate coniferous species such as Japanese cedar.
The objective of this study was to characterize the features of annual rhythms and to identify the respective effects of day length and temperature on annual rhythms in Japanese cedar by analyzing the transcriptome. The genes that exhibited annual expression rhythms under field conditions were presumed to contain the genes influenced by the respective effects of day length and temperature. We had hypothesized that a gymnosperm indeterminate species may have unique regulatory mechanisms compared with other angiosperms and/or determinate tree species reported previously. Therefore, we focused particularly on features at the transition from growth to dormancy, such as growth suppression and preparation for winter, and the regulatory mechanisms of clock genes.
First, annual transcriptome dynamics of field-grown trees were investigated by analyzing shoot samples collected throughout the year. Second, as it can be difficult to determine the effects of individual environmental factors under natural conditions, cuttings grown in a controlled-environment chamber under different day lengths and temperatures were examined using time series microarray analysis. The respective effects of day length and temperature and regulated targets were identified. Finally, the contribution of day length-and temperature-regulated genes to annual transcriptome dynamics were elucidated by comparing the results obtained under field and experimental conditions. The relationship of targets regulated by the respective effects of day length and temperature in annual transcriptome dynamics was predicted via co-expression gene network analysis.

Plant materials and samples
Measuring annual height increase in Japanese cedar. The annual growth of Japanese cedar was measured to estimate the relationship between growth and environmental conditions. The temperature of the controlled-environment chamber was determined based on these data. Because of tree size, 2-year-old trees were used for accurate measurements. The height of six cuttings of a Japanese cedar plus tree clone (Godai-1) planted in 2013 at the Forest Tree Breeding Center (FTBC), Forestry and Forest Products Research Institute (FFPRI, Hitachi, Ibaraki, Japan), was measured 19 times between February 18, 2014, and December 10, 2014. Growth rate was calculated by dividing growth at each time point by the initial value obtained on February 18, 2014.
Annual time series samples for transcriptome analysis. To analyze seasonal transcriptome dynamics under natural conditions, annual time series samples were collected at 8:00 AM on 12 dates over a period of 1 year (from June 2013 to June 2014) from eight cuttings of a clone Godai-1 planted in 1999 at the FTBC (Table 1). Annual changes in day length and temperature at Hitachi are shown in S1 Fig. A 10-cm portion of the lateral branch apex was collected from three individual trees in random order at each time point (36 samples total).
Experimental time series samples for transcriptome analysis. One-year-old cuttings of the Godai-1 clone were potted on March 4, 2014, and grown in a greenhouse at the FTBC. To evaluate the effect of day length, 16 cuttings were repotted on July 7, 2014, and transferred to chambers and cultivated under long-day (LD) conditions (16 h of light and 8 h of darkness) at 25˚C (Table 1). On August 12, 2014, nine pots were transferred to chambers and cultivated under short-day (SD) conditions (8 h of light and 16 h of darkness) at 25˚C, and seven pots were left under the LD conditions. A 10-cm portion of the lateral branch apex was harvested from two cuttings at 0, 7, 21, 35, 49, and 70 days.
To analyze the effect of temperature, the 18 cuttings were transferred from the greenhouse and placed under high-temperature (HT; 25˚C under 16 h of light and 8 h of darkness) conditions on August 5, 2014, and repotted on October 24, 2014. On November 28, 2014, nine pots were transferred to low-temperature (LT; 15˚C under 16 h of light and 8 h of darkness) conditions, and nine pots were left at HT (Table 1). Based on the annual growth pattern (S2 Fig), we selected 25˚C as the control temperature, as this is the average temperature in July, when Japanese cedar exhibits continuous growth, and 15˚C was selected as the low temperature, as this is the average temperature in November, when Japanese cedar have stopped growing. In addition, 15˚C is the temperature at which Japanese cedar starts to acquire cold hardiness [23]. A 10-cm portion of the lateral branch apex was harvested from two cuttings grown at HT and LT at 0, 7, 21, 42, 63, and 90 days after transfer.
The height of three individual cuttings was measured weekly. Lateral crown images of the three cuttings were captured using a WG-II digital camera (Pentax, Tokyo, Japan), and the projected area of the images was analyzed using ImageJ 64 software (http://rsbweb.nih.gov/ij/) to estimate the shoot growth.

RNA extraction and microarray gene expression analysis
All samples were immediately frozen in liquid nitrogen and stored at −80˚C until use. Total RNA was extracted using an RNeasy Plant Mini Kit (Qiagen, Hilden, Germany) according to a previous report [24] from 500 mg of seasonal samples collected from trees planted in the field. DNase digestion was performed on-column using an RNase-free DNase set (Qiagen). Total RNA from cuttings placed in the controlled-environment chamber was first extracted by LiCl precipitation [25], and then extracted again using an RNeasy Plant Mini Kit as described above, since purified RNA was difficult to isolate using only an RNeasy Plant Mini Kit. A NanoDrop 1000 spectrophotometer (Thermo Scientific, Waltham, MA, USA) was used to measure the RNA concentration. RNA integrity was assessed using an Agilent 2100 Bioanalyzer (Agilent Technologies, Mississauga, ON, Canada).
Microarray probes were designed based on isotig sequences from next-generation sequencing (NGS) data collected from analyses of various organs (wood cambium, treetop, shoots, and male strobili) of Japanese cedar at several developmental stages and in different seasons using a Roche GS-FLX system, as reported previously [26][27][28]. The NGS data (total 552.4 Mbp, approximately 1.3 million reads) were assembled into 22,250 isotigs as described in a previous report [28]. SurePrint G3 Gene Expression Custom 8×60K Array probes (Agilent Technologies) consisted of three probe sets corresponding to 22,194 sequences designed using the base composition methodology of the eArray tool (Agilent Technologies) using the default settings. Gene annotations represented the top-scoring BLASTX hits determined using each sequence's predicted protein product as a query against the Arabidopsis protein database TAIR10-pep-20101214 in The Arabidopsis Information Resource (http://www.arabidopsis.org) with a threshold e-value of e-5 using CLC Genomic Workbench software package, version 4.1.1 (CLC bio, Aarhus, Denmark). Total RNA (200 ng) was used for microarray analyses. cRNA was amplified and labeled using a Low Input Quick-Amp Labeling Kit (Agilent Technologies), with hybridization and washing performed according to the manufacturer's instructions. The slides were scanned on a SureScan Microarray Scanner G4900DA (Agilent Technologies), and data from the scans were compiled using Agilent Feature Extraction software, version 11.5.1.1 (Agilent Technologies).

Microarray data analysis
Raw data were normalized using a 75th-percentile shift, and median log 2 -transformed ratios for each time point were normalized to baseline intensity values (normalized intensity values) using GeneSpring software, version 14.5 (Agilent Technologies). A total of 10,439 targets with a raw signal intensity �1,000 in 100% of samples for at least one of the 34 conditions were selected for further analysis.
To obtain a high-level overview of gene expression annual dynamics, principal component analysis (PCA) was performed using the 10,439 targets of the 12 seasonal time series data sets (average normalized intensity value of three replicates) using the prcomp function in R software (R Development Core Team 2015). Based on the composite scores obtained in the first step, the expression of genes under the experimental conditions (test data) was then compared to expression under natural conditions (training data) using the R function predict to estimate the effect of seasonal conditions. To estimate the physiologic difference between growth and dormant periods with respect to the expressed genes, the functional categorizations of the top 1,000 targets with positive and negative component scores in principal component 1 (PC1) were compared to the 10,439 reference targets using the PAN-THER Overrepresentation Test (Fisher's exact test, false-discovery rate [FDR] <0.05) of the PANTHER Classification System (released 20171205, http://pantherdb.org). The unique set of Arabidopsis gene IDs (e-value < e-5) was annotated to 'GO biological process complete'.
To identify genes regulated by day length, the R package maSigPro [29] was used. Significant gene expression profile differences between experimental groups in the time course transcriptome data were identified using maSigPro, based on a regression approach. Targets exhibiting significantly different expression profiles in a time-course experiment (7,21,35,49, and 70 days) under SD and LD conditions were identified using normalized intensity values (quadratic regression, FDR �10 −4 , p-value �10 −4 , R 2 �0.7). The differentially expressed targets were classified into four clusters using the hclust function of maSigPro, and median profiles of resulting clusters were plotted using maSigPro to visualize their expression patterns. The same analysis was performed to identify LT-regulated targets using time series data collected under LT and HT conditions (7,21,42,63, and 90 days). The resulting clusters were categorized based on function and compared to the 10,439 reference targets using the PANTHER Overrepresentation Test, as described above.
To demonstrate the relationship between genes that exhibited dramatic expression changes over the course of a year, a co-expression gene network was estimated from 2,000 targets exhibiting the highest standard deviation in normalized intensity value among 36 annual time series samples. Correlation coefficients among the 2,000 targets were calculated using the R function (cor()), and gene pairs exhibiting a high correlation coefficient (cc) were selected (|cc| �0.8). The co-expression gene network was illustrated using the R package igraph [30] and laid out using 'layout_with_mds' to demonstrate the roles of SD-and LT-regulated genes in the network.

Quantitative RT-PCR
To assess the reliability of the microarray data by comparing the expression patterns determined using both microarray and quantitative RT-PCR (qRT-PCR) techniques, qRT-PCR analyses were carried out for five genes exhibiting differential expression patterns under natural conditions (MFT; mother of FT and TFL1, TEM1; tempranillo 1, RNA-binding (RRM/RBD/RNP motifs) family protein, LHY/CjLHYb, HEAT SHOCK PROTEIN 81-1). The primer pairs were designed using the Oligo software package, version 7 (National Biosciences Inc., Cascade, CO, USA) (S1 Table). First-strand cDNA synthesis and qRT-PCR were performed as described in a previous report [27]. Reaction efficiency was assessed using standard curves based on a 4-fold dilution series of cDNAs synthesized from 1,000 ng of total RNA (1:6 to 1:1,536 dilutions). Each sample was tested independently and in triplicate using all primers. Transcript abundance was normalized to eukaryotic translation initiation factor 4E1 (EIF4E1) and protein phosphatase 2A subunit A2 (PP2AA2), both of which exhibited expression variance <1.5-fold under all 34 conditions in the microarray data, determined using a method described elsewhere [31]. Data obtained for each time point were compared with data obtained for the sample in the seasonal series collected on June 24, 2013. Very similar expression patterns were obtained for transcripts analyzed using both qRT-PCR and microarray techniques (S3 Fig), suggesting that the data obtained in this study were reliable.

Phylogenetic and diurnal analysis of expression of the MFT gene
Further analyses were performed to examine expression of the MFT gene, a member of the phosphatidylethanolamine-binding protein (PEBP) gene family [14], which exhibited an interesting expression pattern in the microarray data. The amino acid sequence of MFT in Japanese cedar was estimated from the nucleic acid sequence of the NGS isotig (reCj16177ex_ne:-LSW_isotig16144, e-value 5.1E-84) using CLC Genomics Workbench software, version 11.0 (CLC bio). Phylogenetic analyses were performed using the amino acid sequences of genes homologous to MFT from plant species registered in the National Center for Biotechnology Information (NCBI) using ClustalW software, version 2.1, on the DNA Data Bank of Japan (DDBJ) website in default mode (https://clustalw.ddbj.nig.ac.jp/) according to the neighborjoining method [32]. The phylogenetic tree was illustrated using MEGA X [33].
As MFT is known to exhibit circadian rhythms in other species, we performed qRT-PCR analyses to elucidate the diurnal expression pattern of MFT in Japanese cedar. Samples were collected every 4 h beginning at 4:00 PM over a period of 2 days in summer (July 30-31), and qRT-PCR methods were the same as described in another report [27]. The primers used for qRT-PCR analysis of MFT are listed in S1 Table.

Annual growth pattern in the field
Japanese cedar began to grow in May, exhibited the most active growth in June and July, and continued growing until mid-October (S2 Fig). The initial height of the trees was 47-53 cm, which increased more than 2-fold by December 10 (108-147 cm). The average weekly temperature and day length on May 26, when the trees began growing, was 16.8˚C and 14:08, respectively, and on October 15, when the trees stopped growing, these values were 18.4˚C and 11:27, respectively. Considering these data, we selected 15˚C and 25˚C as the temperatures for LT and HT conditions.

Annual transcriptome dynamics in the field
Annual transcriptome dynamics was clearly demonstrated by PCA (Fig 1). The 12 sets of annual time series data were plotted in a circle. PC1 explained 77.7% of the total variation in Transcripts making the greatest contribution to PC1 were identified according to the absolute value of their component scores. Fig 2 shows the expression profiles of the top 1,000 targets with positive component score values (PC1 + targets) and the top 1,000 targets with negative component score values (PC1 − targets) (S2 Table). The expression of PC1 + targets tended to increase in winter (December) and decrease in summer (June). In contrast, PC1 − targets exhibited an opposite expression pattern. The PC1 + and PC1 − targets were categorized according to major gene ontology (GO) functional categories for biological processes; GO terms of the lowest hierarchy are listed in Table 2. The PC1 − targets were enriched in cell wall-related categories, such as 'cellulose catabolic process', 'pectin catabolic process', 'plant-type cell wall modification', 'xyloglucan metabolic process', and 'cell cycle'; the PC1 + targets were enriched in genes related to 'starch metabolic process' and 'response to chemical'.

PLOS ONE
Effects of day length-and temperature-regulated genes on annual transcriptome dynamics in Japanese cedar

Day length-and temperature-responsive targets identified based on experimental conditions
In this study, we perform the experiment twice to identify day length-and then temperatureregulated targets. In the two experiments, Japanese cedar trees were cultivated under LD and HT conditions, in which the environmental conditions consisted of 16 h of light and 8 h of darkness at 25˚C, to compare to SD and LT conditions, respectively. Correlation coefficients of the 10,439 targets in the microarray data indicated a degree of correlation of the transcriptome between LD and HT conditions at 0, 7, and 21 days (0.6, 0.6, and 0.7, respectively). The PCA (Fig 1) also demonstrated that the transcriptomes of cuttings from trees grown under LD and HT conditions at 0, 7, and 21 days were similar to the transcriptome of cuttings in June from the field, when trees exhibited the most active growth. These results indicated the reproducibility of tree conditions used in the day length and temperature experimental time series samples. These data also indicated that the experiments accurately represented the day length and temperature responses of Japanese cedar during the growth period.
A total of 2,314 targets (22.2%) exhibited significantly different expression profiles in the time-course experiment between SD and LD conditions (S3 Table). Clustering analyses indicated that 224 targets in Cluster 1 and 815 targets in Cluster 3 tended to be up-regulated under SD conditions, with the expression of Cluster 3 targets increasing over time (Fig 3A). In contrast, a total of 1,127 targets in Cluster 2 and 148 targets in Cluster 4 tended to be down-regulated under SD conditions, with the expression of Cluster 2 targets decreasing over time. The unique set of Arabidopsis gene IDs (e-value <e-5) of the targets of the top 1,000 positive and negative scores for principal component 1 (the analyzed list) was compared to the unique set of Arabidopsis gene IDs of 10,439 targets (the reference list) to investigate the overrepresented categories (FDR �0.05). Only the lowest categories in the GO hierarchy are listed. The designations "all genes" and "gene count" indicate the number of genes in the reference list and the analyzed list for a particular annotation data category, respectively. The designation "expected" indicates the expected number of genes in the analyzed list based on the reference list. The designation "fold enrichment" indicates the enrichment in the analyzed list based on the reference list. The designations "raw P-value" and "FDR" indicate the probability that the number of genes observed occurred by chance (randomly), as determined based on the reference list. https://doi.org/10.1371/journal.pone.0229843.t002

PLOS ONE
Effects of day length-and temperature-regulated genes on annual transcriptome dynamics in Japanese cedar Functional categorization based on biological processes indicated significant GO enrichment in each cluster; the GO terms of the lowest hierarchy are listed in Table 3. Cluster 1 consisted of a significantly higher proportion of genes involved in 'translation'. Cluster 2 consisted of a significantly higher proportion of cell wall-and growth-related categories, such as 'xylan biosynthetic process', 'plant-type secondary cell wall biogenesis', 'cellulose biosynthetic process', 'pectin metabolic process', and 'cell growth'. No statistically significant results were obtained for Clusters 3 and 4. The same statistical analyses were carried out to identify LT-regulated targets using the time series microarray data for LT and HT conditions. A total of 2,045 targets (19.6%) exhibited significantly different expression profiles under LT and HT conditions (S3 Table). Of these targets, 323 targets in Cluster 1 and 692 targets in Cluster 2 tended to be down-regulated under LT conditions, and the expression of targets in Cluster 1 increased over time (Fig 3B). In contrast, 626 targets in Cluster 3 and 404 targets in Cluster 4 tended to be up-regulated under LT conditions, and the expression of Cluster 3 targets decreased over time. GO analyses indicated that Cluster 1 consisted of a significantly higher proportion of genes involved in 'protein phosphorylation', 'transcription, DNA-templated', and 'regulation of transcription, DNA-templated', whereas Cluster 2 consisted of a significantly higher proportion of genes involved in 'response to heat', including putative genes encoding heat shock proteins (Table 3). Cluster 3, which contained LT-up-regulated targets, contained a significantly higher proportion of genes involved in 'iron-sulfur cluster assembly', 'ribosome biogenesis', 'response to cadmium ion', and 'translation'. No statistically significant results were obtained for Cluster 4.

Contribution of SD-and LT-regulated targets to annual transcriptome dynamics
To determine seasonal differences at the transcriptome level in cuttings grown under SD and LD conditions, microarray data obtained from cuttings of trees grown in a controlled-environment chamber were compared to the PCA results for annual series under natural conditions ( Fig 1A). Before transfer to SD conditions (0 days), the transcriptome of the cuttings was similar to the transcriptome observed in June. The PC1 scores increased and PC2 scores decreased over time, and the transcriptome after 70 days under SD conditions was similar to the transcriptome in autumn, between September and October. In contrast, the transcriptome of cuttings from trees grown under LT conditions exhibited a low PC1 score during the 90 days of the experiment and was similar to the transcriptome in August (Fig 1B).
The percentages of SD-and LT-regulated targets within the PC1 + and PC1 − groups were also calculated to investigate their contribution to annual transcriptome dynamics (Fig 4). Overall, PC1 + targets included 336 targets up-regulated under SD conditions and 268 up-regulated under LT conditions (50.2% of total), whereas PC1 − included 512 targets down-regulated under SD conditions and 53 down-regulated under LT conditions (54.3% of total). of SD-regulated targets or for Cluster 4 of LT-regulated targets. The designations "all genes" and "gene count" indicate the number of genes in the reference list and the analyzed list for a particular annotation data category, respectively. The designation "expected" indicates the expected number of genes in the analyzed list based on the reference list. The designation "fold enrichment" indicates the enrichment in the analyzed list based on the reference list. The designations "raw P-value" and "FDR" indicate the probability that the number of genes observed occurred by chance (randomly), as determined based on the reference list. https://doi.org/10.1371/journal.pone.0229843.t003

PLOS ONE
Effects of day length-and temperature-regulated genes on annual transcriptome dynamics in Japanese cedar Gene co-expression network analyses of the 1,976 targets for which the gene pairs exhibited highly correlated expression patterns revealed that 1,953 of the targets formed a large cluster (Fig 5). These 1,953 clustered targets included 892 targets (45.7%) of SD up-and down-regulated targets (SD-regulated targets) and 256 targets (13.1%) of LT up-and down-regulated

PLOS ONE
Effects of day length-and temperature-regulated genes on annual transcriptome dynamics in Japanese cedar targets (LT-regulated targets); a total of 117 targets (6.0%) were both SD-and LT-regulated. The network analysis indicated a mixture of SD-and LT-regulated targets, indicating that these genes contribute significantly to annual transcriptome dynamics.

Effects of day length and temperature on genes related to growth and starch metabolism
PC1 − targets that tended to exhibit increased expression in summer (Fig 2) also exhibited significant enrichment of GO terms related to growth ( Table 2). The effects of day length and temperature on the expression of genes in growth-related categories were investigated by analyzing the ratio of targets down-regulated under SD conditions to those down-regulated under LT conditions. The GO term 'cellulose biosynthetic process' was associated with 11 PC1 − targets, all of which were down-regulated under SD conditions ( Table 4). The terms 'pectin catabolic process', 'plant-type cell wall modification', and 'lignin biosynthetic process' were associated with 13, 17, and 35 PC1 − targets, respectively, and 11 (84.6%), 13 (76.5%), and 20 (57.1%) of these targets were down-regulated under SD conditions. In contrast, only two targets in these four categories were down-regulated under LT conditions, both of which were genes encoding S-adenosylmethionine synthetase family protein (MTO3) in the category 'lignin biosynthetic process'.
The PC1 + targets exhibited significant enrichment in the term 'starch metabolic process' ( Table 2). A total of 42 targets related to starch breakdown and starch synthesis were identified that were described in a previous report from a study involving Arabidopsis [34], and 18 of these targets (42.9%) were PC1 + , and their expression tended to increase in winter (Fig 6,  Table 5). Of the seven starch synthesis-related gene targets, two were up-regulated under SD conditions and three under LT conditions. The expression of two other targets, ADPGLC-PPase large subunit (APL2) and ADP glucose pyrophosphorylase 1 (APS1), which was not up-regulated under SD or LT conditions, peaked in February and then decreased dramatically from April to May (Fig 6A). Of the 11 starch breakdown-related gene targets, three were up-regulated under SD conditions and six under LT conditions, and their expression peaked in December (Fig 6B, Table 5). The expression of other targets that were not up-regulated under either SD or LT conditions tended to be high from December to March (Fig 6B).

Seasonal transcriptome dynamics and contribution of day length-and temperature-regulated genes in Japanese cedar
Japanese cedar, which is as an indeterminate species [4], exhibited continuous growth under LD and HT conditions in the controlled-environment chamber (S4 Fig) and continued growing from May through the middle of October, when the day length and temperature declined (day length 11:27, temperature 18.4˚C) (S2 Fig). Whereas the growth of determinate species is controlled endogenously and day length and temperature play minor roles in growth cessation, the growth of indeterminate species may be more susceptible to both of these environmental factors [4,16,22].
Our microarray data demonstrated dramatic changes in transcripts during the year. To the best of our knowledge, this is the first report describing annual transcriptome dynamics in a gymnosperm indeterminate species. Based on comparisons of all combinations of the 12 sets of annual time series data (Table 1), 28.8% of transcripts (3,004 targets) exhibited expression differences �5-fold in at least one combination. PCA of the microarray data demonstrated continuous changes in the transcriptome throughout the year that appeared as an annual circle when plotted (Fig 1). Among the 12 samples acquired during the year, the December sample exhibited the highest and the June sample the lowest PC1 score, which explained 77.7% of the total variation in gene expression. As day length was greatest in June and least in December (S1 Fig), the PCA results suggest that day length plays an important role in annual transcriptome dynamics. By analyzing trees grown under different experimental conditions (Table 1), we identified 2,314 SD-regulated targets and 2,045 LT-regulated targets, each of which accounted for approximately 20% of the 10,439 targets analyzed. Interestingly, only 429 targets were regulated by both SD and LT. A previous study using transgenic plants incapable of sensing short days indicated that short day and low temperature conditions induced cold acclimation by independent pathways in Populus [35]. Our data also appear to suggest these environmental factors are associated with different regulatory mechanisms. Growth cessation and bud formation were induced by SD conditions, whereas growth cessation and development of reddish-brown needles were induced by LT, in agreement with previous reports on Japanese cedar and other coniferous species (S4 Fig) [3,5,[36][37][38]. These traits enable the trees to survive in harsh winter environments.
The expression patterns of SD-and LT-regulated targets under natural conditions revealed their contribution to annual transcriptome dynamics. In this study, adult trees (15 to 16 years old) were used as annual time series samples, and juvenile trees (1 year old) were used as experimental time series samples due to limitations in chamber size. Although differences in phenology between juvenile and adult trees have not been investigated in Japanese cedar, it has been reported in temperate deciduous forest tree species [39]. The genes that exhibited annual expression in adult trees in the field and regulated in juvenile trees under environmental conditions may be genes universally regulated by day length and temperature regardless of tree age. The PC1 − and PC1 + targets that exhibited dramatic annual cycles in expression are hypothesized to be primarily affected by day length and temperature under field conditions. We found that SD-and LT-regulated targets comprised more than half of the PC1 − and PC1 + targets (Figs 2 and 4). This means that these targets were affected by either day length or temperature and that interactions between day length and temperature were not necessary to exhibit annual expression. Specifically, 336 targets were up-regulated by SD without LT and 268 targets were up-regulated by LT without SD within the PC1 + targets (Fig 4A), and 512 targets were down-regulated by SD without LT and 53 targets were down-regulated by LT without SD within the PC1 − targets (Fig 4B). The remaining targets could primarily include genes that were up-regulated only by the interaction between day length and temperature. 3.25E-55 The targets of GO categories 'pectin catabolic process', 'plant-type cell wall modification', 'cellulose biosynthetic process', and 'lignin biosynthetic process', which constituted the top 1,000 negative scores for principal component 1, are listed. � Indicates up-or down-regulation under SD (day length) or LT (temperature) conditions. https://doi.org/10.1371/journal.pone.0229843.t004 The largest cluster in the annual transcriptome dynamics gene co-expression network consisted of 97.7% of the 2,000 analyzed targets, with the SD-and LT-regulated targets located in close proximity in the network (Fig 5). These observations indicate that annual transcriptome dynamics may be regulated by both the interaction between day length and temperature and the respective effects of day length and temperature. Although day length and temperature appear to regulate different genes, there is a close relationship between day length-and temperature-related genes in annual transcriptome dynamics.
PCA of the microarray data for trees grown in the controlled-environment chamber also indicated that SD-regulated genes contribute significantly to transition to the dormant state.

PLOS ONE
Effects of day length-and temperature-regulated genes on annual transcriptome dynamics in Japanese cedar The PC1 score of the transcriptome increased and the PC2 score decreased over time for SDregulated targets, indicating that the transcriptome was shifted toward a state of dormancy by SD conditions (Fig 1A). In contrast, the LT-regulated targets exhibited a low PC1 score that over the course of the 90-day experiment was similar to the transcriptome in June and August, when Japanese cedar exhibits maximal growth (Fig 1B). These data demonstrated that SD-regulated genes contribute significantly to the transition to initial dormancy at the transcript level. Although growth cessation was affected by both SD and LT conditions (S4 Fig), LT-regulated targets contributed less to the transition to dormancy at the transcript level.

Growth-related genes are down-regulated by short-day rather than lowtemperature conditions
Although growth was suppressed under both SD and LT conditions, the down-regulation of growth-related genes in autumn may be mostly attributable to day length. The genes in cell wall-related GO categories, such as 'pectin catabolic process', 'plant-type cell wall modification', 'lignin biosynthetic process', and 'cellulose biosynthetic process', constituted a significantly higher proportion of PC1 − targets, the expression of which increased during the growth period ( Table 2). Most of the related genes were classified in Cluster 2 of the SD-regulated targets (Table 4), the expression of which decreased under SD conditions (Fig 3A).
Pectin is a structural heteropolysaccharide present in the primary cell walls of plants and contributes to complex physiologic processes such as cell growth and differentiation [40]. Most genes categorized as 'pectin catabolic process' were down-regulated under SD conditions (84.6%), including genes encoding pectin lyase family proteins and the plant invertase/pectin methylesterase inhibitor superfamily (Table 4).
Expansins promote cell wall loosening and extension and are encoded by a superfamily of genes [41]. The 'plant-type cell wall modification' category included 14 expansin gene targets, 13 of which were down-regulated under SD conditions. In addition, all of the genes related to 'cellulose biosynthetic process' and 57.1% of those related to 'lignin biosynthetic process' were down-regulated under SD conditions (Table 4). Shortened day length in autumn may repress expression of these genes related to cell growth, and may even block cell growth, leading to growth cessation. In contrast, only two genes encoding S-adenosylmethionine synthetase family proteins within these four GO categories were down-regulated under LT conditions (Table 4). Growth suppression under LT conditions may involve a different mechanism compared with SD conditions.

Only a few starch-related genes are induced by SD conditions in Japanese cedar
The expression of starch-related genes may induce starch synthesis or breakdown to enhance cold and frost tolerance and contribute energy for bud breaking and shoot growth in spring [42][43][44]. In this study, the expression of 18 of 42 analyzed genes involved in starch synthesis and breakdown (Table 5) increased in winter (December-February) as the day length and temperature decreased in the field (Fig 6). This expression pattern agreed with that reported previously for the cambial region of Japanese cedar [26]. When cambial activity and new xylem formation are activated between March and October, most homologs of starch degradationrelated genes in the cambial region exhibit minimum expression [26]. Seasonal dynamics of carbohydrate storage may be controlled by these starch-related genes. Seasonal dynamics of carbohydrate storage in Japanese cedar seedlings has been studied previously by measuring starch and sugar concentrations in the upper, middle, and lower parts of shoots and rootlets [45]. The starch concentration exhibited a small peak in autumn (November) and a large peak in spring (May), and the sugar concentration exhibited a peak in winter (January) in every part of the seedlings. The starch synthesized in autumn may be converted into sugar during the winter to enable survival in the harsh environment. These phenomena are in agreement with the annual expression pattern of starch-related genes observed in this study.
The annual dynamics of starch-related gene expression may be affected by environmental factors. In Populus and Picea glauca (white spruce), starch synthesis-and breakdown-related genes are induced by shorter day length [8,9,46]. In P. glauca, starch-synthesis genes such as starch synthase, ADP-glucose pyrophosphorylase large subunit, and starch branching enzyme are up-regulated under short day conditions [46]. However, only the starch synthase 3 gene was up-regulated in Japanese cedar under SD conditions in the present study (Table 5). Previous reports indicated that the level of freezing tolerance induced by short day-or low temperature-conditions varies with plant species [47][48][49][50]. As the expression of starch-related genes may be closely related to freezing tolerance, variations between species in the expression levels of starch-related genes under SD conditions may indicate that different processes are involved in induction of freezing tolerance. Only 2 of 18 targets related to starch synthesis and 4 of 24 targets related to starch breakdown were up-regulated under SD conditions (14.3%). Although the temperature under the LT condition was not as low as the winter temperature, more genes were up-regulated under these conditions (33.3%). Low temperature may play an important role in regulating the expression of starch-related genes and the induction of freezing tolerance in Japanese cedar.

Comparing the expression pattern of core clock-related genes to other species
Photoperiod sensing by light receptors is closely related to the circadian clock [5,51]; the diurnal pattern of the transcriptome changes in response to photoperiod variations, and the amplitude of the cycles of clock components is influenced by temperature [52][53][54]. Among the nine clock-related genes analyzed in this study, seven were affected by SD conditions (Table 6). LHY/CjLHYb (late elongated hypocotyl, CjLHYb in [27]) and ZTL (zeitlupe) were up-regulated, and PRR3 (pseudo-response regulator 3) and COL9 (constans-like 9) were down-regulated under LT conditions. Most clock-related genes exhibited high expression in December-January, except CCA1/CjLHYa (circadian clock associated 1, CjLHYa in [27]) and COL4 (constans-like 4) (Fig 7). The high expression of these genes in winter was consistent with our previous report that found dampening of diurnal rhythms and high expression of clock genes in winter [27]. The GI (gigantea), LHY/CjLHYb, and ZTL genes exhibited correlation coefficients >0.8 and were located in close proximity in the co-expression gene network, suggesting that they play similar roles in modulating annual dynamics (Fig 5). In contrast, PRR3 and PRR7 (pseudo-response regulator 7) were located at a distance from these genes.
Interestingly, the expression patterns of the core clock-related genes differed in comparison to other species reported previously. Whereas LHY is reportedly up-regulated under SD conditions in Populus [8], both homologs of LHY (LHY/CjLHYb and CCA1/CjLHYa) were significantly down-regulated under SD conditions in Japanese cedar (Table 6). Moreover, although no annual expression rhythms were reported for LHY and TOC1 (timing of Cab expression 1) in P. menziesii [19], one of the LHY homologs (LHY/CjLHYb) and TOC1 exhibited annual rhythms and high expression in winter in Japanese cedar (Fig 7). Differences in the diurnal expression pattern of core clock-related genes between P. menziesii, Arabidopsis, and Japanese cedar were described in a previous report [19]. These differences in diurnal rhythms may reflect differences in annual expression patterns. These results suggest that different regulatory mechanisms exist between species, which may influence the annual growth pattern.
The expression of two COL genes identified in the conifer P. abies (PaCOL1 and PaCOL2) was shown to decrease significantly in needles and shoot tips under SD conditions prior to growth cessation and bud formation [55]. It was hypothesized that PaCOL1 and PaCOL2 are not functional homologues of Arabidopsis CO [5], which promotes flowering in response to long day conditions [56]. The classification of PaCOL1 and PaCOL2 and their expression pattern were similar to COL of the moss Physycomitrella (PpCOL1) [5,57]. Our results agree with the reported expression pattern of COL genes in P. abies [55], indicating the conservation of regulatory mechanisms among conifer species. Expression of the putative COL9 gene in Japanese cedar was down-regulated under SD conditions (Table 6), decreased from August to November, and increased sharply in January under field conditions (Fig 7). In contrast, no differences in COL4 expression were observed between SD/LD and LT/HT conditions, with stable expression throughout the year ( Table 6, Fig 7).
flowering plants [59]. Based on phylogenetic analyses of plant MFT amino acid sequences, MFT of Japanese cedar is similar to PaMFT1 and PaMFT2 of P. abies, in which their expression is not correlated with bud setting [12,14] (S5A Fig). In P. abies, PaFTL2/PaFT4 expression increases after plants are transferred to SD conditions, in contrast to Populus, and the gene regulates bud setting and annual growth rhythm [12,14]. In wheat and Arabidopsis, MFT promotes seed dormancy [60,61]. A previous study reported that the expression of MFT did not change under SD conditions during a 6-week experiment in Populus [8], whereas another study reported that expression increased for 28 days and subsequently decreased under SD conditions in P. glauca [44]. In contrast to PtFT1, PaFTL2/PaFT4, and AcMFT, MFT in Japanese cedar did not exhibit a diurnal rhythm in July in our study (S5B Fig). These results are suggestive of evolutionary differentiation in the function of PEBP gene family members. Although we collected sequence data from various parts of trees in various developmental stages and seasons (19 libraries, approximately 3 million reads, 34,731 isotigs in total) using an NGS approach [28], no homolog of FT in Japanese cedar was identified.

Conclusion
This study provides considerable insight into the control of phenology in Japanese cedar, a gymnosperm indeterminate species. Our microarray data demonstrated dynamic annual changes in the transcriptome and the significant contribution of SD-and LT-regulated genes. More than half of the top 1,000 up-regulated targets in the growth/dormant period in the field were regulated by day length or temperature solely, indicating that interaction between day length and temperature is not necessary for their expression. Interestingly, we found that SD and LT conditions generally regulate different genes. However, the co-expression network of annual transcriptome dynamics revealed a close relationship between SD-and LT-regulated genes. These results indicate that the respective effects of SD and LT conditions play different roles but interact to regulate annual transcriptome dynamics. Compared with previous reports in other plant species, Japanese cedar exhibited several unique characteristics. Although starchrelated genes are known to be regulated by day length in other species, only a few starch-related genes were up-regulated in Japanese cedar. Upstream signaling pathway genes in the clock and PEBP family also exhibited unique expression patterns. Our data indicate that there are different seasonal regulatory mechanisms among tree species. There are several possible reasons for the variation among species: traits obtained during evolution (i.e., gymnosperms and angiosperms), differences in seasonal growth pattern (i.e., determinate and indeterminate species), and adaptation to environmental conditions in the distribution area. As our results indicate that not all regulatory mechanisms are conserved among species, investigations in other tree species would be of interest. The neighbor-joining method [32] was used to construct the phylogenetic trees. The species names are abbreviated as follows: At, Arabidopsis thaliana (thale cress); Pt, Populus trichocarpa (black cottonwood); Os, Oryza sativa (Japanese rice); Pa, Picea abies (Norway spruce); Ps, Picea sitchensis (Sitka spruce); Cj, Japanese cedar (Cryptomeria japonica); Ac, Adiantum capillus-veneris (pteridophyte); Pp, Physcomitrella patens subsp. patens (moss). The number following the species name indicates its NCBI accession number. (B) Diurnal expression of MFT in summer (July 30-31, 2012), analyzed using qRT-PCR as described in [27]. (TIFF) S1 Table. Primers used for quantitative RT-PCR in this study. � Reference [27]. (XLSX) S2 Table. Top-scoring targets for principal component 1. SD-up, SD-down, LT-up, and LTdown indicate SD up-regulated targets, SD down-regulated targets, LT up-regulated targets, and LT down-regulated targets, respectively. (XLSX) S3 Table. Day length-and temperature-regulated targets.