ABSTRACT: Cold acclimation in conifers is a complex process, the timing and extent of which reflects local adaptation and varies widely along latitudinal gradients for many temperate and boreal tree species. In spite of their ecological and economic importance, little is known about the global changes in gene expression that accompany autumn cold acclimation in conifers. Using three populations of Sitka spruce (Picea sitchensis) spanning the species range, and a Picea cDNA microarray with 21,840 unique elements, we monitored within and among-population gene expression during the fall. Microarray data were validated for selected genes using real-time PCR. Similar numbers of genes were significantly two-fold upregulated (1,257) and downregulated (967) between late summer and early winter. Among those upregulated were dehydrins, pathogenesis-related/antifreeze genes, carbohydrate and lipid metabolism genes, and genes involved in signal transduction and transcriptional regulation. Among-population microarray hybridizations at early and late autumn time points revealed substantial variation in the autumn transcriptome, some of which may reflect local adaptation. Our results demonstrate the complexity of cold acclimation in conifers, highlight similarities and differences to cold tolerance in annual plants, and provide a solid foundation for functional and genetic studies of this important adaptive process in conifers. Keywords: Time course Foliage for RNA extraction was obtained from four-year-old Sitka spruce seedlings, which were grown from seed collected across the species range in a raised-bed outdoor common garden in Vancouver, BC, Canada (49º N) (Mimura and Aitken, 2007, Heredity). Needle tissues from current year upper lateral shoots, separated from stem and bark tissue, were collected at five time points between late summer and early winter, 2004. Three populations were chosen for sampling: Valdez, Alaska, USA (61º N) (AK), Prince Rupert, British Columbia, Canada (54º N) (BC) and Redwood, California, USA (41º N) (CA) (Figure 1a). Needle samples (~1g) were taken from eight individuals in the BC population on each of the five dates (Aug. 30, Oct. 18, Nov. 24, Dec. 1, and Dec. 13). Eight individuals from the AK and CA populations were also sampled on the second and fourth dates, as well as eight additional individuals from the BC population. Tissues were flash frozen in an N2 vapor tank immediately upon collection, and subsequently stored at -80°C until processing. To decrease the effects of biological variance among individual seedlings within populations, equal amounts of foliar tissue were pooled from four individuals prior to RNA extraction. Two pools were collected at each time point for each population, and total RNA was extracted following a previously published protocol (Kolosova et al., 2004). Slides were scanned and spot intensity quantified using ImaGene software (BioDiscovery, Inc., El Segundo, CA). To correct for background intensity, the lowest 10% of median foreground intensities was subtracted from the median foreground intensities. Data were then normalized by variance stabilizing normalization (VSN) to compensate for nonlinearity of intensity distributions (Huber et al., 2002). To identify significant changes in gene expression, a linear mixed-effects model was fit to the normalized intensities in the Cy3 and Cy5 channels of the 32 microarray slides. The model contained an adjustment for dye bias, an array effect indicating which Cy5/Cy3 pair was on each array, a treatment effect indicating sample population and time point, and a random effect to adjust for repeated measures on the same biological sample (Kerr et al., 2000). P-values were computed for each gene-by-treatment effect and Q-values were calculated to adjust for the false discovery rate (FDR) (Storey & Tibshirani, 2003).
