Early phenology and growth trait variation in closely related European pine species

Abstract Closely related taxa occupying different environments are valuable systems for studying evolution. In this study, we examined differences in early phenology (bud set, bud burst) and early growth in a common garden trial of closely related pine species: Pinus sylvestris, P. mugo, and P. uncinata. Seeds for the trial were sourced from populations across the ranges of each species in Europe. Over first 4 years of development, clear differences were observed between species, while the most significant intraspecific differentiation was observed among plants from P. sylvestris populations from continental European locations. Trait differences within P. sylvestris were highly correlated with altitude and latitude of the site of origin. Meanwhile, P. mugo populations from the Carpathians had the earliest bud set and bud flush compared to other populations of the species. Overall, populations from the P. mugo complex from heterogeneous mountain environments and P. sylvestris from the Scottish Highlands showed the highest within‐population variation for the focal traits. Although the three species have been shown to be genetically highly similar, this study reveals large differences in key adaptive traits both among and within species.


| INTRODUCTION
Photoperiod and temperature are important environmental drivers of distribution and local adaptation (Gaston, 2009). In forest trees, annual growth cycles have evolved in response to growing season length variation, making trade-offs between survival and growth (Petit & Hampe, 2006). Across a species' range, many phenotypic traits in trees show steep clinal variation that is usually interpreted as strong evidence for a balance between spatially variable selection and migration (Savolainen, Pyhäjärvi, & Knürr, 2007). In highly outcrossing species such as pines, among-population variation in quantitative traits is often accompanied by very low population structure at neutral genetic loci due to the homogenizing effect of genetic dispersal by pollen and seed (Aitken, Yeaman, Holliday, Wang, & Curtis-McLane, 2008). This indicates that divergent selection is strong enough to maintain population differentiation even in the presence of substantial gene flow (Storz, 2002). Patterns of population differentiation, estimated from quantitative trait variation in samples from across environmental gradients, have been studied in many forest tree species including Quercus petraea (Ducousso, Guyon, & Kremer, 1996), Populus tremula (Luquez et al., 2007), Betula pendula (Viherä-Aarnio, Häkkinen, Partanen, Luomajoki, & Koski, 2005), and Pinus sylvestris (Hurme, Repo, Savolainen, & Paakkonen, 1997;Salmela, Cavers, Cottrell, Iason, & Ennos, 2011). However, comparative common garden studies of multiple closely related tree species have rarely been undertaken and have never previously been conducted simultaneously for the species investigated here. Studying closely related species in this way has the potential to improve our understanding of how the evolutionary process operates on adaptive trait variation across environmental gradients. Furthermore, such phenotypic data are essential to support genomic studies focussed on the genetic architecture of quantitative traits. Considering environmental change predictions, such studies are also important for assessment of mitigation strategies such as assisted migration or marker-assisted breeding, which have been proposed for sustainable management of forest ecosystems.
The three pine species investigated here, Scots pine (Pinus sylvestris L.), dwarf mountain pine (Pinus mugo Turra s.l.), and mountain pine (Pinus uncinata Ramond), differ significantly from each other in phenotype, geographic range, and ecology. The species are phylogenetically closely related, and they are not reproductively isolated and diverged within the last 5 Myr (Wachowiak, Palme, & Savolainen, 2011).
Pinus sylvestris is a tree of up to 40 m height and is socially, ecologically, and economically one of the most important forest tree species in the world. It has the largest range of all pines with populations from western Scotland to eastern Siberia and southern Spain to northern Scandinavia. Pinus mugo grows to a height of about 3 m, and it forms large, shrubby forests above the tree line up to an altitude of about 2,700 m with major populations in mountainous regions of western Europe (the Alps), central and eastern Europe (Sudetes, Carpathians), and southeast through Bosnia and Herzegovina, Montenegro, Serbia, and Romania to the Rila and Pirin mountains of Bulgaria, as well as several outlier populations (Critchfield & Little, 1966). Pinus uncinata is a tree of up to 20-25 m in height found in mountainous regions at altitudes of 1,400-2,700 m in the Pyrenees, the Massif Central, Western Alps as well as several marginal locations in the Iberian Peninsula (Hamerník & Musil, 2007).
Morphological and ecological divergence among these taxa is probably due to adaptation to their different habitats, enforced by the disjunction of their ranges and isolation in different Pleistocene refugia during glacial maxima (Christensen, 1987a,b). The species can be artificially crossed, and they hybridize in contact zones (Jasińska et al., 2010;Wachowiak & Prus-Głowacki, 2008). Given that they currently occupy largely disjunct geographic ranges in ecologically distinct habitats, local adaptation has likely played a major role in their phenotypic diversification.
Although these species are clearly differently adapted and can be recognized based on growth form and size, previous studies have shown that they share a highly similar genetic background.
Given this genetic background, we expect that the genetic architecture, and the set of underlying genes, governing adaptive trait variation within and among species should be similar when considered across similar environmental gradients. Despite their disjunct ranges, the species overlap to some extent in altitude, temperature, and rainfall range. Here, we might expect patterns of adaptive phenotypic variation to be similar in each species. Alternatively, adaptation to specific environmental variables may have involved selection for species-specific or high-frequency divergent alleles. In this case, levels of intraspecific variation for the traits examined would be expected to differ between taxa.
In this study, we looked at patterns of variation in adaptively important quantitative traits across the geographic distribution of the species in Europe. Phenology (bud set, bud burst) and growth (height) data were collected over a 4-year period from a progeny trial of young pines representing 28 populations across three species. We analyzed within-and among-species variation to assess the strength of local adaptation at different ecological and evolutionary scales. We examined the extent of phenotypic differentiation in the traits at intra-and interspecific level. Specifically, we asked what is the amount of phenotypic variation within species considering their similar genetic background but contrasting distribution ranges.

