Changes in annual transcriptome dynamics of a clone of Japanese cedar (Cryptomeria japonica D. Don) planted under different climate conditions

Abstract

Environmental responses are critical for plant growth and survival under different climate conditions. To elucidate the underlying biological mechanisms of environmental responses in Japanese cedar (Cryptomeria japonica D. Don), the annual transcriptome dynamics of common clonal trees (Godai1) planted at three different climate sites (Yamagata, Ibaraki, and Kumamoto Prefectures) were analyzed using microarrays. Both principal component analysis (PCA) and hierarchical clustering of the microarray data indicated the transition to dormant transcriptome status occurred earlier and the transition to active growth status later in the colder region. Interestingly, PCA also indicated that the transcriptomes of trees grown under three different conditions were similar during the growth period (June to September), whereas the transcriptomes differed between sites during the dormant period (January to March). In between-site comparisons, analyses of the annual expression profiles of genes for sites ‘Yamagata vs. Kumamoto’, ‘Yamagata vs. Ibaraki’, and ‘Ibaraki vs. Kumamoto’ identified 1,473, 1,137, and 925 targets exhibiting significantly different expression patterns, respectively. The total of 2,505 targets that exhibited significantly different expression patterns in all three comparisons may play important roles in enabling cuttings to adapt to local environmental conditions. Partial least-squares regression analysis and Pearson correlation coefficient analysis revealed that air temperature and day length were the dominant factors controlling the expression levels of these targets. GO and Pfam enrichment analyses indicated that these targets include genes that may contribute to environmental adaptation, such as genes related to stress and abiotic stimulus responses. This study provided fundamental information regarding transcripts that may play an important role in adaptation to environmental conditions at different planting sites.

Introduction

Based on analyses of tree-height and tree-ring data, many studies have reported that tree growth is affected by environmental conditions [1–4]. As Japanese cedar (Cryptomeria japonica D. Don, also known as ‘Sugi’ in Japan) is a major forestry species in Japan, accounting for 44% of artificial forests [5], and is widely planted throughout the Japanese archipelago, understanding the effects of environmental conditions on this species is very important. The relationship between environmental conditions and phenotypic traits in Japanese cedar has been examined in several studies. A positive correlation between annual ring growth and temperature during February and March was reported by Takahashi [6], and significant correlations between tree age-height and several climatic variables such as warmth, solar radiation, precipitation, and snow depth were reported by Nishizono et al. [7]. Quantitative trait locus (QTL) analysis of several phenotypic traits in three replicated common garden experiments under three different climate conditions identified an average of 53 QTLs [8]. However, only two of these QTLs affected the same traits across all three environments [8]. Moreover, one QTL associated with growth response to drought stress was identified by QTL analysis of trees growing in two contrasting environments [9]. The small number of QTLs common between regions and a QTL with a large contribution to overall climate sensitivity would indicate that Japanese cedar is highly sensitive to environmental differences. Although these previous reports suggested that environmental conditions affect the phenotypic traits of Japanese cedar, there are no published reports of studies examining the underlying biological mechanisms.

Japanese cedar grows continuously until environmental factors (temperature, day length, etc.) or internal factors (aging, etc.) cause growth to cease (indeterminate growth), whereas the degree of annual growth of many other coniferous species, such as those of the genera Picea and Pinus, is regulated endogenously (determinate growth) [10–12]. As such, environmental conditions may exert a greater influence on the timing of growth cessation in Japanese cedar than in determinate growth species. Results from studies using growth chambers indicated that a short photoperiod and low temperature suppressed height growth of Japanese cedar [13, 14]. Prior to growth cessation in autumn, Japanese cedar trees accumulate starch and break it down into sugar in all parts of seedlings (upper, middle, and lower layer shoots, and roots) to enhance cold and frost tolerance as a means of ensuring survival under harsh winter conditions, and this process also contributes energy for bud breaking and shoot growth in the following spring [15, 16]. Soluble nitrogen and total free amino acids in the lateral shoots, stem, and roots increase at the time of growth initiation, but the levels decrease beginning in March in the roots and in April in lateral shoots expanded in the previous year, and in stems [17]. Since new lateral shoots exhibit high accumulation of soluble nitrogen and total free amino acids in May, stored nitrogenous compounds may be transported to new lateral shoots for growth [17].

A few studies of evergreen coniferous species have used time-series transcriptome analyses to investigate transcripts of needles and shoots, the key perennial organs that sense changes in environmental conditions, and these analyses revealed dramatic annual dynamics [14, 18–20]. The annual transcriptome dynamics of Japanese cedar were also clearly demonstrated in our previous report [14]. Principal component analysis of microarray data demonstrated the seasonal cycle of the transcriptome and explained the seasonal phenomena of Japanese cedar. For example, the expression of growth-related genes was up-regulated during the active growth period, and the expression of genes associated with starch metabolic processes that may contribute to cold tolerance was up-regulated during the dormant period. As the respective effects of day length and temperature interact to control this annual transcriptome dynamics [14], the dynamic pattern may change under different environmental conditions, leading in turn to physiological and phenological changes that help trees adapt and survive in different planting sites.

To investigate the changes in annual transcriptome dynamics of Japanese cedar under different environmental conditions, we planted specimens of a single clone at three sites covering a wide range of environmental conditions of natural distribution and analyzed the transcriptome. In this study, microarrays were used for transcriptome analysis. Although microarrays cannot detect novel sequences and splice variants as RNA-seq can, it is a high-throughput, cost-effective, and highly sensitive technology [21] suitable for the purpose of this study. To the best of our knowledge, this is the first report of environmental responses investigated by analyzing annual transcriptome dynamics in multiple sites in a tree species. Our data identified differentially expressed genes that may be related to molecular mechanisms controlling physiological and/or phenological differences in a clone of Japanese cedar between planting sites.

Materials and methods

Test sites and plant materials

One-year-old rooted cuttings of a Japanese cedar plus-tree clone, Godai1, were planted from March to April 2013 at nurseries in three locations, Yamagata (Higashine City, Yamagata Prefecture, Japan [38°23’53"N, 140°22’47"E]), Ibaraki (Hitachi City, Ibaraki Prefecture, Japan [36°41’28"N, 140°41’21"E]), and Kumamoto (Koshi City, Kumamoto Prefecture, Japan [32°52’53"N, 130°44’5"E]) (Fig 1A). The clone ‘Godai1’ was selected from Kimitsu, Chiba (35°13’12.0"N 140°07’48.0"E) as a first-generation plus tree. We have been studying ‘Godai1’ as a model clone of Japanese cedar [14], as it shows average growth and high rooting ability for propagation [22]. A total of 10, 6, and 8 rooted cuttings were planted at more than 25-cm intervals at the respective sites. No obstacles hindered sunlight from reaching the rooted cuttings at the nurseries at the three sites. Data regarding air temperature and day length at the three sites are shown in Fig 1B and 1C. Yamagata is located in northern Japan, which is a cold region with snow cover in winter; Ibaraki is located along the coast of the Pacific Ocean and has a mild climate; and Kumamoto is located in southern Japan, which has a warm climate. Throughout the year of the study, Kumamoto exhibited the highest mean air temperature and Yamagata the lowest mean air temperature among the three sites. The mean air temperature in Ibaraki was intermediate between that of Yamagata and Kumamoto. In Ibaraki, the temperature in summer (June to September) was similar to that in Yamagata, and the temperature in winter (January to March) was similar to that in Kumamoto.

[Figure omitted. See PDF.]

