Constitutive and Cold Acclimation-Regulated Protein Expression Profiles of Scots Pine Seedlings Reveal Potential for Adaptive Capacity of Geographically Distant Populations

Geographically distant Scots pine (Pinus sylvestris L.) populations are adapted to specific photoperiods and temperature gradients, and markedly vary in the timing of growth patterns and adaptive traits. To understand the variability of adaptive capacity within species, molecular mechanisms that govern the physiological aspects of phenotypic plasticity should be addressed. Protein expression analysis is capable of depicting molecular events closely linked to phenotype formation. Therefore, in this study, we used comparative proteomics analysis to differentiate Scots pine genotypes originating from geographically distant populations in Europe, which show distinct growth and cold adaptation phenotypes. Needles were collected from 3-month-old seedlings originating from populations in Spain, Lithuania and Finland. Under active growth-promoting conditions and upon acclimation treatment, 65 and 53 differentially expressed proteins were identified, respectively. Constitutive protein expression differences detected during active growth were associated with cell metabolism and stress response, and conveyed a population-specific adaptation to the distinct climatic conditions. Acclimation-induced protein expression patterns suggested the presence of a similar cold adaptation mechanism among the populations. Variation of adaptive capacity among the genotypes was potentially represented by a constitutive low level of expression of the Ser/Thr-protein phosphatase, the negative regulator of the adaptive response. Also, overall less pronounced acclimation-induced response in seedlings from the Spanish population was observed. Thus, our study demonstrates that comparative proteomic analysis of young conifer seedlings is capable of providing insights into adaptation processes at the cellular level, which could help to infer variability of adaptive capacity within the plant species.


