Norway spruce deploys tissue ‐ specific responses during acclimation to cold

Climate change in the conifer ‐ dominated boreal forest is expected to lead to warmer but more dynamic winter air temperatures, reducing the depth and duration of snow cover and lowering winter soil temperatures. To gain insight into the mechanisms that have enabled conifers to dominate extreme cold environments, we performed genome ‐ wide RNA ‐ Seq analysis from needles and roots of non ‐ dormant two ‐ year Norway spruce ( Picea abies (L.) H. Karst), and contrasted these response to herbaceous model Arabidopsis We show that the main transcriptional response of Norway spruce needles exposed to cold was delayed relative to Arabidopsis, and this delay was associated with slower development of freezing tolerance. Despite this difference in timing, Norway spruce principally utilizes early response transcription factors (TFs) belonging to the same gene families as Arabidopsis, indicating broad evolutionary conservation of cold response networks. In keeping with their different metabolic and developmental states, needles and root of Norway spruce showed contrasting results. Regulatory network analysis identified both conserved TFs with known roles in cold acclimation (e.g. homologs of ICE1 , AKS3 , and of the NAC and AP2/ERF superfamilies), but also a root ‐ specific bHLH101 homolog, providing functional insights into cold stress response strategies in Norway spruce.

Investigations directed at improving frost tolerance in agricultural plants often target the CBF pathway because of its wide conservation in plants (Thomashow, 2010).However, in recent genome-wide transcriptomic analyses, other early cold-induced transcription factors and their role in regulating cold-responsive (COR) genes have received attention, highlighting the need for whole-genome transcriptional analyses to unravel the complexity of the low-temperature gene regulatory network (Park et al., 2015).
In woody perennials, overwintering is initiated after sensing the change of season (Welling et al., 2004), with the shortening of the photoperiod in late summer/early autumn inducing growth cessation, the development of dormancy and cold hardening (Bigras et al., 2001;Chang et al., 2015;Cooke et al., 2012;Guy, 1990;Li et al., 2004;Rossi et al., 2008;Welling et al., 2004).However, maximal frost tolerance is only acquired after exposure to temperatures below 0°C (Beck et al., 2004;Bigras et al., 2001;Greer & Warrington, 1982;Sakai, 1966;Sogaard et al., 2009;Weiser, 1970).The establishment of frost tolerance to extremely low temperatures (<−60°C) in perennial or overwintering tissues makes it possible for above-ground parts to survive above the snow cover and may be a result of greater cellular dehydration than found in herbaceous species (Lang et al., 1994;Rinne et al., 1998;Welling et al., 1997Welling et al., , 2004)), coupled to more intensive accumulation of cryoprotective compounds (Coleman et al., 1991;Kuroda & Sagisaka, 1993;Rinne et al., 1994;Sauter & Wellenkamp, 1998).High concentrations of solutes also increase intracellular viscosity, stabilizing cells when stressed and leading to the formation of aqueous glasses in woody plants (Wisniewski et al., 2003).Nevertheless, woody perennials and herbaceous plant species also share mechanisms of cold regulation and cold-regulated target genes.For example, CBFs have been found in deciduous (poplar, birch) (Benedict et al., 2006;Nanjo et al., 2004;Welling & Palva, 2008) and evergreen angiosperm tree species (eucalyptus) (El Kayal et al., 2006;Navarro et al., 2009) and are involved in the activation of the cold responses of both leaves and winter dormant tissue (Benedict et al., 2006;Welling & Palva, 2008), where specialization after perennial-driven evolution might explain differences in the transcriptomes.
Boreal forests cover about 11% of the earth's surface (Bonan & Shugart, 1989) and are dominated by evergreen conifers of the genera Abies, Picea and Pinus, which can survive extended periods of temperatures below -40°C when fully cold-acclimated (Sakai & Weiser, 1973;Sakai, 1966;Strimbeck et al., 2007Strimbeck et al., , 2008)).In contrast to deciduous trees of temperate and boreal regions, evergreen conifers such as Norway spruce (Picea abies (L.) H. Karst) maintain their photosynthetic tissues for several years.To keep the evergreen foliage alive throughout the winter, the needles have to acquire extreme low-temperature tolerance (Strimbeck et al., 2007(Strimbeck et al., , 2008)).
Similar to angiosperms, changes in carbohydrates (Strimbeck et al., 2008), accumulation of low-molecular-weight cryoprotectant metabolites (Chang et al., 2015;Crosatti et al., 2013) and cryoprotective proteins such as dehydrins (Kjellsen et al., 2013), contribute to acquired freezing tolerance in conifer needles.Frost injuries can also occur in belowground tissues of plants, with soil frost causing fine root dieback, reducing nutrient and water uptake by trees (Groffman et al., 2001).Normally, snow cover insulates the soil from the cold air temperatures and reduces freeze-thaw events in the soil (Campbell et al., 2005), protecting roots from extreme temperature changes.
Although over-wintering conifer roots have been shown to cold acclimate (Burr et al., 2001;Flint et al., 1967) the mechanisms used by roots to develop cold tolerance have received little attention.Arabidopsis has been reported to show as little as 14% overlap in coldinduced leaf and root transcriptomes, based on 8 K microarray chips (Kreps et al., 2002).Furthermore, the limited data available suggests that conifer roots remain metabolically active longer than aboveground tissues and, while they do cold acclimate, they do not develop the same deep frost tolerance as the above-ground tissues (Bigras et al., 2001), suggesting that the molecular responses of roots to changing seasonal temperatures warrant closer investigation.
Climate change models predict that temperatures will increase and become more variable at higher latitudes, particularly during the winter months (Christensen et al., 2007).These changes in seasonal temperatures will not only increase the length of the growing season (Barichivich et al., 2013) but possibly also delay the onset of cold acclimation, impair the development of freezing tolerance in the autumn and lead to early deacclimation during the late winter (Chang et al., 2016;Frechette et al., 2016;Guak et al., 1998;Repo et al., 1996;Stinziano et al., 2015;Wang, 1996).Furthermore, current evidence also suggests that the risk of belowground frost injury is increasing (Campbell et al., 2005), with warmer winter air temperatures reducing the depth and duration of the snow cover, leading to colder soil temperatures (Decker et al., 2003;Groffman et al., 2001).
Conifers have evolved separately from angiosperms for more than 300 million years (Bowe et al., 2000) and large-scale whole-genome transcriptional profiling experiments are needed to understand whether the same gene regulatory networks are involved in cold acclimation in conifers.The assembly of a draft genome of Norway spruce (Nystedt et al., 2013) has made this species an ideal coniferous model.To gain insight into the mechanisms that have enabled conifers to dominate the boreal forest under current climatic conditions, we performed genome-wide RNA-Seq analysis from needles and roots of nondormant 2-year-old Norway spruce seedlings exposed to cold (5°C) and freezing (−5°C) temperatures.This experimental design allowed us to focus only on the temperature responses, without influence from growth cessation and dormancyrelated mechanisms.For direct comparison with the standard herbaceous model, Arabidopsis leaves exposed to cold (5°C) were sampled at equivalent time points.