| Sampling
Open-pollinated seeds of the three pine species were collected from twenty-eight natural populations in Europe covering the extent of each species range, including thirteen populations of P. sylvestris, nine P. mugo, and six P. uncinata (Table 1, Figure 1). Cones were collected from five mature trees from each population. Trees were separated by at least 50 m to minimize sampling of closely related individuals.
Furthermore, the trees were selected to represent typical stands from each location avoiding sampling from population margins or extreme habitats within each stand especially in the mountain areas. Species and populations were sampled across environmental gradients in latitude (temperature, photoperiod), altitude (mountain, lowland), and growing season length to maximize the environmental amplitude that has given rise to adaptive responses in phenology and growth. Trees need to balance growth and dormancy periods to successfully compete for resources within populations while avoiding frost damage.
We focused on these traits as they are under strong selective pressure in natural populations.

| Progeny trial experiments
Genetic variation in phenotypic traits can be identified by growing young pines in a common garden environment in which environmental variation is kept to minimum (e.g., White, Adams, & Neale, 2007).
Seeds from common mother trees were sown on trays (75:25 compost type John Innes #1: sand) in spring 2010. After germination, a provenance-progeny trial was established in an unheated glasshouse at the Centre for Ecology and Hydrology, Edinburgh, UK. Seedlings were transferred to 0.62 l pots (diameter 11 cm, depth 9.6 cm) and grown under natural light conditions (glasshouse was shaded to avoid excess light) with automatic watering applied during the growing season. The glasshouse temperature was set to track outside temperature ±3°C, and watering was automatic at 4 times a day for 3 min 30 s-typically, this was sufficient to keep pots moist; however, a control pot containing a soil humidity probe was used to regulate a minimum soil moisture (70%), below which additional watering would be applied.
For most populations, the trial consisted of five families/population and the trial was divided into 18 randomized blocks (Table 1), with a single seedling per family per block. In total, 1623 young pines were analyzed including 392 from P. sylvestris (continental), 365 from P. sylvestris (Scotland), 399 from P. mugo, and 467 from P. uncinata. The average 2010-2011, 2011-2012, and 2012-2013  T A B L E 1 Geographic location, altitude, and mean annual temperature for the sampled pine species populations from Europe

| Quantitative trait data analysis
Phenology (timing of bud set/ bud burst) and total height were recorded for every seedling to evaluate within-and between-species variation. Bud set was measured from August to October 2010 as the number of days since the first scoring date (i.e., since the date on which the first plant that set terminal bud was observed). Bud set was scored when a visible apical bud with clearly developed scales was formed at the tip of a stem in each seedling. Bud burst was scored as emergence of the new needles around the tip of the apical bud in the main stem. It was measured during spring 2011 and 2012 starting from the day when the first bud burst was observed.
Phenology observations were conducted twice a week. The height of young pines was measured in December 2011, 2012, and 2013.
As the numbers of populations and families within species varied, we employed statistical models that were robust to unbalanced designs. Unbalanced analysis of variance (families nested within populations) was used to test for between species, population, and family differences in phenology and growth performance. A group of six Scottish populations of P. sylvestris were analyzed separately as they cross a strong west-east environmental gradient at a narrow spatial scale (Salmela et al., 2010). This approach allowed us to evaluate adaptation at a local scale without introducing bias in the analysis of populations across the full range due to overrepresentation of Scottish populations. Populations were considered fixed factors, and families within populations and blocks were random.
Variance components due to populations, families, and blocks were estimated using a mixed-model (REML) approach. Similarly, interspecies comparisons were conducted treating species as a fixed factor. Analyses were carried out using GenStat 13.1.0.4470.

| Environmental covariance
We investigated whether differences in climatic conditions have induced phenotypic adaptive variation among natural populations by looking at trait variation associated with altitude and latitude. As variation in the traits analyzed is related strongly to location, and because geographic descriptors capture unmeasured environmental factors that might also be selectively important, we used latitude and altitude as predictors for the investigated traits. Data on mean temperature, which affects species phenology and growth, were also extracted from Met Office and WorldClim databases (Table 1), although, in general, values of mean temperature were correlated with the geographic descriptors of each site (Table S1).

F I G U R E 1 Geographic location of the sampled pine populations across Europe
Mixed-model analyses were performed in R (www.R-project.org/), using a methodology similar to that described by Grueber, Nakagawa, Laws, and Jamieson (2011). For each trait, a global linear mixed-effects model was first specified using the package "lme4" (Bates, Mächler, Bolker, & Walker, 2015) as follows: where the covariates latitude and altitude (m) at sites of origin were fixed effects; population, family within population, and block were treated as random effects with random intercepts. After model fitting, covariates were rescaled to a mean of zero and a standard deviation of 0.5 using the package "arm" (Gelman & Su, 2014). Full model subsets were generated and evaluated by AICc using the package "MuMIn" (Barton, 2015) (R script is provided in Appendix S1). In addition, for each species/trait, a random-effects model was specified to evaluate the contribution of population, family within population, and block to the total variance in the absence of additional covariates. Covariation between selected pairs of traits was examined in a similar manner as described for environmental covariates. The first of each trait pair was defined as the dependent variable; models with and without the second trait as a covariate were then compared in terms of AICc.
Covariates were retained on their original scales. Population, family within population, and block were fitted as random effects.

| Phenotypic differentiation between species
There was a significant difference (p < .0016 for all traits) between species in phenology and growth ( Figure 2, Table S2). On average, the mountain species, P. mugo, had the earliest bud set; however, this still occurred later than the most northern populations of P. sylvestris ( Figure 3, Table S3). Populations of P. sylvestris from northern Sweden set terminal buds earlier than any other population across all species (Figure 3). Most of the P. sylvestris populations set terminal buds later than other species. Pinus mugo was the first to terminate dormancy and had the slowest growth ( Figure 2

| Phenotypic differentiation within species
There was significant between-population variation in phenology among P. sylvestris from mainland locations and among P. mugo populations, but not among Scottish P. sylvestris populations or among P. uncinata populations (Table S4). Continental P. sylvestris and P. mugo had large ranges of variation among population means for bud flush (from about two to 4 weeks between first and last scores) ( Figure 3, Table S3).  (Table S4). In mainland P. sylvestris populations, differentiation between populations was the major component of variation for all traits throughout the experiment ( Figure 4, Table S5).
F I G U R E 3 Bud set and bud burst (a) and growth (height) (b) variation between populations and their standard errors (Table 1) In  Table S5).
Pinus uncinata was the most phenotypically homogenous and did not show any significant variation between populations for phenological traits (Table S4). Significant between-population differences in height in this species were driven by one population from Spain (Castiello de Jaca) that showed the slowest growth compared to other populations.
Large differences between families within populations were observed for bud set in P. sylvestris and P. uncinata and for bud burst in P. uncinata in the second year of measurements (Table S4). There was significant difference in height between families within populations in each species in the second year of measurements (Table S4).

| Environmental associations and covariation between traits
Overall, the populations of origin of our P. sylvestris plants were distributed across a much larger latitudinal range and at lower elevations than the taxa from the P. mugo complex ( Figure S1), and altitude and latitude were highly negatively correlated (r 2 = .82, p < .01) (Table S1).
However, for continental populations of P. sylvestris, evidence was found for an effect of altitude and latitude of population origin on both bud set and height (Table 2) (Table S6).

| Phenotypic variation at interspecific level
The focal pine species diverged from a common ancestor about five million years ago (Wachowiak, Palme, et al., 2011). At present, they have mostly disjunct ranges and are adapted to different environmental F I G U R E 4 Variance components for phenology and growth in the three pine species gradients from northern to southern Europe that experience differing selective pressures in relation to temperature, photoperiod, and water availability. In our study, we grew individuals from seed sampled from populations located across the European range. Trial plants were raised under controlled environmental conditions in a glasshouse to assess variation in key phenotypic traits relevant to species evolution and intraspecific local adaptation. We focused on using geographic descriptors (latitude, altitude) to assess trait variation as these have strong relationships with variation in the phenology and growth traits studied. Common garden experiments are useful for studying genetic variation in phenotypic traits as they minimize variation in environmental conditions. In our case, glasshouse conditions did not reflect those at the home sites of any of the provenances studied and were effectively novel to all test plants; therefore, we expected expression of latent variation and a wider range of variability in phenotype that might be expected at their home sites. However, as conditions were common to all plants, any significant inter-or intraspecific differences observed in our study in the analyzed traits should have a strong heritable component. Although not specifically addressed by our experimental design, it cannot be excluded that some of the observed phenotypic differentiation between populations and species could be affected by epigenetic interactions. Such heritable changes in gene function that do not result explicitly from sequence variation may contribute to phenotypic plasticity and adaptation. However, the influence of epigenetics effects on phenotypic diversity is much less known in pines that may not be affected by epigenetic interaction in a similar fashion like other conifers, e.g., spruce species (Gömöry, Foffová, Longauer, & Krajmerová, 2015;Kohmann & Johnsen, 1994;Yakovlev, Fossdal, & Johnsen, 2010). Also, side effects of the maternal environment are at least partially accounted for by the family component in our models .
Our data showed that these closely related pine species have evolved significant differences in phenology and height. In general, P. sylvestris was phenotypically distinct from the other two species.
However, there was some overlap in phenology between species.
For instance, populations of P. sylvestris from Finland had similar bud set timing to P. mugo populations from the Carpathians, and P. mugo from Austria had similar bud set timing to P. uncinata populations.
We also observed some similar patterns of covariation between traits T A B L E 2 Effects of environmental variables (altitude and latitude; please note that the coefficients are standardized) on phenotypic variation based on mixed-model analysis (REML). ΔAICc values represent the difference in AICc between the null model (no fixed effects) and the best (a value of zero indicates that the null model was the best); Akaike weight represents the probability of that model being the best out of the set considered; df, degrees of freedom in different species. Differences in height among populations were more consistent than differences in phenological traits. Pinus sylvestris showed more rapid growth than the species from the P. mugo complex. Reduced height can be advantageous in mountain regions due to the risk of damage from strong wind and avalanches. This is particularly true for P. mugo, which forms extensive shrubby populations at high altitudes above the tree line, and had the slowest growth compared to other species. The lower level of variation in phenology and growth between taxa from the P. mugo complex than between the latter and P. sylvestris corresponds to adaptation to similar and relatively narrow environmental niche occupied by the species from the P. mugo complex. These patterns are reflected in estimates of evolutionary divergence from biometric, biochemical, and molecular data (Boratynska, Jasinska, & Boratynski, 2015;Lewandowski, Boratyński, & Mejnartowicz, 2000). Despite differences in growth form, the taxa from the P. mugo complex have shown no diagnostic anatomical or morphological needle traits, or molecular markers (Boratynska et al., 2015;Jasińska et al., 2010;Wachowiak, Odrzykoski, Myczko, & Prus-Głowacki, 2006). Therefore, some biometric markers often used in pine taxonomy could not be applied to verify, e.g., hybrid origin of trees from contact zones of the species (Boratynska et al., 2015).

| Patterns of within-and between-population differentiation
At the within-species level, the taxa showed different patterns of within-and between-population variations. In P. sylvestris, we found substantial between-population variation that was highly correlated with latitude and/or altitude of population origin. In general, northern populations from Finland and Sweden set bud and flushed earlier than more southerly populations following general patterns in forest trees related to photoperiod and temperature differences (Hurme et al., 1997;Rohde, Bastien, & Boerjan, 2011). On the other hand, populations from Spain, Scotland, and northern Europe (Finland and Sweden) showed slower growth compared to central European (Poland and Austria) and Italian populations. In previous studies, significant phenotypic differentiation between distant populations of P. sylvestris was observed in timing of growth cessation and frost hardiness (Koski & Sievänen, 1985). When seedlings from different parts of Europe were grown under photoperiods typical of 50° latitude, significant differences were observed in the timing of bud set of young pines from northern versus southern regions (Oleksyn et al., 1992) and height growth cessation of older trees (Oleksyn, et al. 1998;Repo et al., 2000). Furthermore, provenance transfers of P. sylvestris from north to south of Europe resulted in increased survival but reduced growth as compared to local populations (Eriksson et al., 1980;Persson & Ståhl, 1990), supporting the observations of clinal variation of phenotypic traits in the species (Langlet, 1959). In contrast, no clear correlation with latitudinal clines was observed in our study for the two taxa from the P. mugo complex. Some evidence of between-population differentiation in bud burst (third year of measurements) weakly associated with environmental gradients was found in P. mugo. The most southerly populations, from Italy, set terminal buds about 5 days earlier than the most northerly populations from Poland. The most homogenous species in our dataset was P. uncinata. Based on microsatellite loci (Dzialuk et al., 2009), populations of P. uncinata differ genetically across its range, presumably as a result of different postglacial histories. However, our data do not show any differences in phenology or height between P. uncinata populations. With the exception of one outlier population, which showed slower growth, there was no significant between-population differentiation within P. uncinata.
In general, variation between families within populations was much more pronounced in the taxa from the P. mugo complex and P. sylvestris from Scotland compared with the mainland range of the latter. In Between-population variation within species has been observed in previous morphological and biochemical studies. Differences were detected between east and south Carpathian and other populations of P. mugo in a biometric study of needle and cone characters (Boratynska et al., 2015) and allozyme loci (Sannikov, Petrova, Schweingruber, Egorov, & Parpan, 2011). Differentiation was explained as being due to different Holocene histories and Pleistocene isolation that influenced genetic variation between those parts of the species range. Biometric data also indicated differentiation of P. mugo populations from the Balkans, represented in our dataset by the Pirin Mountain population. This was explained as being due to colonization of the Balkans from two different refugia in the Alps and southern Carpathians (Boratynska et al., 2015). However, those populations showed no evidence of significant genetic divergence at neutral markers when compared to putative source populations (Wachowiak et al., 2013). Some differences between Pyrenean and marginal populations of P. uncinata in the Massif Central and Spain were found in previous biometric and genetic studies (Boratynska et al., 2015;Dzialuk et al., 2009). These studies also found some evidence of morphological differences among

| Underpinning future searches for genes of adaptive significance
Recent studies of nucleotide sequence variation showed high genetic similarity between the focal species: they share a high proportion of common, ancestral variation (Wachowiak et al., 2013;Wachowiak, Palme, et al., 2011). The combination of recent evolutionary divergence and substantial ecological differentiation means these species form a valuable model for studies of the genetic basis of local adaptation and speciation. Our data show that the mean and level of variation in important phenotypic traits differ between taxa. Given this evidence for quantitative phenotypic divergence, studies of variation at the genome level can now be better focussed on searches for genes that have experienced selection during the process of speciation. Correspondingly, linking analysis of withinspecies phenotypic variation to analysis of patterns of molecular polymorphism will help to verify which genes are under selection due to adaptation to local environmental conditions. However, as many phenotypic traits are polygenic, searches for genes involved in local adaptation and speciation must involve many genomic regions. Recently developed genomic resources for P. sylvestris and the taxa from the P. mugo complex, comprising whole transcriptome sequences and a large database (over 200,000) of single nucleotide polymorphisms (SNPs), provide a valuable set of markers for such analysis (Wachowiak, Trivedi, Perry, & Cavers, 2015). Such studies should also address epigenetic variation that may be important in local adaptation for long-lived organisms like forest trees, which experience substantial environmental variation during their life spans.

| CONCLUSION
Our data demonstrate significant between-species variation in phenology and height. Distinct intraspecific patterns of phenotypic variation indicate that selection has operated differently on traits important for adaptation to the environments occupied by individual species. Approaches for resolving the genetic architecture of adaptive traits in species derived from heterogeneous environments such as the taxa from the P. mugo complex should take into account the high within-population variation of the studied traits.
Comparison of genetic variation between populations could be effective in species like P. sylvestris which exhibit strong clines in adaptive traits, assuming that they represent populations of shared evolutionary and demographic history (Wachowiak et al. 2014). The study demonstrates genetically driven variation in pine phenotypes across the species distribution in Europe. This variation is important for population resilience, and better knowledge about its distribution and effects can advance management practice as forests are challenged by future environmental changes.