Introduction
Sustainability of future forests is threatened by anthropogenic impacts, mainly due to pollution and climate change [1,2]. Therefore, special attention is being given to understanding of the molecular basis of adaptive traits in tree species (e.g., [3,4]). High genetic diversity, phenotypic plasticity and epigenetic mechanisms are the key factors of the adaptive capacity of forest tree populations (e.g., [5,6]). Trees have a low evolutionary rate in terms of nucleotide substitution and diversification rates as compared to herbaceous species [7]. On the other hand, tree populations are larger, less genetically structured and contain larger genetic diversity [7,8]. Most forest tree species show substantial variability in adaptive traits at levels much higher than is observed at neutral genetic markers and form a collection of highly differentiated populations with contrasting adaptation to local climates [9]. Intra-specific genetic variation of tree species is linked to the high adaptability under constantly changing environments [3] providing diversity and plasticity of genetic variants for evolution, where phenotypic plasticity is a property of a genotype to produce different phenotypes in response to particular environmental conditions that vary both temporally and geographically. Among the main evolutionary forces influencing the variability patterns, natural selection primarily acts to increase genetic differentiation among populations, while gene flow and phenotypic plasticity tend to reduce it [10]. In wide-spread wind pollinated conifers such as Scots pine, the effects of genetic drift and mutations are weak. Although natural variation in epigenetic marks and the relation to phenotypic traits in trees is still an under-explored area, epigenetic variation has been suggested to contribute to the phenotypic plasticity and adaptive potential of forest trees [11].
Gene expression is a dynamic process and depends on the specific tissue, developmental stage and environmental signals. Gene expression analysis at both transcriptome and proteome level could provide details about population variability, and the mechanisms involved in evolution and adaptation [12,13]. Protein expression is the final step in gene expression process, thus, proteomic analysis bridges the gap between the genotype and phenotype, providing a depiction of molecular events that are most directly linked to the phenotype at cellular or organism level. Quantitative variation of protein abundance is the result of a complex network of regulatory mechanisms that integrate genetic variation and interaction with the environment of individual plants [14,15]. A study on Arabidopsis ecotypes showed that the proteomics approach could constitute a powerful tool to mine the biodiversity between ecotypes of a single plant species [16]. Proteomics was used to study phenotypic and molecular response to water deficit of two Eucalyptus genotypes [17] and adaptive variation in Picea abies (L.) H.Karst. ecotypes [18,19].
Scots pine (Pinus sylvestris L.) is the most widespread coniferous tree species occupying an extensive area in Eurasia. The complex biogeographic history of Scots pine in Europe has been well studied [20][21][22]. At the neutral part of the genome, in open pollinated wide-spread conifers such as Scots pine, the phylogenic structure is weak and population differentiation is low (e.g., [23,24]). Based on the maternally inherited mitochondrial DNA markers, the European range of Scots pine is grouped into three major haplotypes: south western (refugia in Pyrenees), central (refugia in Balkan and Italy) and south eastern (refugium in the southern edge of the Ural mountains) [20,[25][26][27]. Gene flow and isolation by distance were other major factors affecting population differentiation in the Scots pine that were defined using neutral DNA markers [28]. Numerous studies based on indoor and field trials revealed a significant clinal variation among the Scots pine populations over latitude gradients in adaptive traits such as phenology [29,30].
Timing of cold acclimation is of adaptive significance for northerly conifers such as the Scots pine [31][32][33]. The seedlings cease active growth and enter the endodormancy stage in response to elongating nights towards the autumn [34], followed by shoot hardiness development in response to chilling temperatures in the early autumn [35,36]. In trees, cold acclimation treatment results in complex structural, biochemical and genetic shifts from elongating meristems to completely dormant and frost hardy tissues [37]. In brief, it causes protoplasmic dehydration, concentration of solutes, reduction of osmotic potential, upregulation of abscisic acid, a change in carbohydrate metabolism along with changes in membrane permeability and fatty acid composition [38,39]. The most prominent changes are induced in the chloroplasts and manifest as structural reorganization of the organelar membranes [40,41]. At protein level, the increased amount of antioxidant enzymes [42][43][44] and accumulation of dehydrins mostly localized in chloroplasts and the mitochondria membrane system was reported during the acclimation [45]. Gene expression studies showed upregulation of enzymes involved with synthesis of sugars and sugar derivatives, dehydration stress and other stress tolerance enzymes, such as heat-shock proteins [46].
Owing to the photoperiod and temperature-triggered adaptation in the northern regions, conifers have population-specific requirement for critical night length to budset and chilling temperatures to induce variable levels of shoot frost hardiness in early autumn [36]. When transferred to a single location, the critical night length for bud set and shoot hardiness development vary markedly between geographically distant populations within a latitudinal gradient [29,32,47]. Therefore, after a cold acclimation treatment, markedly different gene expression profiles could be expected in seedlings from the geographically distant Scots pine populations.
To probe an adaptive capacity within the Scots pine species, our study assessed phenotypic traits and quantitative and qualitative protein expression patterns in the Scots pine seedlings originating from the geographically distant populations of northern and southern Europe. Protein expression patterns were assessed by two-dimensional protein electrophoresis. Firstly, differences among the three populations were assessed during the active growth period. Furthermore, an acclimation treatment of the seedlings was used to extend the analysis of protein expression profiles with specific emphasis on cold adaptation processes in the distinct Scots pine genotypes.

Seedling Production and Cold Acclimation Treatment
We used seeds of Scots pine collected from the three autochthonous populations representing the Mediterranean, temporal and boreal climatic zones in Europe: central Spain (Valsain Forest mountains, Sierra de Guadarrama province, 40 • Figure 1). The seeds were collected in native Scots pine stands naturally-born over generations in large forest tracts, where the influence of non-autochthonous sources is negligible. These populations mainly represent a geographically large latitudinal gradient of genepools adapted to distinct photoperiodic and temperate requirements to initiate dormancy and cold hardiness processes that have a major effect on timing of cold acclimation of the seedlings. The seed lots were mixtures from more than 20 trees within a forest stand, which was considered as representative of a population. Owing to the photoperiod and temperature-triggered adaptation in the northern regions, conifers have population-specific requirement for critical night length to budset and chilling temperatures to induce variable levels of shoot frost hardiness in early autumn [36]. When transferred to a single location, the critical night length for bud set and shoot hardiness development vary markedly between geographically distant populations within a latitudinal gradient [29,32,47]. Therefore, after a cold acclimation treatment, markedly different gene expression profiles could be expected in seedlings from the geographically distant Scots pine populations.
To probe an adaptive capacity within the Scots pine species, our study assessed phenotypic traits and quantitative and qualitative protein expression patterns in the Scots pine seedlings originating from the geographically distant populations of northern and southern Europe. Protein expression patterns were assessed by two-dimensional protein electrophoresis. Firstly, differences among the three populations were assessed during the active growth period. Furthermore, an acclimation treatment of the seedlings was used to extend the analysis of protein expression profiles with specific emphasis on cold adaptation processes in the distinct Scots pine genotypes.

Seedling Production and Cold Acclimation Treatment
We used seeds of Scots pine collected from the three autochthonous populations representing the Mediterranean, temporal and boreal climatic zones in Europe: central Spain (Valsain Forest mountains, Sierra de Guadarrama province, 40°49′ lat., 3°57′ long., 1828 m.a.s.l. alt.), central Lithuania (Kazlu Ruda forest tact, 54°46′ lat., 23°32′ long., 79 m.a.s.l. alt.) and northern Finland (Suomussalmi forest, 65°07′ lat., 28°49′ long., 233 m.a.s.l. alt .), respectively ( Figure 1). The seeds were collected in native Scots pine stands naturally-born over generations in large forest tracts, where the influence of non-autochthonous sources is negligible. These populations mainly represent a geographically large latitudinal gradient of genepools adapted to distinct photoperiodic and temperate requirements to initiate dormancy and cold hardiness processes that have a major effect on timing of cold acclimation of the seedlings. The seed lots were mixtures from more than 20 trees within a forest stand, which was considered as representative of a population.  Seeds were soaked overnight in sterile deionized water and surface sterilized for 30 min with 30% hydrogen peroxide. Floating seeds were discarded, and the remaining seeds were sown in sterile plastic containers (length 50 cm × width 17 cm × height 15 cm) on sterilized peat treated with chalk and soaked with sterile nutritive solution. Seeds were germinated and seedlings were grown in the climatic growth chamber LT-36VL (Percival) during a 3-month period. Constant temperature of 24 • C, 100% humidity and a long-day (16 h, 30 µmol m −2 s −1 ; fluorescent light 5500 K photoperiod) was used to promote active growth. Seedlings were watered at a two-day interval and fertilizer containing 44 g/L nitrogen, 11 g/L P 2 O 5 , 33 g/L K 2 O, 2.2-3.3 g/L MgO, was used every second week. Two containers for each population were used, and the containers were randomly relocated within the chamber once per week to reduce undesired external effects within in the chamber. After germination, seedlings were thinned out to maintain comparable plant density (approx. 100 seedlings per container). In 3 months after germination, samples for the active growth experimental group were collected, and the remaining seedlings were exposed to the 2-week cold acclimation treatment in the same growth chamber as follows: for one week, the plants were maintained at 10 • C 100% humidity and 8 h photoperiod, followed by one week of darkness with the same temperature and humidity conditions.

Seedling Morphology Measurements
Before the cold acclimation treatment, morphometric traits of 40 seedlings per population were measured to the nearest millimeter. To assess growth cessation during the acclimation treatment, the needle weight of pooled sample representing four repeats for the active growth and acclimated experimental groups of each population was measured. The pooled samples were prepared by combining equal amount of material from first crowns of needles from 25 seedlings.

Two-Dimensional Gel Electrophoresis and Protein Identification
Four biological repeats representing the active growth and acclimated experimental groups of each population were prepared as pooled needle samples as described in the previous section, and were frozen in liquid nitrogen and stored at −70 • C. Protein samples were prepared using phenol extraction and ammonium acetate precipitation, as described previously [48]. Briefly, samples were ground in liquid nitrogen and suspended in 500 µL of the extraction buffer (0.5 M Tris-HCl, pH −7.5, 0.1 M KCl, 0.05 M ethylenediaminetetraacetic acid, 0.7 M sucrose, 2% polyvinylpolypyrrolidone, 2% β-mercaptoethanol, 1 mM phenylmethanesulfonylfluoride). Equal volume of Tris-buffered phenol (pH −7.5) was added and incubated with shaking for 30 min at 4 • C. The tubes were centrifugated for 15 min at 15,000× g at 4 • C, the upper phase was recovered, and the phenol extraction was repeated. Proteins were precipitated by addition of equal volume of 3M ammonium acetate in methanol, incubation at −20 • C overnight and centrifugation as before. Protein pellets were washed twice by resuspending in cold methanol and once in cold acetone and centrifugation for 5 min as before. After aspiration of acetone, protein pellets were dried for few minutes and solubilised in two-dimensional electrophoresis lysis buffer (6 M urea, 2 M thiourea, 4% 3-cholamidopropyl dimethylammonio 1-propanesulfonate, 10 mM Tris-HCl, pH−8.5). Protein concentrations were measured using a Bradford assay [49]. Four biological repeats were prepared for each treatment. Internal standards were prepared from a pooled mixture of all protein extracts.
Protein separation and detection was performed using a differential gel electrophoresis procedure as described previously [50]. Briefly, sample aliquots of 50 mg were used for labeling with Cy3 and Cy5 fluorescent dyes, and the internal standard was labeled with Cy2 dye (Lumiprobe). For the preparative gel, 500 mg of unlabeled internal standard was mixed with 50 mg of the labeled internal standard. Proteins were applied to 24 cm immobilized pH gradient strips with a linear gradient of pH 4-7 and isoelectric focusing was performed with an Ettan IPGphor (GE Healthcare, Chicago, IL, USA). Further, the proteins were separated on 1-mm thick 12.5% polyacrylamide gels with Ettan DALTsix electrophoresis apparatus (GE Healthcare) and gels were scanned using a FLA 9000 fluorescence scanner (GE Healthcare). Relative protein quantifications were performed using DeCyder 2-D Differential Analysis Software, v.7.0 (GE Healthcare).
Preparative gel was fixed in 50% methanol and 10% acetic acid. Protein spots were excised manually and subjected to protein digestion with trypsin, according to a method described previously [51]. Protein digests were loaded and desalted on a 100 mm × 20 mm Acclaim PepMap C18 trap column and separated on a 75 mm × 150 mm Acclaim PepMap C18 column using an Ultimate 3000 rapid separation liquid chromatography system (Thermo-Scientific, Waltham, MA, USA), coupled to a Maxis G4 Q-TOF mass spectrometer detector with a Captive Spray nano-electrospray ionization source (Bruker Daltonics, Bremen, Germany). Peptide identification was performed using the MASCOT server (Matrix Science, Boston, MA, USA) against the Pinus taeda L. genome database [52]. The mass error tolerance for peptide matching was limited to 5 ppm. The threshold value for the identification of proteins was a Mascot score >50 and at least two peptides. Where more than one protein was identified from the same spot, the exponentially modified protein abundance index (emPAI) [53] was employed for the relative quantitation of the proteins. Spots were designated as protein mixture if the emPAI indexes were similar for several of the identified proteins.

Data Analysis
A t-test was used to assess to assess the differences in the seedling morphology traits among the populations (40 seedlings per population).
The Biological Variation Analysis module of the DeCyder software was used to match protein spots of the four biological repeats across different gels and to perform statistical analysis of differential protein expression. A two-way analysis of variance (ANOVA) test was used to assess statistically significant (p ≤ 0.01) population and cold acclimation treatment effects. The threshold value of at least 2-fold difference in protein abundance was used.
To find patterns in the expression profiles, a cluster analysis was performed using Euclidian distance matrix and complete linkage method at the Extended Data Analysis module of the DeCyder software.
Blast2GO software [54] was used for the annotation and gene ontology analysis of the identified protein sequences. The identified Pinus taeda sequences were searched with Basic Local Alignment Search Tool (BLAST) against Swissprot database and protein sequences and annotated with Gene Ontology (GO) terms.

Seedling Morphology
After 3-month cultivation under growth-promoting conditions (16 h light, high humidity and constant temperature of 24 • C), the seedlings of the Finnish population had the shortest length of stem and cotyledons (Table 1). There was no significant difference in the seedling height between the Lithuanian and Spanish populations. Seedlings of the Spanish population had more cotyledons and produced the largest needle biomass, as estimated by fresh weight of the sampled needles. During the acclimation treatment, there was no significant change of needle fresh weight for the seedlings from the Spanish population, indicating growth cessation. Meanwhile, the needle mass increased approximately 1,5-fold for the two northern populations. Table 1. Quantitative traits of pine seedlings before and after the cold acclimation treatment.

Differentially Expressed Proteins
The average number of detected protein spots was 1464 ± 145 per gel after alignment ( Figure S1). Analysis of biological variance revealed statistically significant and >2-fold variation in the abundance of 76 proteoforms among the six experimental groups (non-aclimated and acclimated seedlings from the Spanish, Lithuanian and Finnish populations) (Figures 2 and 3). The largest number of unique differentially expressed protein spots were a characteristic of the Finnish population (25 and 20 in the non-acclimated and acclimated experimental groups, respectively), followed by the Spanish population (15 and 19). The Lithuanian population had the lowest number of unique spots (2 and 1) (Figure 2a,b). The number of proteoforms differentially expressed between populations was similar for the non-acclimated and acclimated experimental group. The abundance of 53 proteoforms changed after cold acclimation treatment (Figure 2c), and approx. 72% of the detected changes were independent of the seedling origin.

Differentially Expressed Proteins
The average number of detected protein spots was 1464 ± 145 per gel after alignment ( Figure  S1). Analysis of biological variance revealed statistically significant and >2-fold variation in the abundance of 76 proteoforms among the six experimental groups (non-aclimated and acclimated seedlings from the Spanish, Lithuanian and Finnish populations) (Figures 2 and 3). The largest number of unique differentially expressed protein spots were a characteristic of the Finnish population (25 and 20 in the non-acclimated and acclimated experimental groups, respectively), followed by the Spanish population (15 and 19). The Lithuanian population had the lowest number of unique spots (2 and 1) (Figure 2a,b). The number of proteoforms differentially expressed between populations was similar for the non-acclimated and acclimated experimental group. The abundance of 53 proteoforms changed after cold acclimation treatment (Figure 2c), and approx. 72% of the detected changes were independent of the seedling origin.  Quantitative distribution of differentially expressed protein spots identified by two-dimensional gel electrophoresis analysis of needle protein samples isolated from seedlings originating from the Spanish (SP), Finish (FI) and Lithuanian (LT) Scots pine populations. Number of differentially expressed protein spots between the populations under growth promoting conditions (a) and upon acclimation (b) is shown. Number of upregulated and downregulated protein spots upon acclimation treatment are indicated by numbers separated by slash (c). The portion of the cycles overlapping between two or all three populations indicate a number of protein spots that are differentially expressed between the two or among all three populations, respectively.

Figure 3.
Quantitative differences in protein abundance between non-aclimated and acclimated Scots pine seedlings originating from the geographically separated populations. Representative images of ten identified spots are shown. Complete gel image is presented in Figure S1.

Protein Expression Patterns
Cluster analysis of variation of abundance of the 76 differentially expressed protein spots revealed protein expression patterns associated with pine seedling origin and/or their response to acclimation treatment ( Figure 4). Cold acclimation treatment had little or no significant effect on the expression of the proteins assigned to the clusters 4-5, 7, 9 and 11-12 and these mostly reflected population-depended protein expression differences. The most prominent difference among the three populations was observed for the two proteoforms assigned to the cluster 12. The lactoylglutathione lyase had approximately 5-to 10-fold higher abundance in the seedlings from the Spanish population as compared to the other two populations. Four proteoforms included in the cluster 11 had higher abundance in the Spanish population (none of them could be identified). Quantitative differences in protein abundance between non-aclimated and acclimated Scots pine seedlings originating from the geographically separated populations. Representative images of ten identified spots are shown. Complete gel image is presented in Figure S1.

Protein Expression Patterns
Cluster analysis of variation of abundance of the 76 differentially expressed protein spots revealed protein expression patterns associated with pine seedling origin and/or their response to acclimation treatment ( Figure 4). Cold acclimation treatment had little or no significant effect on the expression of the proteins assigned to the clusters 4-5, 7, 9 and 11-12 and these mostly reflected population-depended protein expression differences. The most prominent difference among the three populations was observed for the two proteoforms assigned to the cluster 12. The lactoylglutathione lyase had approximately 5-to 10-fold higher abundance in the seedlings from the Spanish population as compared to the other two populations. Four proteoforms included in the cluster 11 had higher abundance in the Spanish population (none of them could be identified). Furthermore, six proteoforms (including tubulin beta chain-like and tau class glutathione s-transferase (GST) were assigned to the cluster 9 and had higher abundance in the Spanish and Lithuanian populations.
Forests 2020, 11, x FOR PEER REVIEW 8 of 21 Furthermore, six proteoforms (including tubulin beta chain-like and tau class glutathione s-transferase (GST) were assigned to the cluster 9 and had higher abundance in the Spanish and Lithuanian populations.  Table S1). GO terms associated with the identified unique protein sequences are shown on the right (the terms corresponding to the indicated numbers are listed in the Results section).  Table S1). GO terms associated with the identified unique protein sequences are shown on the right (the terms corresponding to the indicated numbers are listed in the Results section).
Expression pattern of the remaining clusters 1-3, 6, 8, 10 and 13 reflected population-and/or acclimation-dependent differences (Figure 4). Before the acclimation, the majority of the 18 proteoforms separated into the clusters 1, 3 and 13 (including CP29-like, thaumatine-like protein (TLP), cyanate hydratase and two proteoforms of pathogenesis-related 10 (PR-10) were significantly more abundant in the seedlings of the Finnish or both northern populations. These proteins were significantly upregulated in the acclimated seedlings from all the three populations, however, for cluster 3, the response was more prominent in the Finnish and/or in both northern populations. Cyanate hydratase (cluster 13) was upregulated approximately 10-fold in the acclimated seedlings of the Finnish population compared to the non-acclimated control, representing the strongest response to the acclimation treatment.
Among the 10 proteoforms included in cluster 2, only two significant differences (unidentified proteoforms 7 and 8) were detected among the populations before the acclimation treatment. All proteins of this cluster were upregulated by the acclimation treatment, and the response was more prominent in the seedlings for the northern populations, leading to increase of significant differences among the populations (protein disulfide-isomerase (PDI), nucleoside diphosphate kinase 1 (NDK1) and three unidentified proteoforms).
Furthermore, an abundance of 14 proteoforms included in the clusters 6, 8 and 10 was markedly decreased upon the acclimation treatment. The response to acclimation had similar effect on protein expression for all three populations, however, the protein abundance pattern was different. The northern populations had higher abundance of the guanosine-diphosphate (GDP)-mannose epimerase 2 and proteoform 49 among the three proteoforms of the cluster 6. Meanwhile the remaining proteins (clusters 8 and 10) including probable linoleate 9s-LOX 5, aconitase, TIC 62, oxygen-evolving enhancer protein (OEE), nicotinamide adenine dinucleotide (NAD)-binding RF superfamily protein, abscisic acid water deficit stress and ripening-like (ASR-like) and water deficit inducible lipoprotein 3 (LP3)-like, were more abundant in the samples from the Spanish population.
Notably, only a few proteoforms (including four identified proteins) showed a strictly population-specific regulation upon the acclimation treatment. Malate dehydrogenase was upregulated in both northern populations, although the acclimation-induced expression level in the Lithuanian population was still below the constitutive expression level in the Spanish population. Calreticulin-3-like isoform and alcohol dehydrogenase were upregulated in the Finnish and Lithuanian populations, respectively. Meanwhile, the catalytic subunit of Ser/Thr-PP2A was downregulated in the northern populations, and the protein abundance was reduced to the level that was comparable to the constitutive expression observed in the Spanish population.

Discussion
An adaptive capacity of species in response to specific environmental cues could be inferred from mechanistic information at molecular level and its integration with ecologically relevant phenotypes [55]. A previous study with conifer tree species by [19] has demonstrated that proteomic analysis of young seedlings imparts a variation of gene expression that could be of importance for adaptation to contrasting altitude conditions in the closely related Norway spruce ecotypes, and the adaptive capacity could be further explored by heat stress treatment of the seedlings [18]. Our study revealed a similar quantitative variation of protein expression profiles in the needles of young Scots pine seedlings from geographically distant populations, although it reflects a variation of different physiological traits of potentially adaptive significance that could be species specific. Conifer seedlings develop an acclimation capacity in response to a short photoperiod and low temperatures at an early stage of of development [56], therefore cold acclimation treatment of the seedlings have stimulated additional differences in protein expression profiles that reflect a variation in cold adaptation capacity of the pine genotypes.

Stress Response Proteins Upregulated in the Northern Populations under Active Growth Conditions
Under active growth promoting conditions, both northern populations and in particular the Finnish population had higher abundance of proteins involved in protein metabolism (synthesis or degradation), such as heat shock protein (HSP), calreticulin, structural constituents of 40S ribosomes and 26S proteasome (Figure 4), that are also known as constituents of the stress tolerance mechanism. The chaperone protein 1 is a 100 kDa HSP that is critical for thermotolerance [57] and plant development [58]. Studies with Arabidopsis mutants revealed that constitutive expression of Hsp100 leads to preformed adaptation that ensures effective protection under heat stress [59]. Calreticulin is another chaperone protein which folds newly synthesized proteins and is also involved in regulation of Ca 2+ homeostasis in the endoplasmic reticulum [60]. Upregulation of calreticulin has been observed under salt and osmotic stress conditions in Arabidopsis [61], and a function of the signaling component involved in response to cold stress in rice was proposed for the protein [62]. The proteasome subunit beta type-5 is a component of the 26S proteasome that mediates ubiquitin-dependent proteolysis of proteins and plays an important role in the stress-induced degradation of misfolded proteins, as well as the regulation of stress response signaling [63]. The 40s ribosomal protein sa-like is a structural constituent of ribosomes that regulates protein synthesis. A function of the p40 protein homologue is crucial for the active tissue growth [64] and its transcriptional upregulation under salt stress has been described in Arabidopsis [65].
Furthermore, upregulation of the TLP and PR-10, has been unique to the Finnish population. TLP proteins belong to the PR5 protein family and are induced under biotic, osmotic or cold stress in a variety of plants [66][67][68]. It is also well established that the PR-10 proteins are involved in biotic and abiotic stress response [69]. In perennial plants, expression of the PR-10 is induced during winter dormancy in mulbery and it provides either cryoprotectection for freeze-labile enzymes or serves as a nitrogen-storage protein [70].
In addition, the seedlings from the Finnish population have had higher abundance of alcohol dehydrogenase, malate dehydrogenase, GDP-mannose-epimerase 2, and cyanate hydratase enzymes that are involved in the primary cell metabolism, but also play role in the stress tolerance. Upregulation of alcohol dehydrogenase is essential in anaerobic respiration under cold stress conditions [71,72]. Malate is important intermediate and regulator of metabolic processes. Plants have several malate dehydrogenases that catalyze conversion of malate and oxaloacetate coupled to reduction or oxidation of the NAD pool, and this function has been linked to an abiotic stress tolerance in alfalfa and apple [73,74]. Upregulation of the cyanate hydratase upon salt stress was observed in tomato [75] and Suaeda aegyptiaca (Hasselq.) Zohary plants [76]. An accumulation of cyanide compounds is associated with elevated levels of ethylene synthesis as a consequence of plant stress response, therefore cyanate hydratase is considered a key enzyme responsible for the cyanide detoxification [77,78]. GDP-mannose-epimerase converts GDP-D-mannose to GDP-L-galactose, a precursor of ascorbate and cell wall polysaccharides [79]. This branching metabolic pathway plays a dual role in neutralization of the reactive oxygen species [80] as well as plant development [81].
The seedlings of the northern populations showed shorter phenotype and accumulated less needle biomass as compared to the Spanish population under active growth conditions (Table 1). It could be presumed that the reduced growth and upregulation of the stress response proteins might be an indicator of onset of the cold acclimation processes triggered by photoperiod used for the active growth experiments (16 h light and 8 h dark), owing to the short night length critical for growth cessation for the northern populations [34]. However, the seedlings had a relatively high level of the catalytic subunit of PP2A enzyme under active growth conditions that was effectively downregulated by the acclimation treatment. The PP2A is involved in cold stress signaling as a negative modulator of adaptive responses [82], and transcript mRNA levels of the PP2A has been reduced in tomato and alfalfa plants upon cold stress treatment [83,84]. This would suggest that the seedlings of the northern populations maintained active growth under the conditions used in the experiment, and the observed growth and protein expression differences were rather a consequence of physiological differences among the genotypes. Higher abundance of the proteins involved in stress tolerance could be linked to a constitutive cold adaptation mechanism characteristic of the northern genotypes.

Stress Response Proteins Upregulated in the Lithuanian and/or Spanish Population under Active Growth Conditions
The clusters 9 through 12 shown in Figure 4 include 20 proteoforms that have been upregulated mostly in the Lithuanian and/or Spanish populations under active growth conditions. Lactoylglutathione lyase, also known as glyoxalase I, is highly expressed in the Spanish population (approx. 10 and 5-fold compared to the Lithuanian and Finnish populations, respectively). The enzyme detoxifies methylglyoxal, a cytotoxic compound that is produced through glycolysis intermediates and generates free radicals [85]. Levels of methylglyoxal have been shown to rise under salinity stress and are regulated by glutathione and glyoxalase I [86].
As it could be expected due to geographical origin, the Lithuanian population tends to follow a similar protein expression pattern as the Finnish population located further north. However, expression of the proteins assigned to the clusters 9 and 10 (including tubulin beta-chain like, tau class GST, chloroplastic OEE, ASR-and LP3-like) in the Lithuanian population is similar to the Spanish population. The tau class GST, is involved in the reactive oxygen species (ROS) detoxification reactions that conjugate reduced gluthatione with a variety of target compounds. Expression of GSTs in plants is responsive to both biotic and abiotic stresses [87] and overexpression of the enzyme leads to chilling and osmotic stress tolerance in tobacco seedlings [88]. Beta-tubulin is one of the two structural components of cell microtubules and it is intimately involved in regulation of cell morphogenesis. In addition, it has been shown that the ROS signaling-mediated rearrangement in tubulin cytoskeleton plays crucial part in cell adaptation to stress [89]. Further, the abscisic acid induced ASR-like protein plays an important role in drought and salinity stress [90], and overexpression of ASR genes results in an increased cold tolerance [91,92]. LP3-like protein is a water-deficit-induced protein, which is homologous to the ASR and accumulates in loblolly pine roots upon water-deficit stress [93]. OEE is involved in the assembly of the Photosystem II and has a protective effect on the manganese cluster [94]. Increased expression of OEE1 was shown under salt or water deficiency stress [95,96].
The observed variation of stress-related protein expression among the populations from the southern and northern boundaries of the continent indicates presence of the constitutive stress tolerance mechanisms that are pertinent to the climatic conditions of the geographically distant habitats. The Lithuanian population combines a mixture of the traits.

Effect of Acclimation Treatment on Protein Expression
Cold acclimation treatment revealed additional protein expression differences that could be separated into three major groups that include proteins either (i) upregulated or (ii) downregulated in seedlings of all the three populations, or (iii) regulated in a population-specific manner. These protein expression differences suggest the presence of distinct elements of the cold acclimation mechanism and/or variation of their regulation among the geographically distant pine populations.
The function of the first group of proteins was linked to stress response. HSC70-2 is a protein folding catalyst involved in protein stabilization under abiotic stress conditions [97]. PDI belongs to another well-known family of chaperones that mediates formation, isomerization and reduction/oxidation of disulfide bonds during protein synthesis [98]. Ribonucleoprotein chloroplastic-like is required for normal chloroplast development under cold stress conditions by stabilizing transcripts of mRNAs [99]. Nucleoside diphosphate kinase 1 is a cytosolic enzyme that is essential in homeostasis of cellular nucleoside triphosphate pools [100], and its role under cold stress conditions has been described [101]. The stress response related proteins PR-10, TLP and cyanate hydratase have been also highly upregulated upon acclimation treatment, however, their abundance in the Spanish and Lithuanian populations have not exceed the level observed in the Finnish population under active growth conditions. Acclimation treatment downregulates chloroplastic Tic62, a component of Tic translocon that mediated translocation of nuclear-encoded precursor proteins across the inner envelope of chloroplasts, and chloroplastic OEE were downregulated upon the acclimation treatment. Further downregulation of several enzymes, including cytoplasmic aconitase involved in the Krebs cycle, probable linoleate 9S-lipoxygenase 5 of the oxylipin synthesis pathway and GDP-mannose-epimerase 2 enzyme involved in the ascorbate and cell wall polysaccharide synthesis might be the result of reduced metabolic and photosynthetic activity in the seedlings maintained under the acclimation conditions of low temperature and reduced light. However, seasonal changes in chloroplast ultrastructure [102] and depression of photosynthesis [103] is also a part of the winter adaptation process in conifers.
Abundance of ASR-like and homologous LP3-like protein have been also downregulated in the seedling needle samples. This may appear to contradict the previous studies demonstrating that the ASP and LP3 play essential role in stress response [90,91,93]. However, the study by Padmanabhan et al. [93] has shown that the transcript of the LP3 protein is preferentially induced in roots with a constitutive basal level of expression observed in stems and needles, which would explain low expression levels in our experimental set up.
Among the proteins regulated in a population-specific manner, PP2A, the negative modulator of adaptive response [82], has been downregulated, and MDH and calreticulin-3-like isoforms have been upregulated in the northern populations upon the acclimation treatment. In the Spanish population, the response to acclimation treatment has a distinctive reduced intensity compared to the northern populations.

Conclusions
Our study demonstrates that comparative proteomic analysis of young Scots pine seedlings is capable of reflecting on adaptive capacity among the geographically distant populations. Under active growth conditions, quantitative protein expression differences were related to biological processes of cell metabolism and stress response and could represent elements of the adaptive mechanism tailored to the specific growth environment. Furthermore, acclimation treatment revealed an additional stratum of protein expression differences that could provide insight into molecular basis of the cold adaptation capacity among the genotypes. As could be expected, the acclimation-induced protein expression patterns confirmed the presence of a similar cold acclimation mechanism among the pine populations, however the analysis also revealed a discrete variation of the adaptive traits. The most prominent differences were substantiated by the constitutive low levels of the PP2A enzyme expression and the overall less-pronounced response to the acclimation treatment of seedlings from the Spanish population. The identified proteins involved in the constitutive or acclimation-induced expression profiles presented interesting candidates for further studies on cold adaptation physiology in pine.

Supplementary Materials:
The following are available online at http://www.mdpi.com/1999-4907/11/1/89/s1: Figure S1: Two-dimensional electrophoresis gel of the needle protein samples isolated from Scots pine seedlings; Table S1: Proteins differentially expressed in seedlings of the different populations of Scots pine under active growth conditions and upon acclimation treatment that were identified as unknown protein sequences of the Pinus taeda genome or contained several protein mixtures.

Conflicts of Interest:
The authors declare no conflict of interest.   a Sequences that were selected as most abundant sequences from protein mixture using emPAI index. Data presented in graphs as log mean and standard error of the mean of standardized abundance, the same letters indicate statistically significant (p < 0.01) differences between the three populations, the same number of stars indicate significant difference between non-acclimated and acclimated experimental groups within each population. Abbreviations: p. no.-peptide number; SC-sequence coverage; Theoretical MW/pI-theoretical molecular weight and pI values; SP, LT, FI-Scots pine populations from Spain, Lithuania and Finland.
a Sequences that were selected as most abundant sequences from protein mixture using emPAI index. Data presented in graphs as log mean and standard error of the mean of standardized abundance, the same letters indicate statistically significant (p < 0.01) differences between the three populations, the same number of stars indicate significant difference between non-acclimated and acclimated experimental groups within each population. Abbreviations: p. no.-peptide number; SC-sequence coverage; Theoretical MW/pI-theoretical molecular weight and pI values; SP, LT, FI-Scots pine populations from Spain, Lithuania and Finland.