Sampling sites (A), mean air temperature (B), day length (C), and annual height growth of Japanese cedar (D) at the three study sites. We edited the map of Japan obtained from the GSI Maps Vector provided by the Geospatial Information Authority of Japan (https://maps.gsi.go.jp/vector/).

A 10-cm-long apex portion of sunny upper branches (S1 Fig) was randomly collected each month from two of the planted trees at 10:00-11:00 am from February 2014 to February 2015 (S1 Table). All sampling at the three sites was performed within 7 days. Meteorological data for the sampling days are shown in S2 Table. The height of all cuttings at the three sites was also measured at each sampling time point and other times, and growth rate was calculated by dividing the height growth differential from April 2014 by the initial height in April 2014. The significance of differences in annual growth rate (growth rate from April to December 2014) among the three sites was tested using one-way analysis of variance (ANOVA).

RNA extraction and microarray gene expression analysis

All collected samples (3 sites × 13 time points × 2 ramets) were immediately frozen in liquid nitrogen and stored at −80°C until analysis. Total RNA was extracted from samples as reported by Gehrig et al. [23] using an RNeasy Plant Mini kit (Qiagen, Hilden, Germany), and DNase digestion was performed on-column using an RNase-free DNase set (Qiagen). A NanoDrop 1000 spectrophotometer (Thermo Scientific, Waltham, MA, USA) was used to accurately measure RNA concentration. RNA integrity was assessed using an Agilent 2100 bioanalyzer (Agilent Technologies, Mississauga, ON, Canada).

The microarray probes were designed using e-array (Agilent Technologies) with the default Base Composition Methodology based on isotig sequences determined from next-generation sequencing data collected from analyses of various organs (cambium region, sapwood and heartwood, tree-top, shoot, male strobili, roots of seedling) in several developmental stages and seasons using a Roche 454 Genome Sequencer [24]. A SurePrint G3 Gene Expression Custom 8×60K Array (Agilent Technologies) consisting of three probe sets corresponding to 19,304 sequences was used for microarray analyses (GenBank accession number, GPL21366) [25]. Gene annotations shown represent the top-scoring BLASTX hits for each sequence’s predicted protein product as a query against the TAIR Arabidopsis protein database TAIR10-pep-20101214, with a threshold e-value of 10⁻⁵, using the CLC Genomic Workbench software package, version 4.1.1 (CLC bio, Aarhus, Denmark).

A total of 200 μg of total RNA was transcribed to double-stranded (ds) cDNA, and the cRNA was amplified and labeled using a Low-Input Quick-Amp Labeling kit (Agilent Technologies). The cRNA was hybridized to the custom gene expression microarray for 17 h at 65°C and washed using a Gene Expression Hybridization kit (Agilent Technologies) according to the manufacturer’s instructions. The resulting slides were scanned using a SureScan Microarray Scanner G4900DA (Agilent Technologies), and scan data were compiled using Agilent Feature Extraction Software 11.5.1.1 (Agilent Technologies).

Microarray data analysis

Raw microarray data were normalized using a 75th-percentile shift and analyzed using GeneSpring software, version 13.1 (Agilent Technologies). To enable direct comparisons of transcript profiles, median log₂-transformed ratios for each time point were normalized to baseline (normalized intensity value). A total of 13,385 targets with a raw signal value exceeding 100 in each pair of samples under at least one of the 39 conditions (3 sites × 13 time points) were selected for further analysis. To assess the overall trend in annual transcriptome dynamics at the three sites, principal component analysis (PCA) and hierarchical clustering were performed using the 13,385 targets for the 78 samples. PCA was carried out using GeneSpring software. Hierarchical clustering was carried out using the R package “pvclust” (distance: correlation, cluster methods: average, replications: 100) [26].

To detect differences between the sites and identify statistically significant differences in profiles in the time course transcriptome data, the microarray data were analyzed using the R package ‘maSigPro’ [27]. Three between-site comparisons (‘Yamagata vs. Kumamoto’, ‘Yamagata vs. Ibaraki’, and ‘Ibaraki vs. Kumamoto’) were performed using 13 time points of normalized intensity values (cubic regression, false-discovery rate ≤10⁻³, R² ≥0.7). Although a limited number of replicates from each time point were examined (2 replicates), we performed the regression analysis using a sufficient number of time series samples to detect significantly differentially expressed targets.

To identify environmental factors that affect the expression of all 2,505 of the differentially expressed targets among the three comparisons, partial least squares regression (PLSR) analysis was used. In the PLSR analysis, target expression level was used as the response variable, and the meteorological data were used as the predictor variables (S2 Table). As these predictor variables are mutually correlated, we used PLSR to account for the problem of multicollinearity. Meteorological data for the observation sites (Yamagata, observation sites in Higashine-city and Yamagata-city; Ibaraki, observation sites in Hitachi-city; and Kumamoto, observation sites in Kumamoto-city) required to calculate air temperature, precipitation, and sunlight parameters were downloaded from the Japan Meteorological Agency website (https://www.jma.go.jp/jma/indexe.html). To determine the effects of short-term and long-term temperature on the targets, seven parameters of air temperature were calculated across the sampling sites and dates, as follows: air temperature at sampling time (‘temp-0h’); minimum air temperature on the sampling day (‘temp-min’); maximum air temperature on the day before the sampling day (‘temp-max’); and mean air temperature of the previous 7, 30, 60, and 90 days (‘temp-7d’, ‘temp-30d’, ‘temp-60d’, and ‘temp-90d’, respectively). Five parameters of precipitation and five parameters of sunlight were examined: sum of previous 1, 3, 6, 24, and 72 h of precipitation (‘precip-1h’, ‘precip-3h’, ‘precip-6h’, ‘precip-24h’, and ‘precip-72h’, respectively) and sunshine duration (‘sunlight-1h’, ‘sunlight-3h’, ‘sunlight-6h’, ‘sunlight-24h’, and ‘sunlight-72h’, respectively). The day length on the sampling day at the three sites was calculated based on the time of sunrise and sunset obtained from the Ephemeris Computation Office NAOJ website (https://eco.mtk.nao.ac.jp/koyomi/index.html.en). The appropriate number of components was determined based on the explanative performance of the model (R2Y) and predictive performance of the model by cross validation (Q2Y) using the R package ‘ropls’ [28]. Correlation loading plots of the environmental parameters and the targets were created using the appropriate number of components. In addition, Pearson correlation coefficients were calculated for the normalized intensity values of the 2,505 differentially expressed targets and the meteorological data, using data from all three sites (r < −0.7, 0.7 < r, P-value < 0.001).

GO and Pfam enrichment analyses were carried out in order to determine whether the different target lists were enriched for specific GO terms or protein domains. The targets of resulting lists were categorized based on ‘GO biological process complete’ and compared to ‘Arabidopsis thaliana (all genes in database)’ using the PANTHER Overrepresentation Test (released 20220202; http://www.pantherdb.org/). Genes related to cold and UV response were extracted from among the differentially expressed targets according to the GO terms ‘response to cold’ and ‘response to UV’. Prior to the Pfam domain enrichment analysis, the single best ORF of each target was translated into the corresponding amino acid sequence using the tool ‘TransDecoder’, and Pfam domains of each target were identified using the tool ‘rpsblast’ [29–31]. The enrichment P-value for each Pfam domain was calculated using Fisher’s exact test in comparison to the 13,385 targets with a raw signal value exceeding 100 in each pair of samples under at least one of the 39 conditions examined (number of targets > 5, P-value < 0.01).

To further investigate the focused adaptation mechanism, genes related to hormones listed in the ‘plant hormone signal transduction’ pathways of A. thaliana in the Kyoto Encyclopedia of Genes and Genomes (KEGG; https://www.genome.jp/kegg/) and genes related to hormone biosynthesis listed in ‘hormone metabolic pathways and genes in Arabidopsis’ of the RIKEN Plant Hormone Research Network (http://hormones.psc.riken.jp/pathway_hormones.html) were extracted from among the 2,505 significant targets. Genes listed in the ‘starch and sucrose metabolism’ pathway in KEGG and a previous report regarding Arabidopsis [32] were extracted, assuming that those genes are related to starch and sugar metabolism. Similarly, genes related to amino acids were extracted according to the genes listed in the ‘biosynthesis of amino acids’ pathway in KEGG.

Quantitative RT-PCR

To evaluate the reliability of the microarray data, five genes (peroxidase superfamily protein, lipoxygenase5, euonymus lectin S3, cytochrome P450 family 718, and glycosyl hydrolase superfamily protein) that exhibited different seasonal expression patterns in the microarray analysis were analyzed by quantitative RT-PCR (qRT-PCR) using all 78 samples. Primer pairs were designed based on next-generation sequencing isotigs using Oligo version 7.6 (Molecular Biology Insights, Inc., Vondelpark Colorado Springs, CO, USA) and Primer Express version 3.0.1 (Life Technologies) (S3 Table). The first-strand cDNA was synthesized from 500 ng of total RNA using a High Capacity RNA-to-cDNA kit (Life Technologies, Carlsbad, CA, USA). qRT-PCR was performed with Power SYBR Green PCR Master Mix (Life Technologies) and a StepOnePlus Real-time PCR system (Life Technologies), as described in the manufacturer’s instructions. Six microliters of cDNA diluted 1/24 with sterilized water was used in a reaction volume of 20 μl per well. Melting curve analysis was performed from 60 to 95°C, with data captured every 0.3°C to ensure amplification of a single product. Reaction efficiency was assessed using standard curves based on 4-fold dilution series of cDNA synthesized from 2,000 ng of total RNA (1 to 1/256 dilution). Each sample was tested independently and in triplicate using all primers. Transcript abundance was normalized to a putative actin 7 gene, which exhibited a constant microarray value according to the Pfaffl method [33], and the data obtained for each time point were compared with data obtained for shoots collected on February 18 or July 12, 2014, in Ibaraki. Very similar expression patterns were obtained from both the microarray and qRT-PCR analyses, suggesting that the data obtained in this study were reliable (S2 Fig).

Results

Seasonal height growth at the three sites

At Yamagata, Ibaraki, and Kumamoto, the initial tree height (April 2014) was 59.9 cm, 50.1 cm, and 58.6 cm, respectively; the final tree height (December 2014) was 135.2 cm, 130.8 cm, and 146.7 cm, respectively; and the annual growth in 2014 was 75.3 cm, 80.8 cm, and 88.1 cm on average, respectively (S3 Fig). The height growth of the Japanese cedars began in early May in Kumamoto and late May in Yamagata and Ibaraki (Fig 1D). Growth in Ibaraki was more rapid than that in Yamagata during the growth period, from May to October. In Kumamoto, vigorous growth occurred from May to June but was slower in August and September. Height growth ceased in early October in Yamagata and Ibaraki, whereas it continued until late October in Kumamoto. ANOVA of annual growth rate indicated a significant difference in height increase among the three sites (p = 0.05). In Ibaraki and Kumamoto, the cuttings grew approximately 1.6- and 1.5-fold, respectively, during one growing season. The cuttings in Yamagata exhibited less height growth than those in Ibaraki and Kumamoto, with an approximately 1.3-fold growth relative to the initial height.

Annual transcriptome dynamics and site differences

PCA was performed to obtain an overview of the annual transcriptome dynamics at the three sites (Fig 2A). The annual time series data from all sites were plotted along a circle, indicating the seasonal life cycle of Japanese cedar. PC1 explained 69.7% of the total observed variation in gene expression. The PC1 score was high in summer (from June to September) and low in winter (from December to February) (Fig 2A and 2B). PC2 explained 5.9% of the total observed variation in gene expression. The differences in PC2 may be related to the differences in transient seasons between spring (February to April) and autumn (November) (Fig 2C). Between-site variability in annual transcriptome dynamics was also revealed by PCA. The PC1 score was higher in spring (March and April) and autumn (October) at Kumamoto than at Yamagata and Ibaraki, indicating that the transcriptome remained in the growth status for longer at Kumamoto (Fig 2A and 2B). The average PC1 scores in March, April, and October at Kumamoto were 52.3, 34.4, and 54.1 points higher, respectively, than the PC1 scores at Yamagata. The PC2 score at Yamagata was higher compared with Ibaraki and Kumamoto in winter (from December to March), indicating that gene expression differed between the sites, particularly during the dormant period (Fig 2A and 2C). The difference in average PC2 score between Yamagata and Kumamoto was 48.7 during the dormant period (December to March), whereas the difference was 9.1 during the growth period (April to November).

[Figure omitted. See PDF.]

(A) The first two principal components obtained from PCA of annual time series samples of microarray data are shown. The sample names ‘15-Jan’ and ‘15-Feb’ indicate samples collected in January and February 2015, respectively. (B) Relationship between sampling day and PC1 score. (C) Relationship between sampling day and PC2 score.

Hierarchical clustering demonstrated two large clusters, indicating a difference between the growth period (May to October) and dormant period (November to March) (Fig 3). The samples from April divided into the two clusters; the samples from Yamagata and Ibaraki were classified in the cluster with the samples from the dormant period, and the samples from Kumamoto were classified in the cluster with the samples from the growth period. In the cluster containing growth period samples, the samples from Yamagata and Ibaraki collected in October created a sub-cluster, and the samples from Kumamoto created a cluster with other samples from the growth period. In the cluster of dormant-period samples, all six samples from November created a sub-cluster within the dormant-period cluster.

[Figure omitted. See PDF.]

Genes differentially expressed between sites

Analyzing annual expression patterns of the microarray data resulted in the identification of 1,473, 1,137, and 925 differentially expressed targets in the ‘Yamagata vs. Kumamoto’, ‘Yamagata vs. Ibaraki’, and ‘Ibaraki vs. Kumamoto’ comparisons, respectively (Fig 4). The comparisons revealed a total of 2,505 differentially expressed targets (S4 Table). A total of 644 differentially expressed targets were shared between the ‘Yamagata vs. Kumamoto’ and ‘Yamagata vs. Ibaraki’ comparisons, and 308 differentially expressed targets were shared between the ‘Yamagata vs. Kumamoto’ and ‘Ibaraki vs. Kumamoto’ comparisons. In contrast, only 153 differentially expressed targets were shared between the ‘Yamagata vs. Ibaraki’ and ‘Ibaraki vs. Kumamoto’ comparisons. A total of 75 differentially expressed targets were common among all three comparisons (Fig 4, Table 1).

[Figure omitted. See PDF.]

[Figure omitted. See PDF.]

The results of PLSR analyses examining the correlations between the 2,505 differentially expressed targets and meteorological parameters are shown in Fig 5. The correlations between the 2,505 targets (response variables) and latent vectors and the correlations between environmental factors (predictor variables) and latent vectors are presented. Three predictive components were used for the analyses (S5 Table), as the model had high R2Y and Q2X values (0.73 and 0.67, respectively). Parameters closed to targets were positively correlated with the nearby targets, and the meteorological parameters plotted on the opposite side (origin symmetry) were negatively correlated with the targets. As shown in Fig 5, parameters related to air temperature and day length were correlated with most of the analyzed targets. In contrast, few targets were correlated with parameters related to sunlight and precipitation.

[Figure omitted. See PDF.]

Shown is a PLSR correlation-loadings plot of factor 1 versus factor 2, modeled by environmental factors (black rhombi) in the X-matrix and target expression levels (blue dots) in the Y-matrix. The first two PLS components describe 62.1% and 63.7% of the variation in the predictor variables and response variables, respectively.

To identify the most effective meteorological parameter for each target, correlations between the meteorological parameters and expression levels of the 2,505 differentially expressed targets were calculated (S4 Table). Most of the differentially expressed targets showed the highest correlations with air temperature and day length, with 1,842 targets (73.5%) and 222 targets (8.9%), respectively (Table 2). The parameter air temperature at sampling time (‘temp-0h’) correlated with 997 targets. The parameters of long-term air temperature more than 30 days before sampling (‘temp-30d’, ‘temp-60d’, and ‘temp-90d’) showed high correlations with 293, 104, and 35 targets, respectively. The parameters precipitation and sunlight did not show strong correlations with any of the targets. In addition, 441 targets did not show a strong correlation with any of the parameters examined in this study.

[Figure omitted. See PDF.]

GO enrichment analyses revealed that the 2,505 differentially expressed targets were enriched in genes associated with ‘response to stress’ and ‘response to abiotic stimulus’ (Table 3). A total of 571 targets corresponded to 485 Arabidopsis genes associated with ‘response to stress’, and 431 targets corresponded to 381 Arabidopsis genes associated with ‘response to abiotic stimulus’. Among the genes related to ‘response to stress’, 24 targets corresponding to 22 Arabidopsis genes associated with the GO term ‘response to cold’ exhibited a strong correlation with air temperature parameters, including PRP31 and F-BOX PROTEIN 7 (FBP7) (S6 Table). The 11 targets corresponding to 10 Arabidopsis genes related to the GO term ‘response to UV’ exhibited a strong negative correlation with the parameter temperature, including ATCSA-1, root UVB sensitive 1 (RUS1), and UV repair deficient 7 (UVR7) (S7 Table). Genes associated with the GO terms ‘cellular biosynthetic process’ and ‘fatty acid oxidation’ were also enriched among the 2,505 targets. Pfam enrichment analyses indicated that six domains were enriched among the 2,505 targets, including ‘Microtubule binding’ and ‘Polysaccharide biosynthesis protein’ (Table 4). All targets within the ‘Microtubule binding’ domain showed a positive correlation with air temperature (S8 Table).

[Figure omitted. See PDF.]

[Figure omitted. See PDF.]

The 824 targets positively correlated with the parameters of air temperature were enriched in genes associated with ‘microtubule cytoskeleton organization’, ‘plant-type cell wall biogenesis’, ‘developmental growth’, and ‘cell division’, which may be related to organization and growth (Table 5). In contrast, 1,018 targets negatively correlated with air temperature were enriched in genes associated with ‘organic substance metabolic process’, ‘fatty acid beta-oxidation’, and ‘gene expression’. No statistically significant results were obtained for targets correlated with day length.

[Figure omitted. See PDF.]

We extracted 40 plant hormone–related gene targets from the 2,505 differentially expressed targets (Table 6). Five targets were related to abscisic acid (ABA); three of these targets were negatively correlated with ‘temp-0h’, and two of them were negatively correlated with ‘day length’. Four targets were related to ethylene, and all four showed a strong negative correlation with short-term temperature (‘temp-0h’ and ‘temp-max’). Five, six, three, and seven targets were related to growth-related hormones, brassinosteroid, cytokinin, gibberellin, and auxin, respectively. Two brassinosteroid-related genes (brassinosteroid 23-O-glucosilase [UGT73C5], and putative brassinosteroid hydroxylase [CPD/CYP90A1]), two cytokinin-related genes (adenine phosphoribosyl transferase [APT3], and response regulator 1 [RR1]), and four auxin-related genes (tryptophan synthase beta subunit [TSB2], SAUR family protein, transport inhibitor response 1 [TIR1], tryptophan synthase beta subunit homolog [TSB2 homolog]) showed a strong negative correlation with the parameters of air temperature.

[Figure omitted. See PDF.]

Twenty-three targets corresponding to 20 starch- and sugar-related genes were extracted from among the 2,505 differentially expressed targets (Table 7), and these targets showed variations in expression. Expression levels of the starch synthesis–related genes phosphoglucomutase (PGM), AGPase large subunit 2 (APL2), and soluble starch synthase 2 (SS2) were strongly positively correlated with the air temperature parameters ‘temp-0h’ and ‘temp-30d’. Expression levels of the starch degradation–related genes like sex four 1 (LSF1) and beta-amylase 4 (BAM4) were negatively correlated with ‘temp-0h’, and expression levels of the gene glucan water dikinase 1 (GWD1) were negatively correlated with ‘day length’.

[Figure omitted. See PDF.]

Thirty-one targets corresponding to 28 amino acid–related genes were extracted from among the 2,505 differentially expressed targets (Table 8). Of these 31 targets, 19 exhibited a strong correlation with the parameters of air temperature, such as glutamate synthesis–related genes (alanine aminotransferase 2 [ALAAT2] and aspartate aminotransferase 2 [ASP2]), citrulline synthesis–related genes (ornithine carbamoyltransferase [OTC], and peptidase M20/M25/M40 family protein), and a proline synthesis–related gene (pyrroline-5-carboxylate reductase [P5CR]). Among the 31 targets, the gene encoding peptidase M20/M25/M40 family protein exhibited the highest correlation with ‘temp-7d’ (r = −0.97, Fig 6).

[Figure omitted. See PDF.]

[Figure omitted. See PDF.]

Discussion

Annual transcriptome dynamics vary between sites from autumn to spring but not in summer

Our 1-year period microarray data revealed dramatic changes in transcripts over the time course and varied expression patterns between the three sites (Fig 2). Interestingly, differences in the transcriptome dynamics were revealed from November to March. During the dormant period (from January to March), the transcriptome of Japanese cedar planted in the colder site (Yamagata) exhibited a higher PC2 score than the trees planted in the warmer sites (Ibaraki and Kumamoto) (Fig 2C). As annual ring growth of Japanese cedar is positively correlated with temperature in February and March [6], differences in the transcriptome during this season may reflect growth. The transition to the growth period (from March to April) began earlier at the transcriptome level in Kumamoto as compared with Yamagata and Ibaraki, and the transition to the dormant period (from November to December) occurred earlier in Yamagata than in Ibaraki and Kumamoto (Fig 2). This result agreed with the results of hierarchical clustering (Fig 3). The April samples from Kumamoto belonged to the same cluster as growth period samples, whereas the April samples in Yamagata and Ibaraki belonged to the same cluster as the dormant period samples. Also, the October samples from Yamagata and Ibaraki divided from the other growth period samples and created a sub-cluster. These transcriptome differences during the transition period in spring and autumn indicated that the growth period in Kumamoto was longer than that in Yamagata and Ibaraki; thus, the transcriptome differences could have contributed to the longer growth period in Kumamoto (Fig 1D). PCA of the microarray data also demonstrated continuous changes in the transcriptome throughout the dormant period, from November to March, at all three sites. Although Japanese cedar does not grow during the dormant period, these changes may reflect preparations for growth in the following spring.

Only small differences in the transcriptome were observed between the sites during the growth period. The transcriptome of all three sites exhibited higher PC1 scores from June to September, and the data plotted to a similar position (Fig 2B), even though the weekly mean air temperature on the sampling date at the three sites ranged from 19.3 to 27.6°C. Once the growth period began, the Japanese cedar transcriptome may not have been significantly influenced by temperature. In an analysis of Japanese cedar seedlings, Ujino-Ihara [34] reported that heat shock treatment (45°C, 3 h) induced approximately 3,000 differentially expressed unigenes. She also reported that heating pre-treatment reduced the number of differentially expressed genes induced. Thus, these data suggest that the transcriptome changes dramatically when the temperature exceeds a certain threshold and/or that Japanese cedar acclimatize to heat stress during the early stages of the active growth period under field conditions as the temperature gradually rises.

Genes differentially expressed among the three sites

The comparison of ‘Yamagata vs. Kumamoto’ revealed the largest number of significant targets (1,473 targets) among the three comparisons (Fig 4), perhaps because these sites exhibited the largest difference in temperature throughout the year (Fig 1B). The comparisons of ‘Yamagata vs. Ibaraki’ and ‘Ibaraki vs. Kumamoto’ revealed 1,137 and 925 significant targets, respectively. The lower number of significant genes could be attributed to the intermediate temperature in Ibaraki.

The 75 common significant targets among the three comparisons may play an especially important role in environmental responses (Table 1). Among the 75 common targets, the expression of eIF1A showed the highest negative correlation (r = −0.95) with ‘temp-30d’. The mechanisms of translation initiation are fundamentally similar among all organisms, and two eukaryotic translation initiation factors are universally conserved, including eIF1A [35–37]. eIF1A is an important determinant of tolerance to NaCl stress in yeast and plants, and overexpression of eIF1A is associated with a slow-growth phenotype [37]. The higher expression of eIF1A in Japanese cedar at lower temperatures may be related to the response to cold stress and growth cassation at lower temperatures. A positive correlation (r = 0.77) was observed between long-term air temperature (‘temp-90d’) and expression of the gene encoding SAM-Mtase superfamily protein, which is a key enzyme in phenylpropanoid, flavonoid, and other metabolic pathways of biotechnological importance [38]. A negative correlation (r = −0.72) between ‘day length’ and ABA1 expression was also identified (Table 1). The gene encoding highly ABA-induced PP2C protein 3 (HAI3) also exhibited a strong negative correlation with ‘day length’ (Table 6). Expression of ABA-related genes in accordance with changes in day length during xylem formation was also observed in the cambium region of Japanese cedar [39]. These data suggest that the expression of ABA1 is involved in controlling endogenous ABA responses to day length and regulating the timing of xylem formation.

A total of 2,505 targets were differentially expressed in the three comparisons (Fig 4 and S4 Table). GO enrichment analysis indicated that the 2,505 differentially expressed targets included a significantly higher proportion of genes associated with the terms ‘response to stress’ and ‘response to abiotic stimulus’ (Table 3). These genes may play an important role in environmental adaptation. PLSR analysis indicated that the parameters of air temperature and day length subsequently or directly influenced the expression of most of the 2,505 common targets (Fig 5). To estimate the effects of the environment on the expression level of each target, we also calculated Pearson correlation coefficients, and 2,065 targets (82.4%) showed a strong correlation (r > 0.70) with the meteorological parameters examined in this study (Table 2). Air temperature–related parameters were correlated with 73.6% of the 2,505 differentially expressed targets, indicating that temperature is the dominant factor controlling the different transcriptomes between sites. This result was consisted with a study in rice [40, 41], which found that seasonal temperature changes are the dominant factor controlling transcriptome dynamics in the field. Interestingly, long-term air temperature (average air temperature over a 30-day period) was correlated with the expression of 17.2% of the 2,505 differentially expressed targets. These genes may be more affected by long-term temperature trends compared with short-term fluctuations in temperature, which could be more stable. This phenomenon indicates that Japanese cedar may be influenced by temperature trends over the previous 2 to 3 months.

The 441 targets that did not exhibit a strong correlation with meteorological parameters in this study may be influenced by other environmental parameters, such as soil conditions (e.g., temperature, moisture content, and nutrients), or the lack of correlation could perhaps be explained by examining other calculated parameters, such as cumulative temperature values as reported in an analysis of FLC gene expression in Arabidopsis [42]. Although we calculated the correlation to a single environmental parameter in this study, several environmental parameters may affect the expression level of a single target, as reported in rice [40]. Calculations of the effects of individual genes of interest could perhaps provide more details regarding the impact of environmental conditions on gene expression.

Growth-related genes up-regulated as temperature increased

GO enrichment analysis indicated that the 824 targets positively correlated with parameters of air temperature included a significantly higher proportion of genes related to ‘microtubule cytoskeleton organization’, ‘cell wall organization or biogenesis’, and ‘growth’ (Table 5). Pfam enrichment analysis demonstrated a significantly higher proportion of targets including ‘microtubule binding’ domain, and all nine targets associated with ‘microtubule binding’ domain exhibited a strong positive correlation with the parameters of air temperature (Table 4 and S8 Table). These results indicate that Godai1 may promote formation of the microtubule cytoskeleton and growth in summer, especially in warmer sites. Nineteen targets corresponding to 17 Arabidopsis genes were related to ‘microtubule cytoskeleton organization’, including the TON2 gene, which encodes a putative novel protein phosphatase 2A regulatory subunit essential for control of the cortical cytoskeleton in Arabidopsis [43]. Expression of TON2 was positively correlated with ‘temp-60d’; thus, this gene may play an important role in controlling cytoskeleton-related genes. A total of 73 targets corresponding to 65 Arabidopsis genes were related to ‘growth’, including cellulose-related (cellulose synthase A catalytic subunit 1, cellulose synthase interactive 1, and cellulose synthase interactive 3) and expansin-related genes (expansin-A4). This phenomenon suggests that growth differences between sites could be influenced by the length of the growth period.

Acclimation to harsh winter conditions at the colder site

Among the 2,505 differentially expressed targets in this study, 24 targets related to cold response exhibited a strong correlation with the parameters of air temperature, including PRP31 and FBP7 (S6 Table). PRP31 is necessary for pre-mRNA splicing and regulation of the expression of cold-responsive genes [44], and FBP7 is required for protein synthesis during temperature stress in Arabidopsis [45]. These genes may contribute to cold acclimation in winter, especially in colder regions. Our data also indicated the possibility of an important role for chitinases in over-wintering, especially in colder regions. Among the 2,505 differentially expressed targets, a sequence of six targets exhibited high homology to the genes encoding chitinase A, homolog of carrot EP3-3 chitinase, class V chitinase, and chitinase family protein. The expression of these genes exhibited a strong negative correlation with the parameters of air temperature (S4 Table). Plant chitinase has been implicated in defense responses against both biotic and abiotic stressors [46–49]. Diverse chitinases have been identified in spruce, and these enzymes act in concert to protect against freezing injury, store nitrogen, and induce growth cessation by promoting cell maturation [50].

Strong light in low temperature creates an imbalance between light energy absorption and energy use, leading to light stress in winter [51]. Eleven targets related to UV response exhibited a strong negative correlation with the parameters of air temperature, including ATCSA-1, RUS1, and UVR7 (S7 Table). In Arabidopsis, ATCSA-1 plays an important role in UV-B tolerance and genomic integrity [52], and RUS1 and UVR7 mutants are sensitive to UV light [53, 54]. These genes may contribute to photoprotection in winter, especially in colder areas. The photoprotective role of rhodoxanthin in Japanese cedar during cold acclimation is well known [55–58]. Several photoprotective processes that prevent severe damage during exposure to strong light in low-temperature conditions may function in Japanese cedar.

Stress-related hormones may play important roles in acclimating to harsh winter conditions, especially in colder areas. Among the 2,505 differentially expressed targets, five targets related to ABA and four targets related to ethylene showed a strong negative correlation to the parameter of air temperature or day length, indicating up-regulation of the expression of these targets in winter (Table 6). Both ABA and ethylene are known to play roles in regulating bud dormancy [59]. ABA is a major regulator of bud dormancy and promotes starch accumulation in the dormant phase by regulating the expression of starch- and sugar-related genes in grapevine (Vitis vinifera L.) buds [59, 60].

Starch- and sugar-related genes exhibit varied annual expression patterns and differences between sites

Starch breakdown not only enhances cold and frost tolerance, it also contributes energy for bud breaking and shoot growth [61–66]. A previous report indicated that sugar concentration peaks in winter (January) in every part (upper, middle, and lower layer shoots and roots) of the seedlings of Japanese cedar [15]. Frost hardiness is also increased in winter, especially in colder sites [16]. Pfam enrichment analysis indicated a significantly higher proportion of targets associated with ‘Polysaccharide biosynthesis protein’ among the 2,505 differentially expressed targets (Table 4). In addition, 23 targets encoding starch- and sugar-related genes were differentially expressed between sites (Table 7). Among the starch degradation–related genes described in a previous report [32], the BAM4 and LSF1 genes were negatively correlated with ‘temp-0h’, and GWD1 was negatively correlated with ‘day length’. These genes may promote starch degradation at lower temperatures in northern-latitude sites and may contribute to the difference in frost hardiness between sites. In contrast, three target genes related to starch biosynthesis [32] exhibited a positive correlation with parameters related to air temperature (PGM, SS2, and APL2), and these genes may contribute to starch biosynthesis in summer. The varied expression patterns and between-site differences in starch- and sugar-related genes may indicate that Godai1 biosynthesizes sugars and starches that are optimal for the environmental conditions.

Lower temperature may induce energy accumulation in winter for growth in spring

In winter, especially in colder areas, Godai1 may prepare for growth in the upcoming spring by accumulating energy. Brassinosteroid, cytokinin, and auxin are well known as growth-related hormones. A total of eight targets related to these hormones showed a strong negative correlation with the parameters of air temperature in this study (Table 6). We also observed higher expression of genes related to fatty acid beta-oxidation and amino acid synthesis in lower temperatures, which could be related to generation of energy for growth. These results indicate the importance of low temperature to growth in the coming spring in a clone of Japanese cedar.

Beta-oxidation is the main pathway of fatty acid degradation [67]. GO enrichment analysis indicated that the 1,018 targets that were negatively correlated with parameters of air temperature included a significantly higher proportion of genes related to ‘fatty acid beta-oxidation’ (Table 5). This result appears to be consistent with those of previous reports of peanut (Arachis hypogaea L.) [68] and rice (Oryza sativa L.) [69], which exhibited up-regulation of genes related to fatty acid beta-oxidation under cold stress. Up-regulated expression of multiple beta-oxidation genes during germination was also reported in Tung tree (Vernicia fordii, Hemsl.) [70] and Upland cotton (Gossypium hirsutum L.) [71]. As fatty acid beta-oxidation generates energy for a growing embryo [71, 72], it may be an important mechanism to store lipids in winter as an energy source for growth in the coming spring, especially in colder regions.

Amino acids play multiple physiological roles in plants. For example, amino acids function as nitrogen carriers in transport systems, precursors for important metabolites, nitrogen storage molecules, stress response molecules, and signaling molecules [73]. Mori [17] reported that seedlings of Japanese cedar express approximately 20 amino acids, with citrulline, glutamate, and proline being the predominant amino acids in this species. Citrulline, a hydroxyl radical scavenger, accumulates in response to drought stress and nitrogen status in watermelon [74–76]. Mori [77] also reported that citrulline is a major compound involved in nitrogen translocation from roots to the shoot via the xylem in Japanese cedar. Glutamate plays a very important role in plant growth and development and in the response and adaptation to abiotic stress [78–80]. Proline plays a highly beneficial role in plants and also provides the cells sufficient energy to sustain rapid growth [81–83]. With respect to the high level of citrulline and low level of arginine accumulation, Japanese cedar may differ from other gymnosperm tree species that accumulate high levels of arginine [17, 84–86]. Seasonal changes in amino acids have been observed in coniferous trees [17, 87, 88]. In Japanese cedar, citrulline and proline accumulate during the active growth period and then decrease significantly, whereas glutamate exhibits less-marked seasonal changes [17]. In this study, the expression of genes involved in glutamate biosynthesis (ALAAT2 and ASP2) [89], citrulline biosynthesis (OCT and peptidase M20/M25/M40 family protein), and proline biosynthesis (P5CR) were up-regulated as the temperature decreased (Table 8). The gene encoding peptidase M20/M25/M40 family protein demonstrated the highest negative correlation with 7-day mean air temperature among the six glutamate biosynthesis genes, regardless of season and site (Fig 6). This result indicated the importance of low temperature for amino acid biosynthesis. Interestingly, the expression patterns of amino acid biosynthesis genes in this study and the accumulation dynamics of amino acids reported by Mori [17] were not in accordance. The amino acid biosynthesis genes were up-regulated in winter when the temperature decreased, and the amino acids accumulated in the coming spring. It is known that some mRNAs are transcribed during cold acclimation and stored for the subsequent cold de-acclimation [90]. The up-regulation of amino acid biosynthesis–related genes at low temperature (during the dormant period) may represent preparation for growth in the coming spring. The reason for the higher expression of these genes in colder sites remains to be elucidated, however. Godai1 planted in colder regions may accumulate higher amounts of amino acids to promote primary growth in spring. Warmer winters due to climate change could reduce the expression of amino acid–related genes and affect growth in spring. Further study of amino acids in Japanese cedar could increase understanding of the adaptation mechanism during the transition from the dormant to active growth periods.

Conclusions

We revealed changes in the annual transcriptome dynamics in a clone of Japanese cedar by analyzing rooted cuttings planted in three sites differing in terms of climate conditions. We identified a total of 2,505 differentially expressed targets among the three sites that may play important roles in enabling cuttings to adapt to local environmental conditions. The expression of 2,064 targets (82.4%) was affected by air temperature and day length, suggesting that these targets are directly or indirectly (subsequently) regulated by these environmental factors. Although short-term temperature was correlated with most of the genes, long-term temperature was also correlated with some genes, indicating that temperature trends may influence the expression of transcripts over several months. Targets that exhibited a strong negative correlation with air temperature included genes possibly related to cold tolerance and energy accumulation for the coming spring, such as amino acid biosynthesis, starch degradation, and fatty acid beta-oxidation genes. Transcriptome differences between sites observed from autumn to spring may reflect these differentially expressed targets. These results indicate the importance of adaptation to winter climate conditions in planting regions and the importance of low temperature to growth in the coming spring in a clone of Japanese cedar. Understanding how Japanese cedar adapt to various climate conditions is becoming more important due to changes in global climate in recent years. The results of this study may help elucidate the underlying biological mechanisms of environmental response in Japanese cedar.

Geographically, Japanese cedar natural forests range extensively from northern Honshu (30N, 130E) to Yakushima (40N, 140E) [91], and the genome-wide genetic diversity of natural populations and breeding core collections have been reported [92–96]. Tsumura et al. [93] studied the relationships between genotypes of natural populations and environmental variables and identified loci associated with the studied environmental variables. The results of the present study together with those of previous studies suggest the possibility of divergence of environmental responses in Japanese cedar. Although only the Godai1 clone was examined in the present study, other clones may show different environmental responses. Future studies examining additional clones selected from different regions could enhance our understanding of how Japanese cedar adapt to various climate conditions.

Supporting information

S1 Fig. Image of shoot sample, ‘A 10-cm portion of the lateral branch apex’.

https://doi.org/10.1371/journal.pone.0277797.s001

(TIF)

S2 Fig. Expression patterns of five selected targets analyzed by microarray and qRT-PCR.

Blue, green, and pink lines represent the average normalized intensity values in Yamagata, Ibaraki, and Kumamoto.

https://doi.org/10.1371/journal.pone.0277797.s002

(TIF)

S3 Fig.

Annual tree height (A) and annual growth ratio (B) at the three sites.

https://doi.org/10.1371/journal.pone.0277797.s003

(TIF)

S1 Table. Sampling dates in Yamagata, Ibaraki, and Kumamoto.

https://doi.org/10.1371/journal.pone.0277797.s004

(XLSX)

S2 Table. Temperature and day length parameters used in PLSR.

*Sunlight data in Yamagata were obtained at the Yamagata observation site, and other meteorological data were obtained at the Higashine observation site.

https://doi.org/10.1371/journal.pone.0277797.s005

(XLSX)

S3 Table. Primers used in this study for qRT-PCR.

https://doi.org/10.1371/journal.pone.0277797.s006

(XLSX)

S4 Table. The 2,505 significant differentially expressed targets identified by comparisons of the three sites.

https://doi.org/10.1371/journal.pone.0277797.s007

(XLSX)

S5 Table. Three predictive components of PLS analysis.

https://doi.org/10.1371/journal.pone.0277797.s008

(XLSX)

S6 Table. Targets related to the GO term ‘response to cold’ that exhibited a strong negative correlation with parameters of air temperature.

https://doi.org/10.1371/journal.pone.0277797.s009

(XLSX)

S7 Table. Targets related to the GO term ‘response to UV’ that exhibited a strong negative correlation with parameters of air temperature.

https://doi.org/10.1371/journal.pone.0277797.s010

(XLSX)

S8 Table. Targets enriched in the Pfam domain among the 2,505 differentially expressed targets.

https://doi.org/10.1371/journal.pone.0277797.s011

(XLSX)

Acknowledgments

We thank the members of the FTBC, Forestry and Forest Products Research Institute, for assistance with cultivation and greenhouse support of plant materials.

Citation: Nose M, Hanaoka S, Fukatsu E, Kurita M, Miura M, Hiraoka Y, et al. (2023) Changes in annual transcriptome dynamics of a clone of Japanese cedar (Cryptomeria japonica D. Don) planted under different climate conditions. PLoS ONE 18(2): e0277797. https://doi.org/10.1371/journal.pone.0277797

About the Authors:

Mine Nose

Roles: Data curation, Formal analysis, Investigation, Writing – original draft, Writing – review & editing

Affiliation: Forest Tree Breeding Center, Forestry and Forest Products Research Institute, Forest Research and Management Organization, Hitachi, Ibaraki, Japan

ORICD: https://orcid.org/0000-0001-8724-3352

So Hanaoka

Roles: Formal analysis, Investigation

Affiliation: Hokkaido Regional Breeding Office, Forest Tree Breeding Center, Forestry and Forest Products Research Institute, Forest Research and Management Organization, Ebetsu, Hokkaido, Japan

Eitaro Fukatsu

Roles: Formal analysis, Investigation

Affiliation: Forest Tree Breeding Center, Forestry and Forest Products Research Institute, Forest Research and Management Organization, Hitachi, Ibaraki, Japan

Manabu Kurita

Roles: Investigation

Affiliation: Forest Tree Breeding Center, Forestry and Forest Products Research Institute, Forest Research and Management Organization, Hitachi, Ibaraki, Japan

Masahiro Miura

Roles: Investigation, Resources

Affiliation: Kansai Regional Breeding Office, Forest Tree Breeding Center, Forestry and Forest Products Research Institute, Forest Research and Management Organization, Katsuta, Okayama, Japan

Yuichiro Hiraoka

Roles: Investigation, Resources

Affiliation: Faculty of Agricultural Production and Management, Shizuoka Professional University of Agriculture, Iwata, Shizuoka, Japan

Taiichi Iki

Roles: Investigation, Resources

Affiliation: Tohoku Regional Breeding Office, Forest Tree Breeding Center, Forestry and Forest Products Research Institute, Forest Research and Management Organization, Takizawa, Iwate, Japan

Osamu Chigira

Roles: Resources

Affiliation: Iriomote Tropical Tree Breeding Technical Garden, Forest Tree Breeding Center, Forestry and Forest Products Research Institute, Forest Research and Management Organization, Yaeyama-gun, Okinawa, Japan

Kentaro Mishima

Roles: Investigation

Affiliation: Tohoku Regional Breeding Office, Forest Tree Breeding Center, Forestry and Forest Products Research Institute, Forest Research and Management Organization, Takizawa, Iwate, Japan

Makoto Takahashi

Roles: Investigation, Supervision

Affiliation: Forest Tree Breeding Center, Forestry and Forest Products Research Institute, Forest Research and Management Organization, Hitachi, Ibaraki, Japan

Atsushi Watanabe

Roles: Conceptualization, Project administration

E-mail: [email protected]

Affiliation: Department of Forest Environmental Sciences, Faculty of Agriculture, Kyushu University, Nishi-ku, Fukuoka, Japan

References

1. Guan BT, Chung CH, Lin ST, Shen CW. Quantifying height growth and monthly growing degree days relationship of plantation Taiwan spruce. For. Eco. Manage. 2009; 257(11):2270–6.

2. Salminen H, Jalkanen R, Lindholm M. Summer temperature affects the ratio of radial and height growth of Scots pine in northern Finland. Ann. For. Sci. 2009;66(8):810.

3. Kempes CP, West GB, Crowell K, Girvan M. Predicting maximum tree heights and other traits from allometric scaling and resource limitations. PLOS ONE. 2011;6(6):e20551. pmid:21695189

4. Korpela M, Nöjd P, Hollmén J, Mäkinen H, Sulkava M. Hari P. Photosynthesis, temperature and radial growth of Scots pine in northern Finland: identifying the influential time intervals. Trees. 2011;25:323–32.

5. Forestry agency [internet]. the status report of forest resources; 2017 [cited 2021 May 13]. Available from: https://www.rinya.maff.go.jp/j/keikaku/genkyou/h29/3.html. (in Japanese)

6. Takahashi H. Studies on the growth of annual ring of Sugi (Cryptomeria japonica) stands (IV) correlation between the annual ring growth in terms of diameter of some trees and the meteorological factors. J. Jpn. For. Soc. 1971;53(10):315–8.

7. Nishizono T, Kitahara F, Iehara T, Mitsuda Y. Geographical variation in age-height relationships for dominant trees in Japanese cedar (Cryptomeria japonica D. Don) forests in Japan. J. For. Res. 2014;19:305–16.

8. Mori H, Ueno S, Ujino-Ihara T, Fujiwara T, Yamashita K, Kanetani S, et al. Mapping quantitative trait loci for growth and wood property traits in Cryptomeria japonica across multiple environments. Tree Genet. Genomes. 2019;15:43.

9. Mori H, Yamashita K, Saiki S, Matsumoto A, Ujino-Ihara T. Climate sensitivity of Cryptomeria japonica in two contrasting environments: Perspectives from QTL mapping. PLOS ONE. 2020;15(1): e0228278.

10. Kushida T. Studies in dormancy in woody plants. Bull. Mie Univ. For. 2000;24:1–105.

11. Tanino KK, Kalcsits L, Silim S, Kendall E, Gray GR. Temperature-driven plasticity in growth cessation and dormancy development in deciduous woody plants: a working hypothesis suggesting how molecular and cellular function is affected by temperature during dormancy induction. Plant Mol. Biol. 2010;73:49–65. pmid:20191309

12. Cooke JEK, Eriksson ME, Junttila O. The dynamic nature of bud dormancy in trees: environmental control and molecular mechanisms. Plant Cell Environ. 2012;35(10):1707–28. pmid:22670814

13. Itoh T, Shimaji K. Tree growth under controlled environment and its annual rhythm. I Shoot growth in the young stage of Sugi, Karamatsu, Keyaki and Kunugi. Wood Res. Tech. Notes. 1985;20:1–21.

14. Nose M, Kurita M, Tamura M, Matsushita M, Hiraoka Y, Iki T, et al. Effects of day length- and temperature-regulated genes on annual transcriptome dynamics in Japanese cedar (Cryptomeria japonica D. Don), a gymnosperm indeterminate species. PLOS ONE. 2020;15(3):e0229843.

15. Watanabe A. Seasonal changes of carbohydrate in needle and root in Japanese cedar seedlings. Trans. Jpn. For. Soc. 1972;83:238–40. (In Japanese).

16. Horiuchi T, Sakai A. Effect of temperature on the frost hardiness of Cryptomeria japonica. J. Jpn. For Soc. 1973;55:46–51. (In Japanese with English summary)

17. Mori T. Seasonal changes of soluble nitrogen and free amino acids in Cryptomeria japonica seedlings. J. Jpn. Foerst. Soc. 1974;56(10):347–52.

18. Holliday JA, Ralph SG, White R, Bohlmann J, Aitken SN. Global monitoring of autumn gene expression within and among phenotypically divergent populations of Sitka spruce (Picea sitchensis). New Phytol. 2008;178(1):103–122.

19. Hess M, Wildhagen H, Junker LV, Ensminger I. Transcriptome responses to temperature, water availability and photoperiod are conserved among mature trees of two divergent Douglas-fir provenances from a coastal and an interior habitat. BMC Genomics. 2016;17(1):682. pmid:27565139

20. Cronn R, Dolan PC, Jogdeo S, Wegrzyn JL, Neale DB, Clair JBS, et al. Transcription through the eye of a needle: daily and annual cyclic gene expression variation in Douglas-fir needles. BMC Genomics. 2017;18:558. pmid:28738815

21. Mortazavi A, Williams BA, McCue K, Schaeffer L, Wold B. Mapping and quantifying mammalian transcriptomes by RNA-Seq. Nat. Methods. 2008; 5:621–628 pmid:18516045

22. Forest Research and Management Organization [internet]. The characteristics of Sugi plus-tree clones; 2016 [cited 2022 June 8]. Available from: https://www.ffpri.affrc.go.jp/ftbc/business/sinhijnnsyu/seieijyutokuseihyo.html. (in Japanese)

23. Gehrig HH, Winter K, Cushman J, Borland A, Taybi T. An improved RNA isolation method for succulent plant species rich in polyphenols and polysaccharides. Plant Mol. Biol. Rep. 2000;18(4):369–76.

24. Mishima K, Hirao T, Tsubomura M, Tamura M, Kurita M, Nose M, et al. Identification of novel putative causative genes and genetic marker for male sterility in Japanese cedar (Cryptomeria japonica D.Don). BMC Genomics. 2018;19:277.

25. Fukuda Y, Hirao T, Mishima K, Ohira M, Hiraoka Y, Takahashi M, et al. Transcriptome dynamics of rooting zone and aboveground parts of cuttings during adventitious root formation in Cryptomeria japonica D. Don. BMC Plant Biol. 2018;18(1):1–14.

26. Suzuki R, Shimodaira H. Pvclust: an R package for assessing the uncertainty in hierarchical clustering. Bioinformatics. 2006;22(12):1540–1542. pmid:16595560

27. Conesa A, Nueda MJ, Ferrer A, Talon M. maSigPro: a method to identify significantly differential expression profiles in time-course microarray experiments. Bioinform. 2006;22(9):1096–102. pmid:16481333

28. Eriksson I, Johansson E, Kettaneh-Wold N, Wold S. Multi- and megavariate data analysis. Principles and applications, Umetrics Academy, Umeå, 2001.

29. Camacho C, George C, Vahram A, Ning M, Jason P, Bealer K, et al. BLAST+: architecture and applications. BMC Bioinformatics 2009;10:421. pmid:20003500

30. Cock PJA, Chilton JM, Grüning B, Johnson JE, Soranzo N. NCBI BLAST+ integrated into Galaxy. GigaScience.2015;4:39. pmid:26336600

31. Afgan E, Baker D, Batut B, van den Beek M, Bouvier D, Čech M, et al. The Galaxy platform for accessible, reproducible and collaborative biomedical analyses: 2018 update. Nucleic Acids Research. 2018;46:W537–W544. pmid:29790989

32. Streb S, Zeeman SC. Starch metabolism in Arabidopsis. Arabi. Book. 2012;10:e0160.

33. Pfaffl MW. A new mathematical model for relative quantification in real-time RT-PCR. Nucleic Acids Res. 2001;29(9):e45. pmid:11328886

34. Ujino-Ihara T. Transcriptome analysis of heat stressed seedlings with or without pre‐heat treatment in Cryptomeria japonica. Mol. Genet. Genomics. 2020;295:1163–72.

35. Kyrpides NC, Woese CR. Universally conserved translation initiation factors. Proc. Natl. Acad. Sci. 1998;95(1):224–8. pmid:9419357

36. Roll-Mecak A, Shin BS, Dever TE, Burley SK. Engaging the ribosome: universal IFs of translation. Trends Biochem. Sci. 2001;12:705–9. pmid:11738593

37. Rausell A, Kanhonou R, Yenush L, Serrano R, Ros R. The translation initiation factor eIF1A is an important determinant in the tolerance to NaCl stress in yeast and plants. Plant J. 2003;34(3):257–67. pmid:12713533

38. Joshi CP, Chiang VL. Conserved sequence motifs in plant S-adenosyl-L-methionine-dependent methyltransferases. Plant Mol. Biol. 1998;37:663–74. pmid:9687070

39. Mishima K, Fujiwara T, Iki T, Kuroda K, Yamashita K, Tamura M, et al. Transcriptome sequencing and profiling of expressed genes in cambial zone and differentiating xylem of Japanese cedar (Cryptomeria japonica). BMC Genomics. 2014;15:219.

40. Nagano AJ, Sato Y, Mihara M, Antonio BA, Motoyama R, Itoh H, et al. Deciphering and prediction of transcriptome dynamics under fluctuating field conditions. Cell. 2012;10:048. pmid:23217716

41. Nagano AJ, Kawagoe T, Sugisaka J, Honjo MN, Iwayama K, Annual transcriptome dynamics in natural environments reveals plant seasonal adaptation. Nature Plants. 2019;5:74–83. pmid:30617252

42. Aikawa S, Kobayashi MJ, Satake A, Shimizu KK, Kudoh H. Robust control of the seasonal expression of the Arabidopsis FLC gene in a fluctuating environment. PNAS. 2010;107(25):11632–11637. pmid:20534541

43. Camilleri C, Azimzadeh J, Pastuglia M, Bellini C, Grandjean O, Bouchez D. The Arabidopsis TONNEAU2 gene encodes a putative novel protein phosphatase 2A regulatory subunit essential for the control of the cortical cytoskeleton. Plant Cell. 2002;14:833–45.

44. Du JL, Zhang SW, Huang HW, Cai T, Li L. The splicing factor PRP31 is involved in transcriptional gene silencing and stress response in Arabidopsis. Mol. Plant. 2015;8:1053–1068.

45. Calderón-Villalobos LIA, Nill C, Marrocco K, Kretsch T, Schwechheimer C. The evolutionarily conserved Arabidopsis thaliana F-box protein AtFBP7 is required for efficient translation during temperature stress. Gene. 2007;392:106–116. pmid:17240087

46. Gatschet MJ, Taliaferro CM, Porter DR, Anderson MP, Anderson JA, Jackson KW. A cold-regulated protein from bermudagrass crowns is a chitinase. Crop Sci. 1996;36(3):712–18.

47. Yeh S, Moffatt BA, Griffith M, Xiong F, Yang DSC, Wiseman SB, et al. Chitinase genes responsive to cold encode antifreeze proteins in winter cereals. Plant Physiol. 2000;124(3):1251–64. pmid:11080301

48. de los Reyes BG, Taliaferro CM, Anderson MP, Melcher U, McMaugh S. Induced expression of the class II chitinase gene during cold acclimation and dehydration of bermudagrass (Cynodon sp.). Theoretical Appl. Genet. 2001;103(2-3):297–306.

49. Takenaka Y, Nakano S, Tamoi M, Sakuda S, Fukamizo T. Chitinase gene expression in response to environmental stresses in Arabidopsis thaliana: chitinase inhibitor allosamidin enhances stress tolerance. Biosci. Biotech. Biochem. 2009;73:1066–71.

50. González LM, EL Kayal W, Morris JS, Cooke JEK. Diverse chitinases are invoked during the activity-dormancy transition in spruce. Tree Genet. Genom. 2015;11(3):1–21.

51. Öquist G, Huner N. Photosynthesis of overwintering evergreen plants. Annu. Rev. Plant Biol. 2003;54:329–355. pmid:14502994

52. Biedermann S, Hellmann H. The DDB1a interacting proteins ATCSA-1 and DDB2 are critical factors for UV-B tolerance and genomic integrity in Arabidopsis thaliana. 2010;62(3):404–415.

53. Jiang CZ, Yen CN, Cronin K, Mitchell D, Britt AB. UV- and gamma-radiation sensitive mutants of Arabidopsis thaliana. Genetics. 1997;147:1401–1409.

54. Tong H, Leasure CD, Hou X, Yuen G, Briggs W, He ZH. Role of root UV-B sensing in Arabidopsis early seedling development. PNAS. 2008;105(52):21039–21044.

55. Han Q, Mukai Y. Cold acclimation and photoinhibition of photosynthesis accompanied by needle color changes in Cryptomeria japonica during the winter. J. For. Res. 1999;4(3):229–234.

56. Han Q, Shinohara K, Kakubari Y, Mukai Y. Photoprotective role of rhodoxanthin during cold acclimation in Cryptomeria japonica. Plant Cell Environ. 2003;26:715–723.

57. Katahata S, Katoh M, Iio A, Mukai Y. Photoinhibition and pigment composition in relation to needle reddening in sun-exposed Cryptomeria japonica at different altitudes in winter. J. For. Res. 2021;27(2):148–157.

58. Katahata S, Katoh M, Iio A, Mukai Y. Relationship among light intensity, pigment compositions and photoinhibition in relation to rhodoxanthin accumulation in the gymnosperm Cryptomeria japonica. J. For. Res. 2022.

59. Pan W, Liang J, Sui J, Li J, Liu C, Xin Y, et al. ABA and Bud Dormancy in Perennials: Current Knowledge and Future Perspective. Genes. 2021;12:1635. pmid:34681029

60. Rubio S, Noriega X, Pérez FJ. ABA promotes starch synthesis and storage metabolism in dormant grapevine buds. J. Plant Physiol. 2019;234–235:1-8. pmid:30639992

61. Sakai A. The frost-hardening process of woody plant. VI. Seasonal variations in sugars (I). J. Jpn. Forest. Soc. 1960;42(3):97–102.

62. Sakai A. Studies on the frost-hardiness of woody plants; 1.The causal relation between sugar content and frost-hardiness. Contr. Insi. Low Temp. Sci. 1962;B11:1-40.

63. Greer DH, Robinson LA, Hall AJ, Klages K, Donnison H. Frost hardening of Pinus radiata seedlings: effects of temperature on relative growth rate, carbon balance and carbohydrate concentration. Tree Physiol. 2000;20:107–14.

64. Schaberg PG, Snyder MC, Shane JB, Donnelly JR. Seasonal patterns of carbohydrate reserves in red spruce seedlings. Tree Physiol. 2000;20:549–55. pmid:12651436

65. Welling A, Palva ET. Molecular control of cold acclimation in trees. Physiologia Plantarum. 2006;127:167–81.

66. Regier N, Streb S, Zeeman SC, Frey B. Seasonal changes in starch and sugar content of poplar (Populus deltoides × nigra cv. Dorskamp) and the impact of stem girdling on carbohydrate allocation to roots. Tree Physiol. 2010;30:979–87.

67. Goepfert S, Poirier Y. β-Oxidation in fatty acid degradation and beyond. Curr. Opin. Plant Biol. 2007;10(3):245–251.

68. Zhang H, Jiang C, Ren J, Dong J, Shi X, Zhao X, et al. An advanced lipid metabolism system revealed by transcriptomic and lipidomic analyses plays a central role in peanut cold tolerance. Front. Plant Sci. 2020;11:1110. pmid:32849684

69. Liu H, Xin W, Wang Y, Zhang D, Wang J, Zheng H, et al. An integrated analysis of the rice transcriptome and lipidome reveals lipid metabolism plays a central role in rice cold tolerance. BMC Plant Biol. 2022;22:91. pmid:35232394

70. Cui P, Lin Q, Fang D, Zhang L, Li R, Cheng J, et al. Tung tree (Vernicia fordii, Hemsl.) genome and transcriptome sequencing reveals co-ordinate up-regulation of fatty acid β-oxidation and triacylglycerol biosynthesis pathways during eleostearic acid accumulation in seeds. Plant Cell Physiol. 2018;59(10):1990–2003.

71. Sanchez J, Mangat PK, Angeles-Shim RB. Weathering the cold: modifying membrane and storage fatty acid composition of seeds to improve cold germination ability in upland cotton (Gossypium hirsutum L.). Agronomy. 2019;9:684.

72. Hayashi M, Toriyama K, Kondo M, Nishimura M. 2,4-dichlorophenoxybutyric acid-resistant mutants of Arabidopsis have defects in glyoxysomal fatty acid β-oxidation. Plant Cell. 1998;10:183–195.

73. Okumoto S, Funck D, Trovato M, Frolani G. Editorial: Amino acids of the glutamate family: Functions beyond primary metabolism. Front. Plant Sci. 2016;7:318. pmid:27047503

74. Kawasaki S, Miyake C, Kohchi T, Fujii S, Uchida M, Yokota A. Responses of wild watermelon to drought stress: Accumulation of an ArgE homologue and citrulline in leaves during water deficits. Plant Cell Physiol. 2000;41(7):864–73. pmid:10965943

75. Akashi K, Miyake C, Yokota A. Citrulline, a novel compatible solute in drought-tolerant wild watermelon leaves, is an efficient hydroxyl radical scavenger. FEBS Letters. 2001;508:438–42. pmid:11728468

76. Song Q, Joshi M, DiPiazza J, Joshi V. Functional relevance of citrulline in the vegetative tissues of watermelon during abiotic stresses. Front. Plant Sci. 2020;11:512. pmid:32431723

77. Mori T. The role of citrulline in nitrogen translocation of Cryptomeria japonica seedlings. J. Jpn. For. Soc. 1975;57(4):117–20.

78. Walch-Liu P, Liu LH, Remans T, Tester M, Forde BG. Evidence that l -glutamate can act as an exogenous signal to modulate root growth and branching in Arabidopsis thaliana. Plant Cell Physiol. 2006;47(8):1045–57.

79. Yoshida R. Possible amino acid signaling in higher plants. Reg. Plant Growth Develop. 2019;54(1):26–32.

80. Qiu XM, Sun YY, Ye XY, Li ZG. Signaling role of glutamate in plants. Front. Plant Sci. 2020;10:1743. pmid:32063909

81. Nanjo T, Kobayashi M, Yoshiba Y, Sanada Y, Wada K, Tsukaya H, et al. Biological functions of proline in morphogenesis and osmotolerance revealed in antisense transgenic Arabidopsis thaliana. Plant J. 1999;18(2):185–93.

82. Mattioli R, Costantino P, Trovato M. Proline accumulation in plants. Not only stress. Plant Signaling Behavior. 2009;4(11):1016–18.

83. Hayat S, Hayat Q, Alyemeni MN, Wani AS, Pichtel J, Ahmad A. Role of proline under changing environments. A review. Plant Signaling Behavior. 2012;7(11):1456–66.

84. Hatano K. On the amino acids in the seeds and seedlings of Ginkgo biloba. Tokyo Univ. Forests Bull. 1953;44:221–227.

85. Durzan DJ, Steward FC. The nitrogen metabolism of Picea glauca (moench) voss and pinus banksiana lamb. as influenced by mineral nutrition. Can. J. Bot. 1967;45:695–710.

86. Durzan DJ. Nitrogen metabolism of Picea glauca (I). Seasonal changes of free amino acids in buds, shoot apices, and leaves, and the metabolism of uniformly labelled 14C-L-arginine by buds during the onset of dormancy. Can. J. Bot. 1968;46:909–19.

87. Sagisaka S. Plant metabolisms in cold conditions. Kagaku to Seibutsu 1974;12(1):15–29. (In Japanese)

88. Näsholm T, Ericsson A. Seasonal changes in amino acids, protein and total nitrogen in needles of fertilized Scots pine trees. Tree Physiol. 1990;6(3):267–81. pmid:14972938

89. Diab H, Limami AM. Reconfiguration of N metabolism upon hypoxia stress and recovery: Roles of Alanine Aminotransferase (AlaAT) and Glutamate Dehydrogenase (GDH). Plants. 2016;5(2):25. pmid:27258319

90. Nakaminami K, Matsui A, Nakagami H, Minami A, Nomura Y, Tanaka M, et al. Analysis of differential expression patterns of mRNA and protein during cold-acclimation and de-acclimation in Arabidopsis. Mol. Cell. Proteomics. 2014;13(12):3602–11.

91. Hayashi Y. Taxonomical and phytogeographical study of Japanese conifers. Norin-Shuppan, Tokyo. 1960 (in Japanese).

92. Tsumura Y, Kado T, Takahashi T, Tani N, Ujino-Ihara T, Iwata H. Genome scan to detect genetic structure and adaptive genes of natural populations of Cryptomeria japonica. Genetics 2007;176:2393–2403.

93. Tsumura Y, Uchiyama K, Moriguchi Y, Ueno S, Ihara-Ujino T. Genome scanning for detecting adaptive genes along environmental gradients in the Japanese conifer, Cryptomeria japonica. Heredity 2012;109:349–360.

94. Tsumura Y, Uchiyama K, Moriguchi Y, Kimura MK, Ueno S, Ujino-Ihara T. Genetic differentiation and evolutionary adaptation in Cryptomeria japonica. Genes Genomes Genetics 2014;4(12):2389–2402.

95. Uchiyama K, Ujino-Ihara T, Ueno S, Taguchi Y, Futamura N, Shinohara K, et al. Single nucleotide polymorphisms in Cryptomeria japonica: their discovery and validation for genome mapping and diversity studies. Tree Genetics & Genomes. 2012;8:1213–1222.

96. Uchiyama K, Miyamoto N, Takahashi M, Watanabe A, Tsumura Y. Population genetic structure and the effect of historical human activity on the genetic variability of Cryptomeria japonica core collection, in Japan. Tree Genetics & Genomes. 2014;10:1257–1270.