| Experimental design and sampling
Two-year-old Norway spruce (Picea abies (L.) Karst.)seedlings of the seed provenance Lilla Istad (56°30ʹN) were planted in peat in 3 L pots and grown at 18°C/15°C light/dark in a 16-h light period.
Plants had been potted up from 1 L pots 10 weeks before the start of the experiment so that the roots had explored the soil volume.
At the beginning of the experiment, 2 h into the light period, 40 seedlings were shifted to controlled environment rooms with constant 5°C and the 16-h light period was maintained to ensure that only temperature was changed.Needles and roots were collected immediately before transfer and then 6 h, 24 h, 3 days and 10 days after the start of the cold treatment.After this 10-day cold acclimation period at 5°C, the seedlings were transferred to subzero freezing temperature at −5°C (16-h light period) and samples were again collected after 6 h, 24 h, 3 days and 10 days.
Samples were always taken from previously unsampled plants.
Needles were sampled from the most recently expanded mature shoots.Actively growing fine roots (>1 mm in diameter) were sampled, rinsed with distilled water, and blotted dry.Every sampling time point was replicated by needles from five seedlings, three of which were also sampled for roots.Samples were collected onto dry ice and then stored at −80°C until further processing.
Seeds of Arabidopsis thaliana (Col-0) were sown into soil and kept in an 8 h photoperiod (150 μmol photon m −2 s −1 ) at 23°C.After 14 days seedlings were transplanted into individual pots.After a further 30 days plants were shifted to 5°C, 2 h into the photoperiod and the most recently expanded mature leaves of these plants were sampled after 0 h, 6 h, 24 h, 3 days and 10 days at 5°C.Single leaves were harvested from 5 to 10 randomly chosen plants and were pooled.Three pools were collected per time point.

| RNA preparation and sequencing
Norway spruce samples were prepared by the CTAB method (Chang et al., 1993) with the following modifications: Addition of warm extraction buffer, including polyvinyl pyrrolidinone (PVP) 40, and vortex mixing of the ground sample material was followed by an incubation step for 5 min at 65°C.Precipitation with ¼ volume 10 M LiCl took place at −20°C for 2 h and RNA was then harvested by centrifugation at 14 000 rpm for 20 min and 4°C.The RNA was further purified using the RNeasy mini kit (Qiagen).A DNase Digestion with the RNase-free DNase set (Qiagen) was included in the procedure.Highquality total RNA with an RIN ≥ 7.5, OD 260/280 ratio of ≥2.0 and concentrations ≥50 ng/µl were sequenced by SciLifeLab for pairedend (2 × 125 bps).The sequencing library preparation included enrichment for poly-adenylated messenger RNAs and all samples yielded >9.4 million read pairs.Arabidopsis thaliana (Col-0) leaf samples were isolated using the Plant RNA Mini kit (E.Z.N.A.) according to the manufacturer's instructions.DNase treatment was performed after RNA extraction using DNA-free™ DNA Removal Kit (Ambion, Life Technologies).
Samples were sequenced by BGI (Beijing Genomics Institute) and all yielded >23 million paired-end reads.
All RNA samples' integrity was analysed by Agilent RNA 6000 Nano kit (Agilent Technologies) on a Bioanalyzer 2100 (Agilent Technologies) and purity was measured with a NanoDrop 2000 spectrophotometer (Nanodrop Technologies).All samples that passed quality controls were sequenced by Illumina HiSeq.2000 platform.

| Freezing tolerance, cold index injury and H 2 O 2 measurements
To analyse the freezing tolerance of Arabidopsis and Norway spruce tissues we followed the protocol described by Strand et al. (2003).Briefly, for Arabidopsis, two leaf discs (1 cm 2 ) from the most recently fully expanded leaves of 6-week-old plants (nonacclimated as well as 3-and 10-days 5°C cold-acclimated) were put into individual glass tubes containing 200 μl high-performance liquid chromatography (HPLC)-grade double deionized water.Two-year-old Norway spruce plants growing at 18°C were cold acclimated at 5°C for 3 and 10 days and needles and roots samples were collected.Additionally, seedlings that had been cold acclimated for 10 days at 5°C were subsequently exposed to freezing at -5°C for 3 days after which needles and roots were sampled.
For the freezing tolerance assay, five needles from twigs of four seedlings were collected and cut into 10 5-mm pieces and placed in glass vials with 1 ml of HPLC-grade double deionized water (4°C).Additionally, fine roots from four seedlings were collected and three sections (10 mm) were placed in glass vials with 500 μl of HPLC-grade double deionized water.All the root samples were rinsed with 4°C deionized water before being placed in the vials.
All samples were kept on ice before moving to the freezing bath.
For both species, the samples were incubated in a programmable bath (Julabo FP45) first at -2°C for 1 h, at which point ice formation was initiated using a metal probe cooled by liquid nitrogen.
The samples were then subjected to decreasing temperatures at a cooling rate of −2°C/h.At designated intervals, the samples were removed from the freezing bath and kept on ice at 4°C overnight.
On the following day, 2 and 1.3 ml of double deionized water were added to Norway spruce and Arabidopsis sample vials, respectively, and placed on a shaker at room temperature (RT) overnight.
Electrolyte leakage was measured, followed by flash freezing the samples in liquid Nitrogen and shaking again overnight at RT. Final conductivity was then determined.For both species, ion leakage measurements were made using a conductivity metre (CDM210; Radiometer).
For the comparison of different Norway spruce tissues, freezing injury was expressed as an index of injury.Index of injury accounts for the ion leakage present in the tissues aside from that caused by freezing treatment by converting the percentage release of electrolytes to a scale from 0 (unfrozen samples) to 100 (damaged samples).Roots contain higher amounts of residual electrolytes than needles due to the soil around them and the usage of the index of injury to express the freezing damage addresses this issue (Burr et al., 2001;Flint et al., 1967).
Freezing injury of Norway spruce tissues was expressed as index of injury and calculated using the following equation (Driessche, 1976): where L i is the initial conductivity of a sample, L f is the final conductivity of a sample, B i is the Initial conductivity of a blank and B f is the final conductivity of a blank.
For H 2 O 2 measurements, 2-year-old Norway spruce seedlings grown and subjected to the same treatments used for the RNA extraction were sampled and analysed using Amplex ® Red Hydrogen Peroxide/Peroxidase Assay Kit (Molecular Probes).Tissues were harvested and ground by using liquid Nitrogen.For each tissue four biological replicates were analysed.

| Data processing
Sequence read preprocessing and quality assessment of the raw data was performed as in Delhomme et al. (2014) following the standard guidelines using FastQC-0.10.1 https://www.bioinformatics.babraham.

| Functional analysis
The predicted TF genes and their associated families in Norway spruce were obtained from PlantTFDB 3.0 (Jin et al., 2014).Norway spruce sequences and GO annotations were obtained from Conifer Genome Integrative Explorer (ConGenIE; http://congenie.org) (Nystedt et al., 2013).Gene Ontology and GOSlim enrichment analyses were performed with a standalone version of GeneMerge (Castillo-Davis & Hartl, 2003) using as background population a transcriptome of 43 398 genes expressed in the current data set.
To compare globally the functions against Arabidopsis thaliana the GOSlim tags were assigned to Norway spruce genes using the plant slims subset available from the Gene Ontology Consortium website (Ashburner et al., 2000;Beike et al., 2015) using the Map2Slim function from Owltools https://github.com/owlcollab/owltools/wiki/Map2Slim.To identify previously characterized COR genes in Norway Spruce, we performed a bidirectional BlastP (Altschul et al., 1990) using Arabidopsis thaliana (TAIR 10) (Garcia-Hernandez et al., 2002;Rhee et al., 2003) and Picea abies v1.0 references proteomes (Nystedt et al., 2013), which are available in PlantGenIE (Sundell et al., 2015).COR genes were obtained from Park et al. (2015).The best BLAST matches were used to assign gene aliases.Peptide sequences of 10 or fewer amino acids were removed.A total amount of 58 587 amino acid sequences were analysed in Norway spruce after excluding low-confidence sequences.Orthology analysis was performed using the PLAZA resource (Van Bel et al., 2018).

| Regulatory network analysis to find pivotal TF of cold stress response
The upstream regions of Norway spruce gene models containing start and stop codons were obtained from the FTP website available in ConGenIE and dropped all of them at the same lengths of 1 Kb using a Perl script.Genes with upstream regions shorter than 1 Kb were excluded from the analysis.In total 37 621 upstream regions were analysed.Counts of the different analysed motifs in the promoters regions (Table 1) were obtained with the stand-alone version of PatMatch (Yan et al., 2005).A motif was considered overrepresented in a target gene promoter if was present in an upstream region more than the respective confidence interval (CI) upper bound of the abundance of the different analysed motifs (Chawade et al., 2007).The analysed motifs were consensus target sequences for ERF, NAC, MYB, bZIP, bHLH, AP2, DRE, LBD and WRKY families.
The upper limit of the (CIs was below one, with an exception for bHLH where the upper limit was over four.In these cases, just a single occurrence was considered as a case of over-representation because these sequences are not sporadically distributed on the promoter regions.For bHLH five or more times was considered as over-representation of the motif in a promoter.We got the coexpression values between all the identified COR genes and identified  approximately equal numbers up-and down-regulated (Figure 1b).
Comparing the homology relationship of the differentially expressed genes (DEGs) between species using orthology information from the Gymno PLAZA resource (Van Bel et al., 2018) (Figure 1c) it was found that of the 1949 Arabidopsis cold-induced DEGs, 340 had an identified ortholog that also responded to cold in Norway spruce (validated orthologs).A further 47 had a potential ortholog, representing homologous genes identified by best match BLAST (Pearson, 2013), which were also differentially regulated by cold in Norway spruce but not classified as orthologs by PLAZA (Figure 1c and Table S1).1395 (72%) of the Arabidopsis DEGs had identified homologs in Norway spruce that were not cold-responsive in Norway spruce.An additional 167 cold-induced Arabidopsis DEGs had no identifiable homolog in Norway spruce (Arabidopsis singletons).
Similarly, of the 3251 cold-induced DEGs in Norway spruce, 401 had orthologs in Arabidopsis and an additional 74 had potential Arabidopsis orthologs.There were 2068 (64%) Norway spruce coldinduced DEGs with homologs that were not cold-responsive in Arabidopsis and another 708 that had no identifiable homolog in Arabidopsis (Norway spruce singletons).These data revealed a core cold-induced regulon of 387 orthologous genes in Arabidopsis and 475 orthologous COR genes in Norway spruce, which represented 20% and 15% of the cold-induced DEGs, respectively.The orthologous genes were associated with the Gene Ontology (GO) categories 'stress response', 'ion binding' and 'nucleic acid binding transcription factor activity', consistent with the known importance of these processes in cold acclimation and demonstrating conservation of the transcriptional response between these two divergent species (Tables S2 and S3).
To compare the functional response of the two species in more detail, we performed DEG-enrichment analysis using GOSlim categories.We found that the number of enriched categories increased over time in both species, indicating an extensive and progressive remodulation of the transcriptome with the progression of cold exposure (Figure 2 and Tables S4 and S5).Both species showed significant enrichment of genes involved in 'response to stress', 'nucleic acid binding transcription factor activity', 'plasma membrane' and 'carbohydrate metabolic processes'.A distinct response in Arabidopsis resulted in significant enrichment of up-regulated genes in the 'lipid metabolic process', 'secondary metabolic process', 'cell wall organization or biogenesis' and 'catabolic process' that were not significantly enriched in Norway spruce.On the contrary, only Norway spruce showed significant enrichment for 'vacuole', 'transport', 'transmembrane transporter activity', 'ion binding', 'cellular protein modification process', 'signal transduction' and 'Golgi apparatus' and these processes became significant only after 10 days.DEGs within these categories were also induced in Arabidopsis but did not become significantly enriched at any time (Figure 2 and Tables S4   and S5).Among down-regulated genes, Arabidopsis showed strong  Bowden, 1985).S1 were analysed (Figure 3).First, transcriptional regulators known to respond to cold in Arabidopsis were analysed and showed an induction at 6 h with a fast decline observed by 24 h (Figure 3a), as shown previously (Chinnusamy et al., 2003;Kim et al., 2013;Kim et al., 2017;Shen et al., 2015;Vogel et al., 2005).Overall, all of the 120 up-regulated TF-DECs in Arabidopsis responded to cold between 6 and 24 h at 5°C (Figure 3b, S1, and S2) and most belong to ERF, bHLH and MYB families, which were previously reported as families containing members involved in plant cold and dehydration responses (Chinnusamy et al., 2003;Fowler & Thomashow, 2002;Vogel et al., 2005;Yamaguchi-Shinozaki & Shinozaki, 2006).These Arabidopsis TF-DEC fell into two main clusters; one, including ERF6 and MYB7, with a transient expression peak at 6 h (the top cluster in Figure 3b); a second larger cluster with a broader expression peak in which expression was induced by 6 h and remained high until 24 h after exposure to cold (bottom cluster, Figure 3b).The second cluster included several AP2/ERF superfamily members, such as TINY2, DEAR2, DEAR4, ERF3 and AP2 (Figure 3b).This contrasted sharply with what was observed for the 107 up-regulated TF-DEC from Norway spruce, where only a few TFs were induced early (at 6 or 24 h) and these were six ORA47 homologs, one AP2, one MYB33, one ERF53, one NAC025 and one TINY2-like homolog (Figure 3b).The main response of TFs in Norway spruce was delayed relative to Arabidopsis and formed two subclusters, one showing induction between 24 h and 3 days and a second subcluster with an even later expression maximum at 10 days that was completely absent in the Arabidopsis response (Figure 3b).This delayed response in Norway spruce was associated with TFs that belong to AP2/ERF superfamily such as ERF1, ERF2, ERF9 and TEM1 homologs, along with other TFs such as ZFP4; MYB123 and MYB101; anac028 and anac078 and WRKY20-like genes whose role in the cold stress response has not been directly reported in Arabidopsis (Figure 3b).This delayed response in cold-induced gene expression was also reflected in the slower acquisition of frost tolerance by Norway spruce needles relative to Arabidopsis leaves (Figure S12).Thus, both the freezing tolerance and related gene expression data indicated that while Norway spruce shared an early transcriptional response with Arabidopsis at 6 to 24 h, the more extremophile conifer mounted a later transcriptional response following 3 and 10 days of exposure to cold.

| Characterizing Norway spruce responses in needles and roots to cold and freezing
The effect of a progressive exposure to cold (5°C) on the transcriptomes of developmentally and metabolically distinct needles and roots of Norway spruce revealed prominent transcriptional changes in both tissues after 3-10 days at 5°C (Figure 4).Although the root response lagged behind at 6 h, presumably due to the effect of soil acting as a temperature buffer, after 24 h the number of DEGs was largely similar between the two tissues.Subsequent exposure of both tissues to freezing at −5°C had no noticeable effect on the total number of DEGs in needles but did result in a further increase in root DEGs after 3 and 10 days at −5°C, giving a total of 4324 and 4407 up-regulated DEGs in needles and roots respectively (Figure 4b).This progressive increase in DEGs was reflected in increased freezing tolerance in both tissues at 3 and 10 days at 5°C.However, following subsequent exposure to −5°C for 3 days, both tissues showed evidence of freeze injury indicative of only partial cold acclimation, where needles showed higher injury than roots (Figure S13).
GOSlim enrichment analysis showed that while both needles and roots deployed similar numbers of up-regulated DEGs and most GOSlim categories were common between the tissues, in needles, more categories were significantly enriched (Figure 4c).For example, genes in the 'transport' category were enriched in up-regulated DEGs common to both tissues and also in needle-specific genes, showing that modulating the transport of solutes and ions across membranes is an important response to cold in both tissues, but the response by needles is stronger and more diverse.This stronger response by needles to exposure to cold is further demonstrated by the finding that 'plasma membrane', 'transmembrane transporter activity' and 'ion binding' categories were only enriched in the up-regulated needle-specific DEGs.
Unsurprisingly, in the down-regulated DEGs photosynthesis and plastid categories were the main gene classes significantly enriched in the needle-specific list.In addition, down-regulation of 'nucleic acid binding transcription factor activity' genes at 6 and 24 h under freezing, reflected deactivation of bHLH, bZIP, MYB, NAC and AP2/ ERF superfamily TFs in needles following exposure to freezing temperatures (Figure 4c and Table S6 and S7).On the contrary, rootspecific and common down-regulated DEGs were enriched in 'plasma membrane', 'kinase activity' and 'signal transduction' categories, indicating that the silencing of some components of these categories occurred in both tissues but that this suppression response was more pronounced in roots exposed to subzero temperatures.indicating those that were Root-specific (eight clusters) (Figure 5b).

|
Correlation and cluster analysis was then performed with all identified clusters and seven main responses were found, termed Super Clusters (SC) (Figure 5c and Table S9).SC-6 identified a needlespecific behaviour enriched in the GO term 'response to chitin', characterized by early induction that peaked at 6 h, followed by expression falling below control levels by 10 days of cold.SC-7 also showed early response behaviour but was composed of DEGs from both tissues, and reached an expression peak after 24 h of cold.This cluster contained GO enrichments in 'response to biotic stimulus', 'response to wounding', 'protein serine/threonine phosphatase activity' and in terms related to regulation of transcription (Table S9).
On the contrary, expression of SC-1 showed a progressive increase under cold starting at 6 h, and showed a more pronounced increase in expression following exposure to freezing.Similarly, SC-2 and SC-3 showed a progressive increase in gene expression until 10 days at 5°C, but displayed a flat profile under freezing.These three SCs (SC-1, SC-2 and SC-3), represent gene expression profiles from both tissues and included GO enrichments related to oxidation-reduction process, response to biotic stress, transport and plasma membrane (Table S9).SC-4 was identified as a root-specific response that showed a progressive increase in gene expression until 10 days at 5°C, and then a continuous drop during freezing.This cluster showed GO enrichments in 'cysteine-type peptidase activity', 'response to nonionic osmotic stress', 'response to lipid hydroperoxide', 'apoplast', 'pyridoxine biosynthetic process' (vitamin B6) and unexpectedly in 'chlorophyll metabolic process' (Table S9).Thus, early and late response clusters were identified in both tissues, but with maximum induction levels at different time points.In addition to this different timing, while the needle-specific responses identified were composed of genes with biotic and abiotic stress annotations, the root-specific responses were enriched for genes related to protein degradation, vitamin B6 and redox metabolism, possibly reflecting the different metabolism of these contrasting tissues and demonstrating the need to independently investigate the stress responses of these different tissues.

| Identification and tissue specificity of transcription factors significantly induced by cold
The transcription factors associated with these different cold acclimation responses in Norway spruce needles and roots were examined and 215 TFs significantly induced in response to cold (TF-DEC) were identified from needles and 183 from roots, of which 107 (37%) were needle-specific, 75 (26%) root-specific and 108 (37%) common to both tissues (Figure 6 and Table S10).In addition, 20% (57 genes) of these Norway spruce genes lacked orthologs in Arabidopsis (Figure S4).
F I G U R E 5 Identification of coordinated responses to cold.(a) Venn diagram representing needle-specific, root-specific and common up-regulated differentially expressed genes (DEGs).A Gene Ontology (GO) enrichment analysis of these gene lists is available in Table S8.(b) Heat maps and main clusters of tissue-specific and common up-regulated DEGs.Common DEGs were analysed from both tissues and each identified cluster (C) was separated according to cluster trends in needles (N) or in roots (R).(c) Super Clusters (SC) were defined by comparing clusters from (b) using Pearson correlation analysis.Common clusters (C) were separated by tissues to compare against tissue-specific clusters.Thus, a cluster trend for a common cluster number 'i' was called 'C i in needles' or 'C i in roots' according to the tissue analysed.For each SC, scaled gene expression mean values are represented with dotted lines.Grey areas represent variability by two standard deviations.An analysis of SC distribution in the network is available in Figure S3.A GO enrichment analysis for these SC is available in Table S9 [Color figure can be viewed at wileyonlinelibrary.com]I G U R E 6 Transcription factors differentially regulated by cold (TF-DEC).TF-DEC were analysed by heat maps and hierarchical clustering (corrected p ≤ 0.01 and fold change ≥2) crossing normalized needle and root gene expression data.Family members were obtained from the Plant Transcription Factor Database (Jin et al., 2014) and normalized VST expressions were scaled by row means.A file including gene descriptions for each cluster is available in Table S10 [Color figure can be viewed at wileyonlinelibrary.com]Needle-specific TF-DEC formed two main clusters (Figure 6); cluster II TFs were induced between 6 h and 3 days after exposure to cold and returned to pre-cold exposure expression levels or lower by Day 10, cluster I was comprised of TFs induced from 3 to 10 days after exposure to cold.No freezing-specific induced TFs were found in the needle-specific group.Of the root-specific group, cluster IV was comprised of TFs induced early by cold while cluster III was comprised of TFs induced after 10 days of cold and reaching induction peaks under freezing.Unique to roots, clusters V and VI were strongly induced only following exposure to freezing at −5°C (Figure 6 and Table S10).Within the group of significant induced TFs common to both tissues, cluster X was comprised of a group of ERF TFs that were induced sharply in needles at 6-24 h in the cold and then returned to control levels.
However, these same genes were markedly induced in roots only after freezing.Cluster VIII included TFs that reached an induction peak in needles at 10 days in the cold and thereafter expression levels remained high during freezing.However, in roots members of this cluster were induced much earlier, from 24 h exposure to cold, and stayed at high levels under freezing.Cluster IX and XI showed stronger inductions in roots than in needles.Thus, even though there were a large number of common TFs significantly induced in both tissues by exposure to cold and freezing, the majority of these showed differential regulation between the two tissues, either in timing or the strength of their induction and therefore when this is added to the evidence of clusters of tissuespecific TFs, it demonstrates that these dissimilar tissues display very different transcriptional control during acclimation to cold.

| Cold response regulatory network of Norway spruce
To analyse in detail how COR genes were regulated and to identify TF 'hubs' playing a key role in COR gene regulation, we performed a regulatory network analysis of TF-DEC and COR genes combining coexpression and promoter motif analysis of all COR genes.A network model of inferred transcriptional regulations of COR genes was built using an approach previously used to identify putative plant cold acclimation pathways and key TFs in Arabidopsis and rice (Chawade et al., 2007;Lindlof et al., 2009).The analysed motifs were consensus target sequences for ERF, NAC, MYB, bZIP, bHLH, AP2, DRE, LBD and WRKY families, which have previously been associated with plant cold stress (Chawade et al., 2007;Peng et al., 2015).For each motif, (CI were obtained (Table 1) and used as a criterion of overrepresentation.Thus, regulatory interaction between two nodes (genes) was the result of a combination between TF-target coexpression and the over-representation of the recognized motif for the respective TF in the target COR gene promoter.Using a strict correlation threshold, the resulting regulatory network included 2135 links of TF and target genes and 910 nodes (genes).784 were COR target genes and 126 were TFs, with mainly bHLH, MYB, NAC, and ERF family members (Figure 7a).Analysis of the connectivity of the nodes (degree) in the network showed that the network followed a scale-free behaviour (Figure 7b), with only a few genes highly connected.These node degree measurements were then used as a gene prioritization criterion (Moreau & Tranchevent, 2012) and the top 10 most connected nodes (hubs) were analysed, all of which were TFs (Figure 7c).The most connected hub corresponded to the gene MA_448849g0010; a bHLH TF (ICE1-like) that putatively regulated (on the basis of TF-COR co-expression and TF binding motif presence) 159 COR genes.Interestingly, the third most-connected hub (MA_68586g0010) was a bHLH101-like that was only differentially regulated in roots with 88% of the target genes having root-specific expression (Figure 7c and Table S13).These results demonstrated that the cold regulatory network in Norway spruce is highly interconnected, that most of the cold response circuits were common to both tissues and that these include homologs of TFs previously reported as regulators of stress response genes in other species such as Arabidopsis (Figure S5 and Table S13).A notable exception to this was the root hub bHLH101 homolog (MA_68586g0010), which suggests that tissue-specific responses are present in this species, although the function of this hub and its downstream genes remain to be elucidated.

| DISCUSSION
A qualitative comparative analysis of Norway spruce needles with leaves of the temperate herbaceous model plant Arabidopsis showed that, despite their long evolutionary separation, both species share a core of orthologous genes that respond to cold (Figure 1), suggesting that this core response predated the split of the lineages.This is in common with the presence of a conserved, core drought response (Haas et al., 2021).The orthologous responses include the upregulation of several TFs of the ERF3 family, including TCP2 and HB13, that are known to play a role in plant stress responses (Cabello et al., 2012;Chawade et al., 2007;Liu et al., 1998;Nakano et al., 2006;Stockinger et al., 1997) (Table S3).In addition, construction of a gene regulatory network demonstrated that an ICE1like TF was, like its counterpart in Arabidopsis, a central player in mediating the cold response in needles and roots of Norway spruce.
However, Norway spruce demonstrated different transcriptional dynamics relative to Arabidopsis, and this contrasted delayed response was associated with the induction of a large cohort of TFs, including ERF1, ERF3, anac028 and anac078 homologs, that have not previously been described to have a function in cold acclimation.In addition, a bHLH101-like TF was shown to be co-expressed with a root-specific subset of genes in the gene-regulatory network.No central regulators have previously been identified as root-specific and this finding indicates that tissue-specific responses are important in the cold hardiness response, at least in perennial species such as conifers .This is in agreement with similar findings for drought response in Norway spruce (Haas et al., 2021).temperatures (Figure S6) but they must also be protected from freezing damage and from anoxia stemming from ice encasement.To date, metabolomic and proteomic analyses have shown that, as for herbaceous leaves, carbohydrate and lipid metabolism and the accumulation of dehydrins play key roles in survival at extremely low temperatures in perennial evergreen needles (Angelcheva et al., 2014;Kjellsen et al., 2013;Strimbeck et al., 2015) and that both needles and roots increased cold acclimation in response to exposure to low  temperatures, as measured by electrolyte leakage (Burr et al., 2001;Flint et al., 1967).In our analysis, we found that both tissues mounted a progressive transcriptional response to cold and that following subsequent freezing, most common cold acclimation responses were maintained in both tissues.GO enrichment analysis performed on all DEGs that were induced in Norway spruce revealed that regardless of the tissue, common genes induced in response to cold were overrepresented for genes within the core stress categories of 'regulation of transcription', 'transport' and 'response to wounding' (Figure 5b and Table S8).However, the responses of the two tissues did vary (Figure 4c) and GO enrichment analysis showed that needle-specific induced genes were enriched for 'plasma membrane', 'ATP-binding', 'DNA binding', 'integral component of membrane' and 'Golgi apparatus'; while root-specific induced DEGs were enriched for the categories 'peroxidase activity', 'oxidation-reduction process', 'metal ion binding' and 'structural constituent of ribosome' and 'translation' (Table S8).These contrasting responses indicated that while needles are exposed to extreme cold and desiccation, demanding membrane reorganization and protection; roots, which remain more metabolically active during winter (Law & Anthoni, 1999;Martz et al., 2016), may be under higher oxidative stress.Consistent with this conclusion, the Arabidopsis ortholog of the novel bHLH101-like TF identified in Norway spruce roots has previously been reported to be involved in iron homeostasis (H.Y. Wang et al., 2007) and photo-oxidative stress responses (Noshi et al., 2018).While we cannot yet assign any specific function for this TF in roots, the root-specific transient accumulation of H 2 O 2 (Figure S11) coupled with our evidence that the first-degree neighbourhood of this TF is enriched in oxidoreductase activity genes (Table S12), suggests it may play a central role in the oxidative stress response of roots to cold stress.However, further studies will be required to establish the role of this TF and its root-associated regulon.
The best studies of the cold response pathways in herbaceous species is that regulated by the CBF/DREB1 (CRT-binding factor/ DRE-binding protein) family of TFs, which belong to the ERF family (Cook et al., 2004;Fowler & Thomashow, 2002;Jin et al., 2014;Shinozaki & Yamaguchi-Shinozaki, 1996).In Arabidopsis, the CBFs show a transient early induction by cold in roots and leaves (Hruz et al., 2008;Kilian et al., 2007).A single ortholog of CBF1 and CBF3 was found in Norway spruce (MA_20987g0010), which was induced between 6 and 24 h at 5°C in roots, although it did not pass the statistical filters to be classified as an up-regulated DEG or COR gene and it was not responsive to cold in needles (Figure S7).Cluster X of the common TF-DEC (Figure 6 et al., 2003;Ding et al., 2015;Kim et al., 2015;Tang et al., 2020).
These data demonstrate that early response TFs and COR genes can be induced by ICE1 homologs under cold in several plant species such as Zea mays (Lu et al., 2017), rice (Zhang et al., 2017), tomato (Feng et al., 2013), Vitis amurensis (Xu et al., 2014), Pyrus ussuriensis (Huang, Li, et al., 2015) and Poncirus trifoliata (Huang, Zhang, et al., 2015), in line with our finding that in Norway spruce this ICE1-like protein is an important potential regulator of many COR genes, where 150 of these genes appear as targets of ICE1 in our regulatory network (Table S13), and whose progressive late induction suggests that its role may become even stronger following freezing (Figure S10).The recent finding that an ICE2 homolog has been associated with local adaptions to shade and cold at northern latitudes in Norway spruce (Ranade & García-Gil, 2021) further supports a key role for ICE-like proteins in cold tolerance in Norway spruce.
In conclusion, our study provides a comprehensive overview of the transcriptome of Norway spruce involved in cold stress and reveals potential regulators of COR.The TF-DECs and the hubs identified from our regulatory network can be used to identify the most suitable candidate genes for targeted genetic modifications or for directed breeding to generate high-yielding cold stress-tolerant trees to aid the forest industry as it confronts new environmental challenges resulting from climate change.
TF from Norway spruce using Spearman correlations values, performing a significance test of the correlation coefficients and FDRcorrected p values were obtained for them.A threshold of Spearman correlation ≥0.8 and adjusted p ≤ 0.01 was used to filter the correlations.A regulatory interaction in our network was obtained by the combination of a co-expression data and over-representation of binding motifs.So, if a pair TF-target gene overcome a correlation threshold and the target gene has an over-representation of the motif, which recognize the TF in the promoter, then these TF-target gene pair make up a regulatory interaction in our regulatory network.Pajek software was used to obtain network parameters and visualizations(de Nooy et al., 2005).

3. 1 |
Comparing the cold acclimation transcriptional response of Norway spruce needles and Arabidopsis leaves RNA-Seq analysis of needles and leaves exposed to cold (5°C) over a time course of 10 days revealed that the bulk of the transcriptional response of Norway spruce needles was slower than for Arabidopsis leaves (Figure1a,b).Overall, both species differentially expressed approximately 15% of all transcribed genes in response to cold, with 'M' = A or C, 'N' = any nucleotide A or T or C or G, Y = C or T. THE NORWAY SPRUCE COLD TRANSCRIPTOME | 431 enrichment for genes related to 'structural constituent of ribosome', 'translation' and 'RNA binding' that were all largely absent in the Norway spruce response, suggesting a strong reorganization of the translational machinery in the herbaceous leaf that was not present in the Norway spruce needle.Similarly, Arabidopsis strongly and rapidly down-regulated 'plastid' and 'photosynthesis' genes but Norway spruce did so much less and much later (Figure 2).To identify the transcription factors (TFs) that drove the transcriptional responses in both species, the composition and expression profiles of TFs Differentially Expressed in response to Cold (TF-DEC) Comparing the response of Arabidopsis thaliana leaf and Picea abies (Norway spruce) needle transcriptomes exposed to 5°C.(a) Hierarchical clustering using normalized data (see methods).The red numbers correspond to approximately unbiased (AU) values and the green ones to bootstrap probability (BP) values.(b) Analysis of transcriptome progression in response to cold.Differentially expressed gene lists (DEGs) were obtained at each point in the time series, compared against the control, and then represented as a percentage of the transcriptome.DEGs significantly induced in Arabidopsis (green) and Spruce (dark green) and significantly repressed DEGs in Arabidopsis (yellow) and Spruce (dark yellow) were obtained by filtering the data by corrected p ≤ 0.01 and fold change ≥2.(c) Orthologs, homologs and species-specific DEGs for both species (down and up-regulated).Validated orthologs correspond to orthologous genes that are differentially regulated by cold in both species.Gene lists for each group and functional information are available in Table

F
I G U R E 2 GOSlim functional analysis of species comparison.Aerial tissue transcriptomic responses of Arabidopsis thaliana (ATH) and Picea abies (Spruce) to chilling (5°C) were compared functionally by GOSlim analysis.Differentially expressed gene (DEGs) lists were analysed along the time series treatments with different GOSlim tags assigned in the up and down-regulated DEGs.The number of genes in each category is represented by bubble sizes (counts).GOSlim enrichments with Bonferroni corrected p ≤ 0.05 are represented with a red colour scale.Nonsignificant p values (corrected p > 0.05) are in grey.Here only a selected group of GOSlim categories are represented.A full list including all the categories and p values is included in TableS4and a description of the genes and their homologs in TableS5[Color figure can be viewed at wileyonlinelibrary.com]U R E 3 (See caption on next page) Analysis of common and tissue-specific responses to cold Despite being developmentally and metabolically distinct, needles and roots shared 2002 up-regulated DEGs, with a further 2322 being needle-specific and 2405 root-specific (Figure 5a).To identify coordinated positive responses to cold and contrast the tissue-specific responses, we performed hierarchical cluster analysis on the upregulated DEGs.The identified clusters were designated as N prefix clusters for Needle-specific (seven clusters), C prefix clusters for those Common to both tissues (11 clusters) and R prefix clusters

FF
Transcription factor analysis.(a) Expression profiles of previously characterized transcription factors (TF) involved in the cold response in Arabidopsis thaliana: CAMTA4 (AT1G67310), ICE1 (AT3G26744), CBF1 (AT4G25490), CBF2 (AT4G25470), CBF3 (AT4G25480) and ZAT12 (AT5G59820) are shown from our experiment (n = 3) using variance stabilizing transformation (VST) gene expression values.Errors bars represent SD.(b) TF differentially expressed by cold (TF-DEC) were analysed in both Arabidopsis and Norway spruce.TF with positive changes relative to control are shown (corrected p ≤ 0.01 and fold change ≥2).VST data were scaled by row means.For each heatmap, zoom versions including all the identifiers are available in Figures S1 and S2 [Color figure can be viewed at wileyonlinelibrary.com]I G U R E 4 The response of needle and root transcriptomes from Norway spruce to exposure to cold (5°C) and freezing (−5°C).(a) Hierarchical clustering of gene expression data from needle and root samples (see Section 2).Red numbers correspond to approximately unbiased (AU) values and green ones to Bootstrap Probability (BP) values.(b) Analysis of transcriptome progression in response to cold (5°C) and freezing (−5°C).Differentially expressed gene lists (DEGs) were obtained at each point in the time series, compared against the control.DEGs in needles (dark green represents induced genes and light green represents repressed genes) and roots (brown represents induced genes and tan represents repressed genes) were obtained by filtering the data by corrected p ≤ 0.01 and fold change ≥2.(c) The number of DEGs (counts) with different GOSlim tags assigned are represented by bubble size.GOSlim enrichments with Bonferroni corrected p ≤ 0.05 are represented with a red colour scale.Nonsignificant p values (corrected p > 0.05) are in grey.A full list of GOSlim enrichments is included in Table S6 and a complete description of each GOSlim category and the gene functions in Table S7 [Color figure can be viewed at wileyonlinelibrary.com]I G U R E 5 (See caption on next page) THE NORWAY SPRUCE COLD TRANSCRIPTOME | 437 Evergreen plants such as Norway spruce maintain their needles above the snowpack during winter and thus require mechanisms to protect the needles from extreme low temperatures and the associated desiccation.On the contrary, roots face less extreme THE NORWAY SPRUCE COLD TRANSCRIPTOME | 439 Regulatory network analysis.(a) Network representation of predicted regulatory interactions between transcription factors (TF) and cold-responsive (COR) genes.TFs are represented by diamonds and their family by colours.COR genes are represented by circles and coloured according to the tissue in which they are differentially regulated.(b) Network Degree distribution in Log10/Log10 scale.(c) Subnetwork of the 10 genes with the highest centrality.Gene Ontology (GO) enrichments in the hub neighbourhoods are available in ) contained nine ERF TFs showing early induction in needles, indicating that other genes of the same TF family have a role in cold stress regulation in aerial tissues of Norway spruce and that this stress response may be supported by an expansion of the ERF family of TFs in Norway spruce (FiguresS8 and S9).On the basis of co-expression and cis-element overrepresentation analysis, we built a regulatory network to identify potential key conserved TFs mediating cold acclimation in Norway spruce (Figure7).We found extensive crosstalk between COR genes and TFs and a hierarchical structure in their connections, being COR genes regulated by upstream TFs and also a module with a root-specific bias.Interestingly, the most connected TF in the network was an ICE1-like gene (MA_448849g0010), orthologous to the Arabidopsis bHLH transcription factor INDUCER OF CBF EXPRES-SION1 (ICE1).Although recently a direct role of ICE1 in the regulation of CBFs has been questioned by data showing that an ice1-1 mutation does not impair cold induction of the CBF genes(Kidokoro et al., 2020), other studies provide evidence that ICE1 does play a role in cold response transcriptional regulation by binding to promoters of CBFs following activation by the kinase OST1 (Chinnusamy