# Within and between population comparisons suggest independently acting selection maintaining parallel clines in Scots pine (Pinus sylvestris)

Sonja T Kujala, Komlan Avia, Timo A Kumpula, Hannu Kärkkäinen, Juha

Heikkinen, Katri Kärkkäinen, Outi Savolainen

<figure><img src="image_1.png" /></figure>

## ► To cite this version:

Sonja T Kujala, Komlan Avia, Timo A Kumpula, Hannu Kärkkäinen, Juha Heikkinen, et al. Within and between population comparisons suggest independently acting selection maintaining parallel clines in Scots pine (Pinus sylvestris). Evolution Letters, In press, 10.1093/evlett/qrad054. hal-04232435v2

HAL Id: hal-04232435

https://hal.inrae.fr/hal-04232435v2

Submitted on 15 Apr 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
# Within and between population comparisons suggest independently acting selection maintaining parallel clines in Scots pine (Pinus sylvestris)

Sonja T Kujala, Komlan Avia, Timo A Kumpula, Hannu Kärkkäinen, Juha

Heikkinen, Katri Kärkkäinen, Outi Savolainen

► To cite this version:

Sonja T Kujala, Komlan Avia, Timo A Kumpula, Hannu Kärkkäinen, Juha Heikkinen, et al. Within and between population comparisons suggest independently acting selection maintaining parallel clines in Scots pine (Pinus sylvestris). Evolution Letters, 2023, 8 (2), pp.231-242. 10.1093/evlett/qrad054. hal-04232435v2

HAL Id: hal-04232435

https://hal.inrae.fr/hal-04232435v2

Submitted on 15 Apr 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Within- and between-population comparisons suggest independently acting selection maintaining parallel clines in Scots pine (Pinus sylvestris)

Sonja T. Kujala<fnref n="1" />, Komlan Avia<fnref n="2" />,<fnref n="3" />, Timo A. Kumpula<fnref n="2" />, Hanni Kärkkäinen<fnref n="4" />, Juha Heikkinen<fnref n="5" />, Katri Kärkkäinen<fnref n="1" />, Outi Savolainen<fnref n="2" />

Corresponding author: Natural Resources Institute Finland, Paavo Havaksen tie 3, 90570 Oulu, Finland. Email: sonja.kujala@luke.fi

## Abstract

Parallel clines in traits related to adaptation in a species can be due to independent selection on a pair of traits, or due to selection in one trait resulting in a parallel cline in a correlated trait. To distinguish between the mechanisms giving rise to parallel adaptive population divergence of multiple traits along an environmental gradient we need to study variation, correlations, and selective forces within individual populations along the gradient. In many tree species, budset timing (BST) forms a latitudinal cline, and parallel clinal variation is also found in other seedling traits, such as first-year height (FYH) and fall frost injury (FFI). In this study, we set up a common garden experiment with open pollinated progeny from natural populations of Scots pine (Pinus sylvestris), with one large sample from single population (500 families) and smaller samples from across a latitudinal gradient. BST, FYH and induced FFI were first measured in a greenhouse. The seedlings were then planted in the field, where survival and height were measured at the age of 9 years as fitness proxies. We compared between- and within-population variation and genetic correlations of these three seedling traits, and estimated selection gradients at the family level in our main population, taking into account the potential effects of seed weight. Between-population genetic correlations between seedling traits were high (0.76–0.95). Within-population genetic correlations in the main population were lower (0.14–0.35), as in other populations (0.10–0.39). Within population, extensive adaptive variation persists in the seedling traits, in line with rather weak selection gradients, yet maintaining the clines. Although our sampling does not cover the whole cline equally, the results suggest that the individual clines in these traits are maintained by largely independently acting selection, which results in fewer constraints in adaptation under changing climate.

**Keywords:** budset timing, fall frost injury, first-year height, genetic correlation, local adaptation, selection gradient

## Lay summary

Clines in genetically governed traits result from adaptation to gradually varying environmental gradients. Genetic dependencies between traits can influence adaptation. As the climate warms, northern Scots pine must adapt to new combinations of photoperiod and temperature to time its yearly growth period optimally along the cline. Here, we examined the dependency between adaptive traits clines in budset timing, height, and cold tolerance of first-year Scots pine seedlings. While these clines are highly correlated between populations along the cline, within-population analysis with a large natural population sample revealed much weaker genetic correlations. Further, selection was shown to act mostly independently on each trait. The lack of strong genetic correlations between adaptive traits results in fewer constraints on adaptation.

## Introduction

When environmental conditions change gradually over an area, for example, latitudinally, adaptation can be reflected as clinal patterns where population means of adaptation-related traits track local optima (Felsenstein, 1977). Clinal variation in many traits of plants and animals has been observed throughout centuries of research (Endler, 1977). Where overall genomic variation suggests little or no clinal differentiation, the clines likely have been formed by the action of natural selection (Chen et al., 2014; Kujala & Knürr et al., 2017; Mimura & Aitken, 2007; Savolainen

et al., 2007; Viherä-Aarnio et al., 2005). When populations along the cline are connected through current gene flow, clines reflect a balance between selection and gene flow (Slatkin, 1973, 1978).

Many different traits within a species may display similar clines, raising the question whether selection has acted independently on the traits, or selection on one trait has resulted in correlated clines in other traits through physiological constraints and/or pleiotropy and linkage disequilibrium (indirect selection; Guillaume, 2011; Lande, 1979; Merilä & Björklund, 1999). Kingsolver et al. (2001) did not find much evidence of indirect

<footnote n="1">¹Natural Resources Institute Finland, Oulu, Finland</footnote>

<footnote n="2">²University of Oulu, Ecology and Genetics, Oulu, Finland</footnote>

<footnote n="3">³Université de Strasbourg, INRAE, SVQV UMR-A 1131, Colmar, France</footnote>

<footnote n="4">⁴Natural Resources Institute Finland, Jokioinen, Finland</footnote>

<footnote n="5">⁵Natural Resources Institute Finland, Helsinki, Finland</footnote>
# Within- and between-population comparisons suggest independently acting selection maintaining parallel clines in Scots pine (Pinus sylvestris)

Sonja T. Kujala<fnref n="1" />, Komlan Avia<fnref n="2" />,<fnref n="3" />, Timo A. Kumpula<fnref n="2" />, Hanni Kärkkäinen<fnref n="4" />, Juha Heikkinen<fnref n="5" />, Katri Kärkkäinen<fnref n="1" />, Outi Savolainen<fnref n="2" />

Corresponding author: Natural Resources Institute Finland, Paavo Havaksen tie 3, 90570 Oulu, Finland. Email: sonja.kujala@luke.fi

## Abstract

Parallel clines in traits related to adaptation in a species can be due to independent selection on a pair of traits, or due to selection in one trait resulting in a parallel cline in a correlated trait. To distinguish between the mechanisms giving rise to parallel adaptive population divergence of multiple traits along an environmental gradient we need to study variation, correlations, and selective forces within individual populations along the gradient. In many tree species, budset timing (BST) forms a latitudinal cline, and parallel clinal variation is also found in other seedling traits, such as first-year height (FYH) and fall frost injury (FFI). In this study, we set up a common garden experiment with open pollinated progeny from natural populations of Scots pine (Pinus sylvestris), with one large sample from single population (500 families) and smaller samples from across a latitudinal gradient. BST, FYH and induced FFI were first measured in a greenhouse. The seedlings were then planted in the field, where survival and height were measured at the age of 9 years as fitness proxies. We compared between- and within-population variation and genetic correlations of these three seedling traits, and estimated selection gradients at the family level in our main population, taking into account the potential effects of seed weight. Between-population genetic correlations between seedling traits were high (0.76–0.95). Within-population genetic correlations in the main population were lower (0.14–0.35), as in other populations (0.10–0.39). Within population, extensive adaptive variation persists in the seedling traits, in line with rather weak selection gradients, yet maintaining the clines. Although our sampling does not cover the whole cline equally, the results suggest that the individual clines in these traits are maintained by largely independently acting selection, which results in fewer constraints in adaptation under changing climate.

**Keywords:** budset timing, fall frost injury, first-year height, genetic correlation, local adaptation, selection gradient

## Lay summary

Clines in genetically governed traits result from adaptation to gradually varying environmental gradients. Genetic dependencies between traits can influence adaptation. As the climate warms, northern Scots pine must adapt to new combinations of photoperiod and temperature to time its yearly growth period optimally along the cline. Here, we examined the dependency between adaptive traits clines in budset timing, height, and cold tolerance of first-year Scots pine seedlings. While these clines are highly correlated between populations along the cline, within-population analysis with a large natural population sample revealed much weaker genetic correlations. Further, selection was shown to act mostly independently on each trait. The lack of strong genetic correlations between adaptive traits results in fewer constraints on adaptation.

## Introduction

When environmental conditions change gradually over an area, for example, latitudinally, adaptation can be reflected as clinal patterns where population means of adaptation-related traits track local optima (Felsenstein, 1977). Clinal variation in many traits of plants and animals has been observed throughout centuries of research (Endler, 1977). Where overall genomic variation suggests little or no clinal differentiation, the clines likely have been formed by the action of natural selection (Chen et al., 2014; Kujala & Knürr et al., 2017; Mimura & Aitken, 2007; Savolainen

et al., 2007; Viherä-Aarnio et al., 2005). When populations along the cline are connected through current gene flow, clines reflect a balance between selection and gene flow (Slatkin, 1973, 1978).

Many different traits within a species may display similar clines, raising the question whether selection has acted independently on the traits, or selection on one trait has resulted in correlated clines in other traits through physiological constraints and/or pleiotropy and linkage disequilibrium (indirect selection; Guillaume, 2011; Lande, 1979; Merilä & Björklund, 1999). Kingsolver et al. (2001) did not find much evidence of indirect

<footnote n="1">¹Natural Resources Institute Finland, Oulu, Finland</footnote>

<footnote n="2">²University of Oulu, Ecology and Genetics, Oulu, Finland</footnote>

<footnote n="3">³Université de Strasbourg, INRAE, SVQV UMR-A 1131, Colmar, France</footnote>

<footnote n="4">⁴Natural Resources Institute Finland, Jokioinen, Finland</footnote>

<footnote n="5">⁵Natural Resources Institute Finland, Helsinki, Finland</footnote>
selection. Geber and Griffen (2003), in contrast, concluded that indirect selection effects are common among plant studies. To resolve these issues, it is important that the covariation and the effects of selection are also studied within population. Between- and within-population correlations in many adaptive traits have been reported (see Howe et al., 2003 for review on conifers), but comparisons of these two levels of correlations are needed in natural populations sampled along clines.

Genetic models of clinal variation assume that each population evolves toward a different optimum value of the trait through diversifying selection, while stabilizing selection keeps each local population concentrated around its local optimum (Barton, 1999; Slatkin, 1978). Within population, recombination, gene flow, and mutation bring in new variation that selection acts on (Walsh & Lynch, 2018). Selection gradients for the individual traits—and selection on the correlation itself—can be measured via a multivariate analysis (Arnold et al., 2008; Geyer et al., 2007; Kingsolver & Pfennig, 2007; Lande & Arnold, 1983; Phillips & Arnold, 1989). A prospective study, where individuals of a cohort are characterized for the traits of interest early on, and then monitored for fitness, allows catching part of the lifetime selection.

Scots pine (Pinus sylvestris) is a widespread wind pollinated boreal conifer with a large effective population size (Mirov, 1967), that shows very little differentiation of genomic markers between populations (Pyhäjärvi et al., 2019; Savolainen et al., 2007; Tyrmi et al., 2020). Seedling traits, however, show clinal variation in common garden conditions (Aho, 1994; Kujala & Knürr et al., 2017; Mikola, 1982; Notivol et al., 2007). In many other tree species as well, budset timing (BST) at the end of the first growing season shows a clear latitudinal cline in population means (Alberto et al., 2013; Howe et al., 2003). As fall approaches, northern populations in common gardens set their buds on average earlier (in longer days) than southern populations, an apparent adaptation to avoid cold damage in the fall. A cline is also observed in first-year height (FYH), as southern seedlings grow taller than the northern ones in common gardens (Alberto et al., 2013; Notivol et al., 2007). In natural conditions, the size of the seedling can be an indicator of general vigor and competitive ability in dense populations and among other vegetation (for Scots pine, see Jäghagen & Albrektsen, 1996; Johansson et al., 2015). Fall frost injury (FFI) also shows a clinal pattern (Aho, 1994; Howe et al., 2003; Hurme et al., 1997), and these three traits are highly correlated across divergent populations (Hurme et al., 1997; Mikola, 1980).

The nature of genetic correlations and selection pressures in these traits is highly relevant in a warming climate; constraints due to correlations can influence the ability of trees to adapt via changes in the timing of annual growth (Aitken et al., 2008; Benito Garzón et al., 2011; Bradshaw & Holzapfel, 2001; Ettersson & Shaw, 2001; Oleksyn et al., 1998; Saikkonen et al., 2012; Savolainen et al., 2004). Here we examine the interdependence of clines. We make two distinct predictions: (a) Under the scenario of independently formed clines, we expect to see high between population, but low within-population correlations between traits. Furthermore, selection gradients within population are expected to be significant for each independent trait, and selection for between trait correlation not significant. (b) If parallel clines form as a result of indirect selection, we expect high correlations both between and within populations. For the trait under predominantly indirect selection, the selection gradient is expected to be low. We compared the correlation patterns between and within populations and characterized intrapopulation selection gradients. Although our analysis is mainly from a single large population and further work in other parts of the clines is needed, our results

offer support to the hypothesis that the different seedling clines have not formed solely because of correlated, indirect selection responses.

# Methods

## Material

Our focal population consists of naturally regenerated stands in the Evoltree Intensive Study Site in Punkaharju (EVOLTREE, n.d.), southeastern Finland. Open pollinated seed from 500 trees were sampled in two areas roughly 20 km apart with a minimum distance of 20 m between most pairs of trees to avoid close relatedness. About 380 trees were selected in Ranta-Halola (61°39'N, 29°17'E), and 120 trees in the slightly more northern Mäkrä (61°50'N, 29°23'E) site (Table 1; Figure 1), assuming that two naturally regenerated Scots pine stands 20 km apart represent the same gene pool, and that the trees in both locations have adapted to similar environmental conditions. The age of the mother trees (estimated from the number of tree rings at breast height) varied between 33 and 145 years. For each mother tree, 24 open pollinated seedlings were sown in a greenhouse experiment. We expect open pollinated offspring to be mostly half-sibs (Kärkkäinen & Savolainen, 1993; Muona & Harju, 1989; Squillace, 1974).

In addition, 10 open pollinated families from five natural populations covering latitudes ~50°N–67°N were included, and hereafter called latitudinal populations (Table 1; Figure 1). The number of half-sib seedlings per family in these populations was 12. We settled for a smaller sample size in these latitudinal populations, as the clinal variation has been characterized in many previous studies with larger data.

## Greenhouse experiment and measurements

During June 9th–13th, 2008, a greenhouse experiment was established in Loppi, Finland (60°37'N, 24°25'E; Figure 1). The open pollinated seed of 550 mother trees were sown in disinfected Plantek 49F seedling trays filled with Kekkilä M6 W peat. The greenhouse had ambient temperature and light conditions. Seedlings were watered regularly. The material was fully randomized within six blocks; four seedlings per family for Punkaharju; and two for latitudinal populations were grown within each block. The FYH (in mm) of each seedling was measured with a caliper from the lowest needle to the top of the seedling from August 4th to October 24th until there was no change for three consecutive weeks. BST was recorded as days from sowing to the formation of the terminal bud (scored visually as when the bud was completely visible). Each block was scored for budset and height once a week (blocks 1–5 on consecutive days, block 6 in five parts on 5 consecutive days). Data on seedlings that died already during the greenhouse experiment or showed signs of not being healthy were excluded.

<table><caption>**Table 1.** Origins and sample sizes of the plant material used in this study.</caption><thead><tr><th>Population</th><th>Country</th><th>Latitude</th><th>Longitude</th><th>No. of families</th></tr></thead><tbody><tr><td>PUN</td><td>Finland</td><td>61°39'N–61°50'N</td><td>29°17'E–29°23'E</td><td>500</td></tr><tr><td>FINn</td><td>Finland</td><td>67°11'N</td><td>24°03'E</td><td>10</td></tr><tr><td>FINc</td><td>Finland</td><td>63°45'N</td><td>24°05'E</td><td>10</td></tr><tr><td>FINs</td><td>Finland</td><td>60°52'N</td><td>21°20'E</td><td>10</td></tr><tr><td>SWE</td><td>Sweden</td><td>56°28'N</td><td>15°55'E</td><td>10</td></tr><tr><td>POL</td><td>Poland</td><td>50°41'N</td><td>20°05'E</td><td>10</td></tr></tbody></table>

Note. FINc = Central Finland; FINn = Northern Finland; FINs = Southern Finland; POL = Poland; PUN = Punkaharju; SWE = Sweden.
selection. Geber and Griffen (2003), in contrast, concluded that indirect selection effects are common among plant studies. To resolve these issues, it is important that the covariation and the effects of selection are also studied within population. Between- and within-population correlations in many adaptive traits have been reported (see Howe et al., 2003 for review on conifers), but comparisons of these two levels of correlations are needed in natural populations sampled along clines.

Genetic models of clinal variation assume that each population evolves toward a different optimum value of the trait through diversifying selection, while stabilizing selection keeps each local population concentrated around its local optimum (Barton, 1999; Slatkin, 1978). Within population, recombination, gene flow, and mutation bring in new variation that selection acts on (Walsh & Lynch, 2018). Selection gradients for the individual traits—and selection on the correlation itself—can be measured via a multivariate analysis (Arnold et al., 2008; Geyer et al., 2007; Kingsolver & Pfennig, 2007; Lande & Arnold, 1983; Phillips & Arnold, 1989). A prospective study, where individuals of a cohort are characterized for the traits of interest early on, and then monitored for fitness, allows catching part of the lifetime selection.

Scots pine (Pinus sylvestris) is a widespread wind pollinated boreal conifer with a large effective population size (Mirov, 1967), that shows very little differentiation of genomic markers between populations (Pyhäjärvi et al., 2019; Savolainen et al., 2007; Tyrmi et al., 2020). Seedling traits, however, show clinal variation in common garden conditions (Aho, 1994; Kujala & Knürr et al., 2017; Mikola, 1982; Notivol et al., 2007). In many other tree species as well, budset timing (BST) at the end of the first growing season shows a clear latitudinal cline in population means (Alberto et al., 2013; Howe et al., 2003). As fall approaches, northern populations in common gardens set their buds on average earlier (in longer days) than southern populations, an apparent adaptation to avoid cold damage in the fall. A cline is also observed in first-year height (FYH), as southern seedlings grow taller than the northern ones in common gardens (Alberto et al., 2013; Notivol et al., 2007). In natural conditions, the size of the seedling can be an indicator of general vigor and competitive ability in dense populations and among other vegetation (for Scots pine, see Jäghagen & Albrektson, 1996; Johansson et al., 2015). Fall frost injury (FFI) also shows a clinal pattern (Aho, 1994; Howe et al., 2003; Hurme et al., 1997), and these three traits are highly correlated across divergent populations (Hurme et al., 1997; Mikola, 1980).

The nature of genetic correlations and selection pressures in these traits is highly relevant in a warming climate; constraints due to correlations can influence the ability of trees to adapt via changes in the timing of annual growth (Aitken et al., 2008; Benito Garzón et al., 2011; Bradshaw & Holzapfel, 2001; Ettersson & Shaw, 2001; Oleksyn et al., 1998; Saikkonen et al., 2012; Savolainen et al., 2004). Here we examine the interdependence of clines. We make two distinct predictions: (a) Under the scenario of independently formed clines, we expect to see high between population, but low within-population correlations between traits. Furthermore, selection gradients within population are expected to be significant for each independent trait, and selection for between trait correlation not significant. (b) If parallel clines form as a result of indirect selection, we expect high correlations both between and within populations. For the trait under predominantly indirect selection, the selection gradient is expected to be low. We compared the correlation patterns between and within populations and characterized intrapopulation selection gradients. Although our analysis is mainly from a single large population and further work in other parts of the clines is needed, our results

offer support to the hypothesis that the different seedling clines have not formed solely because of correlated, indirect selection responses.

# Methods

## Material

Our focal population consists of naturally regenerated stands in the Evoltree Intensive Study Site in Punkaharju (EVOLTREE, n.d.), southeastern Finland. Open pollinated seed from 500 trees were sampled in two areas roughly 20 km apart with a minimum distance of 20 m between most pairs of trees to avoid close relatedness. About 380 trees were selected in Ranta-Halola (61°39'N, 29°17'E), and 120 trees in the slightly more northern Mäkrä (61°50'N, 29°23'E) site (Table 1; Figure 1), assuming that two naturally regenerated Scots pine stands 20 km apart represent the same gene pool, and that the trees in both locations have adapted to similar environmental conditions. The age of the mother trees (estimated from the number of tree rings at breast height) varied between 33 and 145 years. For each mother tree, 24 open pollinated seedlings were sown in a greenhouse experiment. We expect open pollinated offspring to be mostly half-sibs (Kärkkäinen & Savolainen, 1993; Muona & Harju, 1989; Squillace, 1974).

In addition, 10 open pollinated families from five natural populations covering latitudes ~50°N–67°N were included, and hereafter called latitudinal populations (Table 1; Figure 1). The number of half-sib seedlings per family in these populations was 12. We settled for a smaller sample size in these latitudinal populations, as the clinal variation has been characterized in many previous studies with larger data.

## Greenhouse experiment and measurements

During June 9th–13th, 2008, a greenhouse experiment was established in Loppi, Finland (60°37'N, 24°25'E; Figure 1). The open pollinated seed of 550 mother trees were sown in disinfected Plantek 49F seedling trays filled with Kekkilä M6 W peat. The greenhouse had ambient temperature and light conditions. Seedlings were watered regularly. The material was fully randomized within six blocks; four seedlings per family for Punkaharju; and two for latitudinal populations were grown within each block. The FYH (in mm) of each seedling was measured with a caliper from the lowest needle to the top of the seedling from August 4th to October 24th until there was no change for three consecutive weeks. BST was recorded as days from sowing to the formation of the terminal bud (scored visually as when the bud was completely visible). Each block was scored for budset and height once a week (blocks 1–5 on consecutive days, block 6 in five parts on 5 consecutive days). Data on seedlings that died already during the greenhouse experiment or showed signs of not being healthy were excluded.

<table><caption>**Table 1.** Origins and sample sizes of the plant material used in this study.</caption><thead><tr><th>Population</th><th>Country</th><th>Latitude</th><th>Longitude</th><th>No. of families</th></tr></thead><tbody><tr><td>PUN</td><td>Finland</td><td>61°39'N–61°50'N</td><td>29°17'E–29°23'E</td><td>500</td></tr><tr><td>FINn</td><td>Finland</td><td>67°11'N</td><td>24°03'E</td><td>10</td></tr><tr><td>FINc</td><td>Finland</td><td>63°45'N</td><td>24°05'E</td><td>10</td></tr><tr><td>FINs</td><td>Finland</td><td>60°52'N</td><td>21°20'E</td><td>10</td></tr><tr><td>SWE</td><td>Sweden</td><td>56°28'N</td><td>15°55'E</td><td>10</td></tr><tr><td>POL</td><td>Poland</td><td>50°41'N</td><td>20°05'E</td><td>10</td></tr></tbody></table>

Note. FINc = Central Finland; FINn = Northern Finland; FINs = Southern Finland; POL = Poland; PUN = Punkaharju; SWE = Sweden.
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Origins of genetic material used in this study, and locations of study sites (greenhouse, G, and field site, F, in inset). Punkaharju population (in inset) consists of two close-by subpopulations in Eastern Finland, Ranta-Halola (PUN<sub>rh</sub>) and Mäkrä (PUN<sub>m</sub>). Latitudinal populations: POL = Poland, SWE = Sweden, FINs = Southern Finland, FINc = Central Finland, FINn = Northern Finland.</figcaption></figure>

To obtain estimates for FFI, a similar greenhouse experiment with the same families and seedlots was repeated in 2010. About 12 seedlings per family from Punkaharju, and six seedlings per family for the latitudinal populations were sown between June 8th and 10th, and grown in a similar manner as in 2008, in three blocks. As the daylength and temperature in the greenhouse decreased (following outdoor conditions), growth cessation and cold acclimation were triggered. In the week of October 4th–8th, the seedlings were taken to a cold chamber. The cold treatment followed the protocol in, for example, Aho (1994). Preliminary testing (data not shown) for finding the exact treatment conditions was first done to ensure that there was variation in freezing damage. Seedlings were first sprayed with distilled water to prevent excessive super-cooling before initiation of cold exposure. The freezing test conditions were: Incubation in +5 °C for 1 hr; cooling at the rate of 3 °C hr⁻¹; freezing at the target temperature of -20° C for 4 hr; warming at the rate of 3 °C hr⁻¹; incubation at +5 °C for 6 hr. After treatment the seedlings were kept in a greenhouse under natural light conditions. Two weeks after the exposure to frost, each seedling was visually scored for damage indicated by chlorophyl breakdown that is, the percentage of browned needles area. Five ordinal classes for FFI assessments were used as follows: 0 = no damage, 1 = 25% damage, 2 = 50% damage, 3 = 75% damage, and 4 = dead. FYH measurements and

BST are not reported here for the 2010 experiment as that data was truncated (seedlings were taken to the cold treatment before all had ceased growth and set bud).

## Field experiment and measurements

We set up a field experiment in which we measured the survival and height of the seedlings at the age of 9 years. Height at this age is not merely a measure of growth capacity, but also of the general health of the tree, perhaps even of future growth and reproductive potential and is often used as fitness proxy (Savolainen et al., 2007; Younginger et al., 2017). Half of the seedlings of the 2008 greenhouse experiment (three blocks) were used to establish this field experiment. All 500 + 50 families were represented. For the first winter, these seedlings were placed in a controlled cold room to overwinter, covered with tarpaulin. After overwintering they were transferred to a nursery garden close to the origin of the focal population, where they were kept for one full growing season. All surviving seedlings were planted in the final study field in Tihvii (61°52'N, 29°17'E) (Figure 1) in early June 2010. The study field covers 2 km² of flat dry boreal forest soil. The seedlings were fully randomized within 10 field blocks and placed roughly 2 m apart. Each family was represented by one seedling within each block. These 10 blocks were surrounded by a minimum of one row of buffer seedlings. Survival in 2010 was recorded at the
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Origins of genetic material used in this study, and locations of study sites (greenhouse, G, and field site, F, in inset). Punkaharju population (in inset) consists of two close-by subpopulations in Eastern Finland, Ranta-Halola (PUN<sub>rh</sub>) and Mäkrä (PUN<sub>m</sub>). Latitudinal populations: POL = Poland, SWE = Sweden, FINs = Southern Finland, FINc = Central Finland, FINn = Northern Finland.</figcaption></figure>

To obtain estimates for FFI, a similar greenhouse experiment with the same families and seedlots was repeated in 2010. About 12 seedlings per family from Punkaharju, and six seedlings per family for the latitudinal populations were sown between June 8th and 10th, and grown in a similar manner as in 2008, in three blocks. As the daylength and temperature in the greenhouse decreased (following outdoor conditions), growth cessation and cold acclimation were triggered. In the week of October 4th–8th, the seedlings were taken to a cold chamber. The cold treatment followed the protocol in, for example, Aho (1994). Preliminary testing (data not shown) for finding the exact treatment conditions was first done to ensure that there was variation in freezing damage. Seedlings were first sprayed with distilled water to prevent excessive super-cooling before initiation of cold exposure. The freezing test conditions were: Incubation in +5 °C for 1 hr; cooling at the rate of 3 °C hr⁻¹; freezing at the target temperature of -20° C for 4 hr; warming at the rate of 3 °C hr⁻¹; incubation at +5 °C for 6 hr. After treatment the seedlings were kept in a greenhouse under natural light conditions. Two weeks after the exposure to frost, each seedling was visually scored for damage indicated by chlorophyl breakdown that is, the percentage of browned needles area. Five ordinal classes for FFI assessments were used as follows: 0 = no damage, 1 = 25% damage, 2 = 50% damage, 3 = 75% damage, and 4 = dead. FYH measurements and

BST are not reported here for the 2010 experiment as that data was truncated (seedlings were taken to the cold treatment before all had ceased growth and set bud).

## Field experiment and measurements

We set up a field experiment in which we measured the survival and height of the seedlings at the age of 9 years. Height at this age is not merely a measure of growth capacity, but also of the general health of the tree, perhaps even of future growth and reproductive potential and is often used as fitness proxy (Savolainen et al., 2007; Younginger et al., 2017). Half of the seedlings of the 2008 greenhouse experiment (three blocks) were used to establish this field experiment. All 500 + 50 families were represented. For the first winter, these seedlings were placed in a controlled cold room to overwinter, covered with tarpaulin. After overwintering they were transferred to a nursery garden close to the origin of the focal population, where they were kept for one full growing season. All surviving seedlings were planted in the final study field in Tihvii (61°52'N, 29°17'E) (Figure 1) in early June 2010. The study field covers 2 km² of flat dry boreal forest soil. The seedlings were fully randomized within 10 field blocks and placed roughly 2 m apart. Each family was represented by one seedling within each block. These 10 blocks were surrounded by a minimum of one row of buffer seedlings. Survival in 2010 was recorded at the
time of planting and subsequent survival check was performed during fall 2011.

The traits measured on the 9-year-old trees in 2017 were survival and tree height. Trees were scored as dead (0) or alive (1). Missing trees were scored as dead. Trees which could not be identified were excluded. Living trees were measured in 5 cm intervals. In branched trees only the longer branch was measured.

## Seedling trait variance, family estimates, heritabilities

Variance in seedling traits was examined both at the level of individual seedlings, and at family level. For Punkaharju, a linear mixed model

$$Y_{ijq} = m + b_i + f_j + e_{ijq} \quad (1)$$

where $Y_{ijq}$ is the trait value of the qth seedling in ith block belonging to family j, m is the overall mean, $b_i$ is the fixed effect of the ith block, $f_j$ is the random effect of the jth family, and $e_{ijq}$ the random error term, was fitted with REML (restricted maximum likelihood) using the function lme in R (R Core Team, 2022, version 4.2.2) in package "nlme" (Pinheiro et al., 2020, version 3.1-148). Family estimates were calculated as family BLUPs (best linear unbiased predictors) from Equation 1, centered around the intercept (mean). Equation 1 was also used for FFI, as the underlying phenotype is continuous in spite of the ordinal scale scoring.

Additive genetic variances ($V_A$) for the traits were estimated from the above model as 4 times the variance attributable to the family effect. Narrow sense heritabilities were calculated as ($V_A$) divided by the total variance. Standard errors for heritabilities were calculated as square root of Dickerson approximation for variance of heritability (Dickerson, 1969). $V_A$ and heritability were estimated only in the Punkaharju population as the latitudinal population samples were considerably smaller. Our estimates of $V_A$ and heritability assume that the progenies consist of half-sibs only, that mother trees are unrelated, and that there are no maternal environmental or genetic effects. Violations of these assumptions could inflate the estimates. However, full sibs are rare in Scots pine and maternal trees of this population have very low relatedness (Harju & Muona, 1989; Niskanen et al. in review). Maternal effects through SW were examined in the Supplementary material.

Student's two-sided t-test was used to test for differences in mean trait values between the subpopulations Ranta-Halola and Mäkrä at the level of family estimates. Effects of early mortality on seedling trait distributions were examined in the Supplementary material.

## Phenotypic and genetic correlations

Phenotypic correlations were calculated as Pearson correlation of the individual seedling values in Punkaharju population. Within-population genetic correlations were estimated in Punkaharju population as Pearson correlation of the family estimates from Equation 1. For the latitudinal populations, a model

$$Y_{ijq} = m + b_i + p_k + f_j + e_{ijq} \quad (2)$$

was fitted. This is similar to Equation 1 except for having population k as an additional fixed factor. Between-population correlations were calculated (a) as Pearson correlations of population wise family estimates $m + f_j + p_k$ (n = 50) and (b) with population means of family estimates (n = 5). We also estimated a second set of within-population correlations using the combined

family estimates (n = 50) of the latitudinal populations, $m + f_j$. Multicollinearity between the traits was checked using the Punkaharju data by computing variance inflation factors with R package "car" (Fox & Weisberg, 2019). Maternal environmental and genetic effects could increase the family component of variation. These genetic correlations can thus not be considered strictly additive.

## Estimation of selection gradients

Estimates of the strength of directional selection (average slope of fitness surface) and quadratic selection (curvature of relative fitness surface near population mean) can be acquired from a multiple regression of relative fitness on the traits of interest (Arnold et al., 2008; Kingsolver & Pfennig, 2007; Lande & Arnold, 1983; Phillips & Arnold, 1989). The method corrects for the possible correlations among the measured traits, and thus selection gradients describe only the direct effects of selection on a particular trait, compared to the selection differential which encompasses both the direct and indirect effects. Here, we analyzed selection gradients in the Punkaharju population at the family level (see Rausher, 1992) with a multiple regression model including all seedling traits (BST, FYH, and FFI) and seed weight (SW). Family estimates from Equation 1 were used as an input for the seedling traits, and the weight of 200 seeds from the respective seedlots as seed weight. These four traits were first standardized by dividing the values with their standard deviations. Absolute fitness for each family was calculated in two ways: (a) as a family estimate of height in 2017 (age 9 at the field site) with Equation 1, also including the dead seedlings with height 0, and (b) as a proportion of surviving seedlings in each family. Relative fitness was then obtained by dividing with mean absolute fitness. Selection analysis was repeated by omitting seedlings that died before planting in the field (Supplementary material).

As our data were not strictly multivariate normal (Madria's test in R package "MVN"; Korkmaz et al., 2014), we calculated the gradients with a two-step process (Lande & Arnold, 1983; Phillips & Arnold, 1989): Estimates for directional selection were first obtained from a model consisting only of the linear terms,

$$w = \alpha + \sum_{i=1}^{4} \beta_i f_i + \epsilon, \quad (3)$$

where w is relative fitness, $f_i$ family estimate of trait i, $\alpha$ an intercept, and $\beta_i$ the directional selection gradient of trait i. The nonlinear selection gradients were then estimated using a model containing also the second-order terms,

$$w = \alpha + \sum_{i=1}^{4} \beta_i^* f_i + \sum_{i=1}^{4} \frac{1}{2} \gamma_{ii} f_i^2 + \sum_{i=1}^{4} \sum_{k>i} \gamma_{ik} f_k + \epsilon, \quad (4)$$

where $\gamma_{ii}$ gives the stabilizing/disruptive selection gradient of trait i, and $\gamma_{ik}$ the correlational selection gradient for traits i and k. In both cases, parameters were estimated using R function lm(). (Note that the function returns the quadratic coefficients as $1/2 \gamma_{ii}$, and hence must be doubled [see Stinchcombe et al., 2008]).

Model in Equation 4 indicated statistically significant correlational selection for some pairs of traits, possibly leading to biased estimates for the stabilizing/disruptive selection gradients. Therefore, we also performed a canonical analysis using mutually independent index traits, following Blows and Brooks (2003); Brooks et al. (2005); Phillips and Arnold (1989); and Simms (1990). The canonical analysis is based on an Eigen decomposition of the symmetric matrix
time of planting and subsequent survival check was performed during fall 2011.

The traits measured on the 9-year-old trees in 2017 were survival and tree height. Trees were scored as dead (0) or alive (1). Missing trees were scored as dead. Trees which could not be identified were excluded. Living trees were measured in 5 cm intervals. In branched trees only the longer branch was measured.

## Seedling trait variance, family estimates, heritabilities

Variance in seedling traits was examined both at the level of individual seedlings, and at family level. For Punkaharju, a linear mixed model

$$Y_{ijq} = m + b_i + f_j + e_{ijq} \quad (1)$$

where $Y_{ijq}$ is the trait value of the qth seedling in ith block belonging to family j, m is the overall mean, $b_i$ is the fixed effect of the ith block, $f_j$ is the random effect of the jth family, and $e_{ijq}$ the random error term, was fitted with REML (restricted maximum likelihood) using the function lme in R (R Core Team, 2022, version 4.2.2) in package “nlme” (Pinheiro et al., 2020, version 3.1-148). Family estimates were calculated as family BLUPs (best linear unbiased predictors) from Equation 1, centered around the intercept (mean). Equation 1 was also used for FFI, as the underlying phenotype is continuous in spite of the ordinal scale scoring.

Additive genetic variances ($V_A$) for the traits were estimated from the above model as 4 times the variance attributable to the family effect. Narrow sense heritabilities were calculated as ($V_A$) divided by the total variance. Standard errors for heritabilities were calculated as square root of Dickerson approximation for variance of heritability (Dickerson, 1969). $V_A$ and heritability were estimated only in the Punkaharju population as the latitudinal population samples were considerably smaller. Our estimates of $V_A$ and heritability assume that the progenies consist of half-sibs only, that mother trees are unrelated, and that there are no maternal environmental or genetic effects. Violations of these assumptions could inflate the estimates. However, full sibs are rare in Scots pine and maternal trees of this population have very low relatedness (Harju & Muona, 1989; Niskanen et al. in review). Maternal effects through SW were examined in the Supplementary material.

Student's two-sided t-test was used to test for differences in mean trait values between the subpopulations Ranta-Halola and Mäkrä at the level of family estimates. Effects of early mortality on seedling trait distributions were examined in the Supplementary material.

## Phenotypic and genetic correlations

Phenotypic correlations were calculated as Pearson correlation of the individual seedling values in Punkaharju population. Within-population genetic correlations were estimated in Punkaharju population as Pearson correlation of the family estimates from Equation 1. For the latitudinal populations, a model

$$Y_{ijq} = m + b_i + p_k + f_j + e_{ijq} \quad (2)$$

was fitted. This is similar to Equation 1 except for having population k as an additional fixed factor. Between-population correlations were calculated (a) as Pearson correlations of population wise family estimates $m + f_j + p_k$ (n = 50) and (b) with population means of family estimates (n = 5). We also estimated a second set of within-population correlations using the combined

family estimates (n = 50) of the latitudinal populations, $m + f_j$. Multicollinearity between the traits was checked using the Punkaharju data by computing variance inflation factors with R package “car” (Fox & Weisberg, 2019). Maternal environmental and genetic effects could increase the family component of variation. These genetic correlations can thus not be considered strictly additive.

## Estimation of selection gradients

Estimates of the strength of directional selection (average slope of fitness surface) and quadratic selection (curvature of relative fitness surface near population mean) can be acquired from a multiple regression of relative fitness on the traits of interest (Arnold et al., 2008; Kingsolver & Pfennig, 2007; Lande & Arnold, 1983; Phillips & Arnold, 1989). The method corrects for the possible correlations among the measured traits, and thus selection gradients describe only the direct effects of selection on a particular trait, compared to the selection differential which encompasses both the direct and indirect effects. Here, we analyzed selection gradients in the Punkaharju population at the family level (see Rausher, 1992) with a multiple regression model including all seedling traits (BST, FYH, and FFI) and seed weight (SW). Family estimates from Equation 1 were used as an input for the seedling traits, and the weight of 200 seeds from the respective seedlots as seed weight. These four traits were first standardized by dividing the values with their standard deviations. Absolute fitness for each family was calculated in two ways: (a) as a family estimate of height in 2017 (age 9 at the field site) with Equation 1, also including the dead seedlings with height 0, and (b) as a proportion of surviving seedlings in each family. Relative fitness was then obtained by dividing with mean absolute fitness. Selection analysis was repeated by omitting seedlings that died before planting in the field (Supplementary material).

As our data were not strictly multivariate normal (Madria's test in R package “MVN”; Korkmaz et al., 2014), we calculated the gradients with a two-step process (Lande & Arnold, 1983; Phillips & Arnold, 1989): Estimates for directional selection were first obtained from a model consisting only of the linear terms,

$$w = \alpha + \sum_{i=1}^{4} \beta_i f_i + \epsilon, \quad (3)$$

where w is relative fitness, $f_i$ family estimate of trait i, $\alpha$ an intercept, and $\beta_i$ the directional selection gradient of trait i. The nonlinear selection gradients were then estimated using a model containing also the second-order terms,

$$w = \alpha + \sum_{i=1}^{4} \beta_i^* f_i + \sum_{i=1}^{4} \frac{1}{2} \gamma_{ii} f_i^2 + \sum_{i=1}^{4} \sum_{k>i} \gamma_{ik} f_k + \epsilon, \quad (4)$$

where $\gamma_{ii}$ gives the stabilizing/disruptive selection gradient of trait i, and $\gamma_{ik}$ the correlational selection gradient for traits i and k. In both cases, parameters were estimated using R function lm(). (Note that the function returns the quadratic coefficients as $1/2 \gamma_{ii}$, and hence must be doubled [see Stinchcombe et al., 2008]).

Model in Equation 4 indicated statistically significant correlational selection for some pairs of traits, possibly leading to biased estimates for the stabilizing/disruptive selection gradients. Therefore, we also performed a canonical analysis using mutually independent index traits, following Blows and Brooks (2003); Brooks et al. (2005); Phillips and Arnold (1989); and Simms (1990). The canonical analysis is based on an Eigen decomposition of the symmetric matrix
$$ \Gamma = \begin{bmatrix} \gamma_{11} & \frac{1}{2}\gamma_{12} & \frac{1}{2}\gamma_{13} & \frac{1}{2}\gamma_{14} \\ \frac{1}{2}\gamma_{12} & \gamma_{22} & \frac{1}{2}\gamma_{23} & \frac{1}{2}\gamma_{24} \\ \frac{1}{2}\gamma_{13} & \frac{1}{2}\gamma_{23} & \gamma_{33} & \frac{1}{2}\gamma_{34} \\ \frac{1}{2}\gamma_{14} & \frac{1}{2}\gamma_{24} & \frac{1}{2}\gamma_{34} & \gamma_{44} \end{bmatrix} $$

consisting of the nonlinear gradients from Equation 4, $\Gamma = M\Lambda M^T$, where $M$ is an orthonormal matrix with the eigenvectors $m_i$ of $\Gamma$ as its columns, and $\Lambda$ a diagonal matrix of the corresponding eigenvalues $\lambda_i$. Now we can write Equation 4 in matrix form, expand with $MM^T = I$, reorganize and rename to obtain

$$ w = \alpha + (f^T M) (M^T \beta^*) + (f^T M) M^T \Gamma M (M^T f) = \alpha + z^T \theta + z^T \Lambda z, \quad (5) $$

where $z = M^T f$ represents values of the new index variable, $\theta = M^T \beta^*$ the vector of directional selection gradients with respect to the index variables and the eigenvalues $\lambda_i$ the stabilizing/disruptive selection gradients for index trait i. The relationship between original and index traits can be ascertained as trait's loadings in the eigenvectors $m_i$. Significance levels for the canonical gradients were obtained by using the index variables in a quadratic regression model similar to Equation 4. The stationary point of the selection surface with respect to the index variables was determined by setting the derivative of Equation 5 to zero and solving for $z^* = -\theta/2\Lambda$.

# Results

## Seedling trait variation

Within-population variation in seedling traits was extensive, especially in Punkaharju population with 500 half-sib families. Population trait means form distinct clines, while population distributions overlapped (Table 2; Figure 2). Differences in trait means between the two Punkaharju subpopulations were very small: Means of Ranta-Halola and Mäkrä differed by 1 day in BST (102.7 and 103.8 days, $p = .00019$), by 1 mm in FYH (60.4 and 61.7 mm, $p = .019$), and by 0.04 in FFI (1.90 and 1.94, $p = .016$), respectively. These differences were very small relative to the within Punkaharju variation and to differences between populations.

## Correlations between early seedling traits

Between-population genetic correlations (with five latitudinal populations) varied between 0.76 and 0.99 (Table 3; Figure 3B). The within-population correlations in Punkaharju were much lower; the genetic correlations ranged from 0.14 to 0.35 (Table 3; Figure 3A). Within-population estimates of genetic correlation in the pooled data of five latitudinal populations were between 0.10 and 0.39 (Table 3). The phenotypic correlation between BST and FYH was 0.23. Low correlations were observed between seedling traits and SW, and between tree age and SW (Supplementary Figure S1). Multicollinearity was low; 1.06 (BST), 1.21 (FYH), 1.17 (FFI), and 1.08 (SW).

<table><caption>**Table 2.** Mean and variances of the seedling traits measured in the greenhouse. For Punkaharju estimates are from model (1). For the latitudinal populations mean and total variance were calculated from the seedling data directly. Note that the proportion of seedlings not forming terminal bud during the greenhouse experiment (i.e., missing data) was largest in POL, followed by SWE, likely biasing the mean and variance of these two populations downward.</caption><thead><tr><th rowspan="2">Pop</th><th colspan="4">BST 2008</th><th colspan="4">FYH 2008</th><th colspan="4">FFI 2010</th></tr><tr><th>n</th><th>mean</th><th>VP</th><th>VA</th><th>n</th><th>Mean</th><th>VP</th><th>VA</th><th>n</th><th>Mean</th><th>VP</th><th>VA</th></tr></thead><tbody><tr><td>PUN</td><td>11,837</td><td>103.0</td><td>83.6</td><td>40.4</td><td>11,861</td><td>60.7</td><td>250.3</td><td>141.5</td><td>5980</td><td>1.91</td><td>0.83</td><td>0.27</td></tr><tr><td>FINn</td><td>117</td><td>86.3</td><td>47.2</td><td>-</td><td>117</td><td>38.4</td><td>191.4</td><td>-</td><td>60</td><td>1.17</td><td>0.55</td><td>-</td></tr><tr><td>FINc</td><td>116</td><td>97.2</td><td>64.0</td><td>-</td><td>116</td><td>36.5</td><td>202.7</td><td>-</td><td>60</td><td>0.95</td><td>0.49</td><td>-</td></tr><tr><td>FINs</td><td>100</td><td>108.7</td><td>83.1</td><td>-</td><td>102</td><td>41.2</td><td>185.5</td><td>-</td><td>56</td><td>1.66</td><td>1.25</td><td>-</td></tr><tr><td>SWE</td><td>100</td><td>113.9</td><td>157.0</td><td>-</td><td>108</td><td>67.2</td><td>228.1</td><td>-</td><td>60</td><td>3.68</td><td>0.22</td><td>-</td></tr><tr><td>POL</td><td>93</td><td>114.3</td><td>210.8</td><td>-</td><td>107</td><td>72.0</td><td>260.5</td><td>-</td><td>60</td><td>3.90</td><td>0.09</td><td>-</td></tr></tbody></table>

Note. FINc = Central Finland; FINn = Northern Finland; FINs = Southern Finland; POL = Poland; PUN = Punkaharju; SWE = Sweden.
BST = budset timing (days); FFI = fall frost injury (score 0-4); FYH = first-year height (mm).
n = number of seedlings; VA = additive genetic variance; VP = total variance.

<figure><img src="image_4.png" /><figcaption>**Figure 2.** Phenotypic variation in seedling traits measured in the greenhouse during the first growing season. Plotted values are the individual seedling measurements. BST and FYH were measured in 2008 trial, FFI in 2010 trial. POL = Poland, SWE = Sweden, FINs = Southern Finland, PUN = Punkaharju, FINc = Central Finland, FINn = Northern Finland. BST = budset timing, FYH = first-year height, FFI = fall frost injury.</figcaption></figure>
$$ \Gamma = \begin{bmatrix} \gamma_{11} & \frac{1}{2}\gamma_{12} & \frac{1}{2}\gamma_{13} & \frac{1}{2}\gamma_{14} \\ \frac{1}{2}\gamma_{12} & \gamma_{22} & \frac{1}{2}\gamma_{23} & \frac{1}{2}\gamma_{24} \\ \frac{1}{2}\gamma_{13} & \frac{1}{2}\gamma_{23} & \gamma_{33} & \frac{1}{2}\gamma_{34} \\ \frac{1}{2}\gamma_{14} & \frac{1}{2}\gamma_{24} & \frac{1}{2}\gamma_{34} & \gamma_{44} \end{bmatrix} $$

consisting of the nonlinear gradients from Equation 4, $\Gamma = M\Lambda M^T$, where $M$ is an orthonormal matrix with the eigenvectors $m_i$ of $\Gamma$ as its columns, and $\Lambda$ a diagonal matrix of the corresponding eigenvalues $\lambda_i$. Now we can write Equation 4 in matrix form, expand with $MM^T = I$, reorganize and rename to obtain

$$ w = \alpha + (f^T M) (M^T \beta^*) + (f^T M) M^T \Gamma M (M^T f) = \alpha + z^T \theta + z^T \Lambda z, \quad (5) $$

where $z = M^T f$ represents values of the new index variable, $\theta = M^T \beta^*$ the vector of directional selection gradients with respect to the index variables and the eigenvalues $\lambda_i$ the stabilizing/disruptive selection gradients for index trait i. The relationship between original and index traits can be ascertained as trait's loadings in the eigenvectors $m_i$. Significance levels for the canonical gradients were obtained by using the index variables in a quadratic regression model similar to Equation 4. The stationary point of the selection surface with respect to the index variables was determined by setting the derivative of Equation 5 to zero and solving for $z^* = -\theta/2\Lambda$.

# Results

## Seedling trait variation

Within-population variation in seedling traits was extensive, especially in Punkaharju population with 500 half-sib families. Population trait means form distinct clines, while population distributions overlapped (Table 2; Figure 2). Differences in trait means between the two Punkaharju subpopulations were very small: Means of Ranta-Halola and Mäkrä differed by 1 day in BST (102.7 and 103.8 days, *p* = .00019), by 1 mm in FYH (60.4 and 61.7 mm, *p* = .019), and by 0.04 in FFI (1.90 and 1.94, *p* = .016), respectively. These differences were very small relative to the within Punkaharju variation and to differences between populations.

## Correlations between early seedling traits

Between-population genetic correlations (with five latitudinal populations) varied between 0.76 and 0.99 (Table 3; Figure 3B). The within-population correlations in Punkaharju were much lower; the genetic correlations ranged from 0.14 to 0.35 (Table 3; Figure 3A). Within-population estimates of genetic correlation in the pooled data of five latitudinal populations were between 0.10 and 0.39 (Table 3). The phenotypic correlation between BST and FYH was 0.23. Low correlations were observed between seedling traits and SW, and between tree age and SW (Supplementary Figure S1). Multicollinearity was low; 1.06 (BST), 1.21 (FYH), 1.17 (FFI), and 1.08 (SW).

<table><caption>**Table 2.** Mean and variances of the seedling traits measured in the greenhouse. For Punkaharju estimates are from model (1). For the latitudinal populations mean and total variance were calculated from the seedling data directly. Note that the proportion of seedlings not forming terminal bud during the greenhouse experiment (i.e., missing data) was largest in POL, followed by SWE, likely biasing the mean and variance of these two populations downward.</caption><thead><tr><th rowspan="2">Pop</th><th colspan="4">BST 2008</th><th colspan="4">FYH 2008</th><th colspan="4">FFI 2010</th></tr><tr><th>n</th><th>mean</th><th>VP</th><th>VA</th><th>n</th><th>Mean</th><th>VP</th><th>VA</th><th>n</th><th>Mean</th><th>VP</th><th>VA</th></tr></thead><tbody><tr><td>PUN</td><td>11,837</td><td>103.0</td><td>83.6</td><td>40.4</td><td>11,861</td><td>60.7</td><td>250.3</td><td>141.5</td><td>5980</td><td>1.91</td><td>0.83</td><td>0.27</td></tr><tr><td>FINn</td><td>117</td><td>86.3</td><td>47.2</td><td>-</td><td>117</td><td>38.4</td><td>191.4</td><td>-</td><td>60</td><td>1.17</td><td>0.55</td><td>-</td></tr><tr><td>FINc</td><td>116</td><td>97.2</td><td>64.0</td><td>-</td><td>116</td><td>36.5</td><td>202.7</td><td>-</td><td>60</td><td>0.95</td><td>0.49</td><td>-</td></tr><tr><td>FINs</td><td>100</td><td>108.7</td><td>83.1</td><td>-</td><td>102</td><td>41.2</td><td>185.5</td><td>-</td><td>56</td><td>1.66</td><td>1.25</td><td>-</td></tr><tr><td>SWE</td><td>100</td><td>113.9</td><td>157.0</td><td>-</td><td>108</td><td>67.2</td><td>228.1</td><td>-</td><td>60</td><td>3.68</td><td>0.22</td><td>-</td></tr><tr><td>POL</td><td>93</td><td>114.3</td><td>210.8</td><td>-</td><td>107</td><td>72.0</td><td>260.5</td><td>-</td><td>60</td><td>3.90</td><td>0.09</td><td>-</td></tr></tbody></table>

Note. FINc = Central Finland; FINn = Northern Finland; FINs = Southern Finland; POL = Poland; PUN = Punkaharju; SWE = Sweden.
BST = budset timing (days); FFI = fall frost injury (score 0-4); FYH = first-year height (mm).
n = number of seedlings; VA = additive genetic variance; VP = total variance.

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Phenotypic variation in seedling traits measured in the greenhouse during the first growing season. Plotted values are the individual seedling measurements. BST and FYH were measured in 2008 trial, FFI in 2010 trial. POL = Poland, SWE = Sweden, FINs = Southern Finland, PUN = Punkaharju, FINc = Central Finland, FINn = Northern Finland. BST = budset timing, FYH = first-year height, FFI = fall frost injury.</figcaption></figure>
## Survival and height in the field

Survival until year 2017 was followed in 4,963 Punkaharju and 295 latitudinal population seedlings. Of the Punkaharju seedlings 93% were still alive in 2010 at the time of field planting, 71% in 2011 inventory, and 64% in 2017. The respective numbers for the latitudinal populations (combined) were 73%, 42%, and 34%. The

<table><caption>**Table 3.** Heritabilities and genetic correlations in seedling traits. Diagonal: Heritabilities within Punkaharju population (in bold, *n* = 500 families). Below diagonal: Genetic correlations within Punkaharju population (*n* = 500 families). Above diagonal: Genetic correlations in latitudinal populations. Between population estimates in italics: 1st line from family estimates (*n* = 50) and 2nd line from population mean of family estimates (*n* = 5). Within-population estimates (overall in five latitudinal populations) in plain text. 95% confidence intervals (for correlations) and SE (for heritabilities) in parenthesis.</caption><thead><tr><th></th><th>BST</th><th>FYH</th><th>FFI</th></tr></thead><tbody><tr><td>BST</td><td><strong>0.48 (0.04)</strong></td><td>0.76 (0.61;0.86)***<br>0.79 (-0.31;0.99)<br>0.39 (0.12;0.60)**</td><td>0.81 (0.68;0.89)***<br>0.83 (-0.21;0.99)<br>0.10 (-0.18;0.37)</td></tr><tr><td>FYH</td><td>0.14 (0.05;0.22)**</td><td><strong>0.57 (0.05)</strong></td><td>0.96 (0.93;0.98)***<br>0.99 (0.92;1.00)***<br>0.39 (0.13;0.60)**</td></tr><tr><td>FFI</td><td>0.19 (0.10;0.27)***</td><td>0.35 (0.27;0.42)***</td><td><strong>0.33 (0.04)</strong></td></tr></tbody></table>

Note. BST = budget timing; FFI = fall frost injury; FYH = first-year height. Significance from the respective Pearson's product moment correlation: \* *p* < .05, \*\* *p* < .01, \*\*\* *p* < .001.

heaviest mortality took place during the first year after planting in the field. The high overall mortality within this year is likely due to random environmental factors that do not correlate strongly with the focal traits (Supplementary Tables S1 and S2).

In 2017, 3,165 Punkaharju seedlings were still alive. Mean height of the survivors at age of 9 years was 2.457 m (min 0.35 m; max 4.2 m). The correlation between FYH and height at 9 years was very weak (*r* = 0.063, *p* = .0005, Punkaharju only). The fitness of latitudinal populations was poorer, as predicted given local adaptation. The more southern populations had poorer survival. The northernmost population had fair survival, but the poor growth of the survivors lowered their combined fitness proxy (Figure 4).

## Selection within population

Selection gradients were calculated at the family level for BST, FYH, FFI, and SW (Table 4) to estimate selection strength and patterns. In the two-step analysis (Lande & Arnold, 1983; Phillips & Arnold, 1989), we observed statistically significant directional selection for greater FYH and stabilizing selection in FFI. Selection for positive genetic correlation between FYH and FFI was statistically significant. Selection for negative correlation between FFI and SW was statistically significant when fitness was estimated as family survival. In analysis where seedlings that died before field planting were omitted, only stabilizing selection for FFI was observed (Supplementary Table S3). In canonical analysis (Table 5), directional multitrait selection was indicated on one of

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Bivariate distributions and regression lines (with coefficient of determination *R*² and associated *p*-values) of seedling traits measured in the greenhouse during the first growing season. BST and FYH were measured in 2008 trial, FFI in 2010 trial. Plotted values are the family estimates (A) within Punkaharju population (*n* = 500 families), (B) in five latitudinal populations (*n* = 50 families total; circles = Poland, down-pointing triangles = Sweden, diamonds = Southern Finland, up-pointing triangles = Central Finland, squares = Northern Finland). BST = budget timing, FYH = first-year height, FFI = fall frost injury.</figcaption></figure>
## Survival and height in the field

Survival until year 2017 was followed in 4,963 Punkaharju and 295 latitudinal population seedlings. Of the Punkaharju seedlings 93% were still alive in 2010 at the time of field planting, 71% in 2011 inventory, and 64% in 2017. The respective numbers for the latitudinal populations (combined) were 73%, 42%, and 34%. The

<table><caption>**Table 3.** Heritabilities and genetic correlations in seedling traits. Diagonal: Heritabilities within Punkaharju population (in bold, *n* = 500 families). Below diagonal: Genetic correlations within Punkaharju population (*n* = 500 families). Above diagonal: Genetic correlations in latitudinal populations. Between population estimates in italics: 1st line from family estimates (*n* = 50) and 2nd line from population mean of family estimates (*n* = 5). Within-population estimates (overall in five latitudinal populations) in plain text. 95% confidence intervals (for correlations) and SE (for heritabilities) in parenthesis.</caption><thead><tr><th></th><th>BST</th><th>FYH</th><th>FFI</th></tr></thead><tbody><tr><th>BST</th><td><strong>0.48 (0.04)</strong></td><td>0.76 (0.61;0.86)***<br>0.79 (-0.31;0.99)<br>0.39 (0.12;0.60)**</td><td>0.81 (0.68;0.89)***<br>0.83 (-0.21;0.99)<br>0.10 (-0.18;0.37)</td></tr><tr><th>FYH</th><td>0.14 (0.05;0.22)**</td><td><strong>0.57 (0.05)</strong></td><td>0.96 (0.93;0.98)***<br>0.99 (0.92;1.00)***<br>0.39 (0.13;0.60)**</td></tr><tr><th>FFI</th><td>0.19 (0.10;0.27)***</td><td>0.35 (0.27;0.42)***</td><td><strong>0.33 (0.04)</strong></td></tr></tbody></table>

Note. BST = budset timing; FFI = fall frost injury; FYH = first-year height. Significance from the respective Pearson's product moment correlation: \* *p* < .05, \*\* *p* < .01, \*\*\* *p* < .001.

heaviest mortality took place during the first year after planting in the field. The high overall mortality within this year is likely due to random environmental factors that do not correlate strongly with the focal traits (Supplementary Tables S1 and S2).

In 2017, 3,165 Punkaharju seedlings were still alive. Mean height of the survivors at age of 9 years was 2.457 m (min 0.35 m; max 4.2 m). The correlation between FYH and height at 9 years was very weak (*r* = 0.063, *p* = .0005, Punkaharju only). The fitness of latitudinal populations was poorer, as predicted given local adaptation. The more southern populations had poorer survival. The northernmost population had fair survival, but the poor growth of the survivors lowered their combined fitness proxy (Figure 4).

## Selection within population

Selection gradients were calculated at the family level for BST, FYH, FFI, and SW (Table 4) to estimate selection strength and patterns. In the two-step analysis (Lande & Arnold, 1983; Phillips & Arnold, 1989), we observed statistically significant directional selection for greater FYH and stabilizing selection in FFI. Selection for positive genetic correlation between FYH and FFI was statistically significant. Selection for negative correlation between FFI and SW was statistically significant when fitness was estimated as family survival. In analysis where seedlings that died before field planting were omitted, only stabilizing selection for FFI was observed (Supplementary Table S3). In canonical analysis (Table 5), directional multitrait selection was indicated on one of

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Bivariate distributions and regression lines (with coefficient of determination *R*² and associated *p*-values) of seedling traits measured in the greenhouse during the first growing season. BST and FYH were measured in 2008 trial, FFI in 2010 trial. Plotted values are the family estimates (A) within Punkaharju population (*n* = 500 families), (B) in five latitudinal populations (*n* = 50 families total; circles = Poland, down-pointing triangles = Sweden, diamonds = Southern Finland, up-pointing triangles = Central Finland, squares = Northern Finland). BST = budset timing, FYH = first-year height, FFI = fall frost injury.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Variation in family means of fitness components measured at the field site at the age of 9 years. Family mean height is calculated only among seedlings alive in 2017. Family absolute fitness is the family mean height including also dead seedlings. POL = Poland, SWE = Sweden, FINs = Southern Finland, PUN = Punkaharju, FINc = Central Finland, FINn = Northern Finland.</figcaption></figure>

<table><caption>**Table 4.** Standardized linear selection gradients (β) from the linear model and standardized quadratic and correlational selection gradients (γ) from the full quadratic model. β represents the strength of directional selection on the traits. Diagonal γ (in bold) describes the strength of stabilizing/disruptive selection on the traits. Off-diagonal γ expresses the selection on the pairwise correlation between traits. Above: gradients when family fitness was defined through both survival and height at field. Below: gradients with mean survival of the family as the fitness.</caption><thead><tr><th rowspan="2"></th><th rowspan="2">β</th><th colspan="4">γ</th></tr><tr><th>BST</th><th>FYH</th><th>FFI</th><th>SW</th></tr></thead><tbody><tr><td>BST</td><td>-1.50e-03</td><td><strong>-2.11e-04</strong></td><td>-1.48e-03</td><td>1.13e-03</td><td>2.06e-03</td></tr><tr><td></td><td>-1.50e-02</td><td><strong>-1.09e-02</strong></td><td>-1.45e-02</td><td>1.28e-02</td><td>1.23e-02</td></tr><tr><td>FYH</td><td>4.21e-03**</td><td><strong>-1.08e-03</strong></td><td>3.34e-03*</td><td>-3.37e-04</td><td></td></tr><tr><td></td><td>2.95e-02*</td><td><strong>-2.62e-02</strong></td><td>3.80e-02**</td><td>-2.00e-03</td><td></td></tr><tr><td>FFI</td><td>-8.58e-04</td><td></td><td><strong>-5.68e-03**</strong></td><td>-2.80e-03</td><td></td></tr><tr><td></td><td>-7.39e-03</td><td></td><td><strong>-5.39e-02**</strong></td><td>-2.88e-02*</td><td></td></tr><tr><td>SW</td><td>5.04e-05</td><td></td><td></td><td><strong>8.38e-06</strong></td><td></td></tr><tr><td></td><td>4.61e-03</td><td></td><td></td><td><strong>-1.02e-02</strong></td><td></td></tr></tbody></table>

Note. BST = budset timing; FFI = fall frost injury; FYH = first-year height; SW = seed weight.

Significance from the respective multiple regression model: *p* < .10, \* *p* < .05, \*\* *p* < .01, \*\*\* *p* < .001.

the rotated axes, dominated by FYH and SW (m₃), and stabilizing multitrait selection in another (m₄), for which FFI contributed the most. m₂ θ was marginally significant. The stationary point z* is within our observed range. The gradients were however very small. At the individual seedling level (with only BST and FYH, Supplementary Table S4; Supplementary Figure S2), directional selection was observed for earlier BST and greater FYH, combined with stabilizing selection for FYH.

# Discussion

## High between-population and modest within-population correlations in seedling traits suggest mostly independent clines

Our data showed similar trend in seedling trait clines as observed in many previous studies; southern Scots pine seedlings set bud later, grew taller, and expressed more frost injuries than the northern seedlings. Concordant with the scenario of mostly independent clines, between-population correlations over five latitudinal populations were high, while the within-population correlations were lower, both in our focal Punkaharju sample, and in the pooled data of latitudinal populations.

Agreeing with our results on correlation between BST and FYH, Notivol et al. (2007) found within-population genetic correlation of 0.07 in a different set of populations of P. sylvestris across Europe. Mikola (1980) reported between-population correlation of 0.86 in a geographically wide sample of Scots pine populations. In Skrøppa et al. (2023), the correlation between budset forming percentage and FYH varied between -0.33 and 0.03 within Norway spruce populations. In Howe et al. (2003), see Table 3 with many tree species, within-population correlations between BST and growth range between 0.19 and 0.77 and between 0.59 and 0.86 across populations.

Hurme et al. (1997) observed very high between-population correlation between BST and frost injuries (up to 0.97) in Scots pine. Hurme et al. (2000) also found mostly independent quantitative trait loci (QTLs) for frost hardiness and BST, supporting the view of genetic independence of the traits. Within a single northern Scots pine population in a similarly designed study, Savolainen et al. (2004) reported a genetic correlation of 0.57 (SD 0.07). Our estimate was much smaller, raising the possibility that this correlation is more important in northern populations.

Generally, high between population correlations have been observed between growth and frost damage in other tree species. Within-population correlations vary from negative to positive (see Table 3 in Howe et al., 2003). In Abrahamsson et al. (2012) genetic correlation in open pollinated progeny was 0.23 (±0.20). In Pinus contorta, MacLachlan et al. (2021) found that many single nucleotide polymorphisms (SNPs) associated with seedling height also affected cold tolerance and phenology traits, indicating significant genetic correlation. Their study included trees from a very large geographical and environmental range.

## Field performance is consistent with local adaptation through mostly independent selection

There is ample evidence of local adaptation (Kawecki & Ebert, 2004) in Scots pine from previous studies (e.g., Eiche, 1966; Savolainen et al., 2007). Although we compared the fitness of different populations only at one site, our results were consistent with local adaptation. For both survival and height, the local population outperformed the nonlocal populations. The poor fitness of the southern populations was dominated by low survival, likely due to lower cold tolerance and later cessation of growth (Supplementary Tables S1 and S2). Northernmost population had high survival, but poorer height.

The Lande-Arnold analysis showed statistically significant directional selection for greater FYH and stabilizing selection for
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Variation in family means of fitness components measured at the field site at the age of 9 years. Family mean height is calculated only among seedlings alive in 2017. Family absolute fitness is the family mean height including also dead seedlings. POL = Poland, SWE = Sweden, FINs = Southern Finland, PUN = Punkaharju, FINc = Central Finland, FINn = Northern Finland.</figcaption></figure>

<table><caption>**Table 4.** Standardized linear selection gradients (β) from the linear model and standardized quadratic and correlational selection gradients (γ) from the full quadratic model. β represents the strength of directional selection on the traits. Diagonal γ (in bold) describes the strength of stabilizing/disruptive selection on the traits. Off-diagonal γ expresses the selection on the pairwise correlation between traits. Above: gradients when family fitness was defined through both survival and height at field. Below: gradients with mean survival of the family as the fitness.</caption><thead><tr><th rowspan="2"></th><th rowspan="2">β</th><th colspan="4">γ</th></tr><tr><th>BST</th><th>FYH</th><th>FFI</th><th>SW</th></tr></thead><tbody><tr><td>BST</td><td>-1.50e-03</td><td><strong>-2.11e-04</strong></td><td>-1.48e-03</td><td>1.13e-03</td><td>2.06e-03</td></tr><tr><td></td><td>-1.50e-02</td><td><strong>-1.09e-02</strong></td><td>-1.45e-02</td><td>1.28e-02</td><td>1.23e-02</td></tr><tr><td>FYH</td><td>4.21e-03**</td><td><strong>-1.08e-03</strong></td><td>3.34e-03*</td><td>-3.37e-04</td><td></td></tr><tr><td></td><td>2.95e-02*</td><td><strong>-2.62e-02</strong></td><td>3.80e-02**</td><td>-2.00e-03</td><td></td></tr><tr><td>FFI</td><td>-8.58e-04</td><td></td><td><strong>-5.68e-03**</strong></td><td>-2.80e-03</td><td></td></tr><tr><td></td><td>-7.39e-03</td><td></td><td><strong>-5.39e-02**</strong></td><td>-2.88e-02*</td><td></td></tr><tr><td>SW</td><td>5.04e-05</td><td></td><td></td><td><strong>8.38e-06</strong></td><td></td></tr><tr><td></td><td>4.61e-03</td><td></td><td></td><td><strong>-1.02e-02</strong></td><td></td></tr></tbody></table>

Note. BST = budset timing; FFI = fall frost injury; FYH = first-year height; SW = seed weight.

Significance from the respective multiple regression model: *p* < .10, \* *p* < .05, \*\* *p* < .01, \*\*\* *p* < .001.

the rotated axes, dominated by FYH and SW ($m_3$), and stabilizing multitrait selection in another ($m_4$), for which FFI contributed the most. $m_2 \theta$ was marginally significant. The stationary point $z^*$ is within our observed range. The gradients were however very small. At the individual seedling level (with only BST and FYH, Supplementary Table S4; Supplementary Figure S2), directional selection was observed for earlier BST and greater FYH, combined with stabilizing selection for FYH.

## Discussion

### High between-population and modest within-population correlations in seedling traits suggest mostly independent clines

Our data showed similar trend in seedling trait clines as observed in many previous studies; southern Scots pine seedlings set bud later, grew taller, and expressed more frost injuries than the northern seedlings. Concordant with the scenario of mostly independent clines, between-population correlations over five latitudinal populations were high, while the within-population correlations were lower, both in our focal Punkaharju sample, and in the pooled data of latitudinal populations.

Agreeing with our results on correlation between BST and FYH, Notivol et al. (2007) found within-population genetic correlation of 0.07 in a different set of populations of *P. sylvestris* across Europe. Mikola (1980) reported between-population correlation of 0.86 in a geographically wide sample of Scots pine populations. In Skrøppa et al. (2023), the correlation between budset forming percentage and FYH varied between -0.33 and 0.03 within Norway spruce populations. In Howe et al. (2003), see Table 3 with many tree species, within-population correlations between BST and growth range between 0.19 and 0.77 and between 0.59 and 0.86 across populations.

Hurme et al. (1997) observed very high between-population correlation between BST and frost injuries (up to 0.97) in Scots pine. Hurme et al. (2000) also found mostly independent quantitative trait loci (QTLs) for frost hardiness and BST, supporting the view of genetic independence of the traits. Within a single northern Scots pine population in a similarly designed study, Savolainen et al. (2004) reported a genetic correlation of 0.57 (SD 0.07). Our estimate was much smaller, raising the possibility that this correlation is more important in northern populations.

Generally, high between population correlations have been observed between growth and frost damage in other tree species. Within-population correlations vary from negative to positive (see Table 3 in Howe et al., 2003). In Abrahamsson et al. (2012) genetic correlation in open pollinated progeny was 0.23 (±0.20). In Pinus contorta, MacLachlan et al. (2021) found that many single nucleotide polymorphisms (SNPs) associated with seedling height also affected cold tolerance and phenology traits, indicating significant genetic correlation. Their study included trees from a very large geographical and environmental range.

### Field performance is consistent with local adaptation through mostly independent selection

There is ample evidence of local adaptation (Kawecki & Ebert, 2004) in Scots pine from previous studies (e.g., Eiche, 1966; Savolainen et al., 2007). Although we compared the fitness of different populations only at one site, our results were consistent with local adaptation. For both survival and height, the local population outperformed the nonlocal populations. The poor fitness of the southern populations was dominated by low survival, likely due to lower cold tolerance and later cessation of growth (Supplementary Tables S1 and S2). Northernmost population had high survival, but poorer height.

The Lande-Arnold analysis showed statistically significant directional selection for greater FYH and stabilizing selection for
<table><caption>**Table 5.** Results of canonical analysis. **M** = matrix of eigenvectors from the canonical analysis, expressing the contribution of each original trait to eigenvectors. **m**₁–**m**₄ = eigenvectors of the response surface. **θ** and **λ** = the linear and quadratic selection gradients along the eigenvectors. **z*** describes the stationary point in the rotated coordinate system, min and max the observed range. Above: family fitness was defined through both survival and height at field. Below: mean survival of the family as the fitness.</caption><thead><tr><th rowspan="2"></th><th colspan="4">M</th><th colspan="5">Selection</th></tr><tr><th>BST</th><th>FYH</th><th>FFI</th><th>SW</th><th>θ</th><th>λ</th><th>z*</th><th>min</th><th>max</th></tr></thead><tbody><tr><td>m<sub>1</sub></td><td>0.495</td><td>-0.474</td><td>-0.322</td><td>0.653</td><td>-2.45e-03</td><td>1.60e-03</td><td>1.954e-06</td><td>-2.378</td><td>2.929</td></tr><tr><td></td><td>0.316</td><td>-0.526</td><td>-0.470</td><td>0.635</td><td>-1.19e-02</td><td>9.44e-03</td><td>5.616e-05</td><td>-2.648</td><td>3.338</td></tr><tr><td>m<sub>2</sub></td><td>0.667</td><td>0.586</td><td>0.438</td><td>0.136</td><td>9.31e-04</td><td>-1.76e-04</td><td>8.196e-08</td><td>-3.847</td><td>3.193</td></tr><tr><td></td><td>0.912</td><td>0.097</td><td>0.389</td><td>-0.086</td><td>-1.89e-02</td><td>-4.09e-03</td><td>-3.855e-05</td><td>-3.383</td><td>3.400</td></tr><tr><td>m<sub>3</sub></td><td>-0.494</td><td>0.522</td><td>-0.156</td><td>0.677</td><td>3.63e-03**</td><td>-5.56e-04</td><td>1.009e-06</td><td>-2.571</td><td>3.538</td></tr><tr><td></td><td>-0.067</td><td>0.688</td><td>0.142</td><td>0.709</td><td>2.91e-02**</td><td>-9.52e-03</td><td>1.385e-04</td><td>-3.139</td><td>4.193</td></tr><tr><td>m<sub>4</sub></td><td>-0.255</td><td>-0.398</td><td>0.825</td><td>0.310</td><td>-2.56e-03</td><td>-4.34e-03**</td><td>-5.554e-06</td><td>-2.706</td><td>2.155</td></tr><tr><td></td><td>-0.252</td><td>-0.490</td><td>0.780</td><td>0.297</td><td>-2.10e-02</td><td>-4.65e-02***</td><td>-4.876e-04</td><td>-2.737</td><td>2.042</td></tr></tbody></table>

Note. BST = budget timing; FFI = fall frost injury; FYH = first-year height; SW = seed weight. Significance from the respective multiple regression model with eigenvectors: *p* < .10, \* *p* < .05, \*\* *p* < .01, \*\*\* *p* < .001.

FFI, although with very modest gradients, and some indications for correlated selection. A closer look at the trait combinations with canonical analysis revealed directional selection on the rotated axis **m**₃ reflecting the higher fitness of tall seedlings, to which large seed weight also contributed. Stabilizing selection was observed on **m**₄, where FFI contributes the most.

Weaker selection gradients from reanalysis without seedlings that died already before field planting highlighted the importance of early selection (Supplementary Table S3). Luoranen et al. (2018) showed that most of the damage observed after 1 year after planting was due to frost. General vigor of seedlings and interactions with other factors such as drought may also contribute. Early selection is important also in *P. contorta*, related to drought (Warwell & Shaw, 2019). However, there is also much early random mortality (Lönnroth, 1925; Sarvas, 1962).

Among plant functional traits, there is often selection for earlier phenology (toward earlier and shorter flowering and rapid germination) and for increased size (Caruso et al., 2020). Liu and El-Kassaby (2019) observed selection for later BST in *Populus* tri-chocarpa. At the family level we did not observe clear signs of selection for BST. In the canonical analysis **m**₂—to which BST had the highest contribution—had a marginally significant **θ**. At the individual seedling level Supplementary Table S4; Supplementary Figure S2) earlier budget was favored but cold damage data were not available for the same seedlings. Thus, the measured selection on BST might have been due to the low correlation with FFI. For FYH, genetic- and phenotypic-level analyses agree on the directional selection, although the selection seemed much stronger at the phenotypic level, as the range of variation was much larger than for family estimates.

Although our selection analysis showed statistically significant multitrait selection, the coefficients were small. This is in line with the rather low within-population genetic correlation we observe; strong selection favoring a specific combination of traits would be expected to strengthen the genetic correlation (Roff & Fairbairn, 2012; Svensson et al., 2021).

## Modest selection allows maintenance of large within-population variation within a latitudinal cline

In genetic models of clinal variation, stabilizing selection operates within populations to keep them close to their local optima (Barton, 1999; Slatkin, 1978). The lack of stabilizing selection on BST and FYH here contrasts with this expectation. However, 9-year-old trees have experienced only a part of the selection as selection pressures can vary within one generation (e.g., Grant &

Grant, 1995; Kelly, 1992; Siepielski et al., 2009). Selection due to climate acts early (Haig et al., 1941; Leck et al., 2008), which could favor cold adapted genotypes, whereas competition, possibly favoring taller trees, acts later when the trees are large enough to interfere with each other, for example, by shading or nutrient competition. In a field study with higher density, canopy closure took place at age of 11 and lead to increased mortality in Scots pine (Koelewijn et al., 1999, see also Egbäck et al., 2012). Observing larger selection gradients when using only family mean survival as fitness measure supports this view.

Earlier studies agree with our data showing that despite the pronounced cline, variation in seedling traits exists within Scots pine populations (Hurme et al., 1997; Kujala et al., 2017; Mikola, 1982; Notivol et al., 2007; Savolainen et al., 2004). Pollen-driven gene flow and recombination are expected to generate the extensive variation. Within-population selection cannot be strong relative to gene flow, congruent with fairly high heritability estimates in this and previous studies (Wheelwright et al., 2014). Nevertheless, earlier reciprocal transplant experiments show that even rather short north/south transfers result in fitness reduction (Berlin et al., 2016; Eiche, 1966; Savolainen et al., 2007). With different optima along the cline, strong within-population stabilizing selection is not a necessary condition for high between-population differentiation ($Q_{ST}$) and cline formation even under strong gene flow (see simulation work of Le Corre & Kremer, 2003, Figure 6).

## Parallel clines in changing climate

In northern latitudes, the climate is warming rapidly. The growing season in Scandinavia has increased by ~1 week between 1951 and 2000 and is expected to increase by 20–30 days within few decades in most areas in Finland and in northern and middle Sweden (Bärring et al., 2017; Linderholm et al., 2008). However, photoperiodic conditions will remain the same. As seedling traits in northern populations of Scots pine are adapted to local combinations of photoperiod and temperature (Koski & Sievänen, 1985; Oleksyn et al., 1992, 1998; Olsen, 2010; Vaartaja, 1959), a genetic change will be needed for growth to continue longer while ensuring correct timing of dormancy (Hänninen et al., 2001; Saikkonen et al., 2012; Savolainen et al., 2004, 2011; Singh et al., 2017).

Genetic correlations can either reinforce or impede evolutionary responses, depending on the direction of correlation in relation to selection (Lande, 1979). In the legume *Chamaecrista fasciculata*, correlations likely will constrain adaptation to climate change (Etterson & Shaw, 2001). MacLachlan et al. (2021)
<table><caption>**Table 5.** Results of canonical analysis. **M** = matrix of eigenvectors from the canonical analysis, expressing the contribution of each original trait to eigenvectors. **m**<sub>1</sub>–**m**<sub>4</sub> = eigenvectors of the response surface. **θ** and **λ** = the linear and quadratic selection gradients along the eigenvectors. **z**<sup>*</sup> describes the stationary point in the rotated coordinate system, min and max the observed range. Above: family fitness was defined through both survival and height at field. Below: mean survival of the family as the fitness.</caption><thead><tr><th rowspan="2"></th><th colspan="4">M</th><th colspan="5">Selection</th></tr><tr><th>BST</th><th>FYH</th><th>FFI</th><th>SW</th><th>θ</th><th>λ</th><th>z*</th><th>min</th><th>max</th></tr></thead><tbody><tr><td>m<sub>1</sub></td><td>0.495</td><td>-0.474</td><td>-0.322</td><td>0.653</td><td>-2.45e-03</td><td>1.60e-03</td><td>1.954e-06</td><td>-2.378</td><td>2.929</td></tr><tr><td></td><td>0.316</td><td>-0.526</td><td>-0.470</td><td>0.635</td><td>-1.19e-02</td><td>9.44e-03</td><td>5.616e-05</td><td>-2.648</td><td>3.338</td></tr><tr><td>m<sub>2</sub></td><td>0.667</td><td>0.586</td><td>0.438</td><td>0.136</td><td>9.31e-04</td><td>-1.76e-04</td><td>8.196e-08</td><td>-3.847</td><td>3.193</td></tr><tr><td></td><td>0.912</td><td>0.097</td><td>0.389</td><td>-0.086</td><td>-1.89e-02</td><td>-4.09e-03</td><td>-3.855e-05</td><td>-3.383</td><td>3.400</td></tr><tr><td>m<sub>3</sub></td><td>-0.494</td><td>0.522</td><td>-0.156</td><td>0.677</td><td>3.63e-03**</td><td>-5.56e-04</td><td>1.009e-06</td><td>-2.571</td><td>3.538</td></tr><tr><td></td><td>-0.067</td><td>0.688</td><td>0.142</td><td>0.709</td><td>2.91e-02**</td><td>-9.52e-03</td><td>1.385e-04</td><td>-3.139</td><td>4.193</td></tr><tr><td>m<sub>4</sub></td><td>-0.255</td><td>-0.398</td><td>0.825</td><td>0.310</td><td>-2.56e-03</td><td>-4.34e-03**</td><td>-5.554e-06</td><td>-2.706</td><td>2.155</td></tr><tr><td></td><td>-0.252</td><td>-0.490</td><td>0.780</td><td>0.297</td><td>-2.10e-02</td><td>-4.65e-02***</td><td>-4.876e-04</td><td>-2.737</td><td>2.042</td></tr></tbody></table>

Note. BST = budget timing; FFI = fall frost injury; FYH = first-year height; SW = seed weight. Significance from the respective multiple regression model with eigenvectors: *p* < .10, \* *p* < .05, \*\* *p* < .01, \*\*\* *p* < .001.

FFI, although with very modest gradients, and some indications for correlated selection. A closer look at the trait combinations with canonical analysis revealed directional selection on the rotated axis **m**<sub>3</sub> reflecting the higher fitness of tall seedlings, to which large seed weight also contributed. Stabilizing selection was observed on **m**<sub>4</sub>, where FFI contributes the most.

Weaker selection gradients from reanalysis without seedlings that died already before field planting highlighted the importance of early selection (Supplementary Table S3). Luoranen et al. (2018) showed that most of the damage observed after 1 year after planting was due to frost. General vigor of seedlings and interactions with other factors such as drought may also contribute. Early selection is important also in *P. contorta*, related to drought (Warwell & Shaw, 2019). However, there is also much early random mortality (Lönnroth, 1925; Sarvas, 1962).

Among plant functional traits, there is often selection for earlier phenology (toward earlier and shorter flowering and rapid germination) and for increased size (Caruso et al., 2020). Liu and El-Kassaby (2019) observed selection for later BST in *Populus* tri-chocarpa. At the family level we did not observe clear signs of selection for BST. In the canonical analysis **m**<sub>2</sub>—to which BST had the highest contribution—had a marginally significant **θ**. At the individual seedling level Supplementary Table S4; Supplementary Figure S2) earlier budget was favored but cold damage data were not available for the same seedlings. Thus, the measured selection on BST might have been due to the low correlation with FFI. For FYH, genetic- and phenotypic-level analyses agree on the directional selection, although the selection seemed much stronger at the phenotypic level, as the range of variation was much larger than for family estimates.

Although our selection analysis showed statistically significant multitrait selection, the coefficients were small. This is in line with the rather low within-population genetic correlation we observe; strong selection favoring a specific combination of traits would be expected to strengthen the genetic correlation (Roff & Fairbairn, 2012; Svensson et al., 2021).

## Modest selection allows maintenance of large within-population variation within a latitudinal cline

In genetic models of clinal variation, stabilizing selection operates within populations to keep them close to their local optima (Barton, 1999; Slatkin, 1978). The lack of stabilizing selection on BST and FYH here contrasts with this expectation. However, 9-year-old trees have experienced only a part of the selection as selection pressures can vary within one generation (e.g., Grant &

Grant, 1995; Kelly, 1992; Siepielski et al., 2009). Selection due to climate acts early (Haig et al., 1941; Leck et al., 2008), which could favor cold adapted genotypes, whereas competition, possibly favoring taller trees, acts later when the trees are large enough to interfere with each other, for example, by shading or nutrient competition. In a field study with higher density, canopy closure took place at age of 11 and lead to increased mortality in Scots pine (Koelewijn et al., 1999, see also Egbäck et al., 2012). Observing larger selection gradients when using only family mean survival as fitness measure supports this view.

Earlier studies agree with our data showing that despite the pronounced cline, variation in seedling traits exists within Scots pine populations (Hurme et al., 1997; Kujala et al., 2017; Mikola, 1982; Notivol et al., 2007; Savolainen et al., 2004). Pollen-driven gene flow and recombination are expected to generate the extensive variation. Within-population selection cannot be strong relative to gene flow, congruent with fairly high heritability estimates in this and previous studies (Wheelwright et al., 2014). Nevertheless, earlier reciprocal transplant experiments show that even rather short north/south transfers result in fitness reduction (Berlin et al., 2016; Eiche, 1966; Savolainen et al., 2007). With different optima along the cline, strong within-population stabilizing selection is not a necessary condition for high between-population differentiation ($Q_{ST}$) and cline formation even under strong gene flow (see simulation work of Le Corre & Kremer, 2003, Figure 6).

## Parallel clines in changing climate

In northern latitudes, the climate is warming rapidly. The growing season in Scandinavia has increased by ~1 week between 1951 and 2000 and is expected to increase by 20–30 days within few decades in most areas in Finland and in northern and middle Sweden (Bärring et al., 2017; Linderholm et al., 2008). However, photoperiodic conditions will remain the same. As seedling traits in northern populations of Scots pine are adapted to local combinations of photoperiod and temperature (Koski & Sievänen, 1985; Oleksyn et al., 1992, 1998; Olsen, 2010; Vaartaja, 1959), a genetic change will be needed for growth to continue longer while ensuring correct timing of dormancy (Hänninen et al., 2001; Saikkonen et al., 2012; Savolainen et al., 2004, 2011; Singh et al., 2017).

Genetic correlations can either reinforce or impede evolutionary responses, depending on the direction of correlation in relation to selection (Lande, 1979). In the legume *Chamaecrista fasciculata*, correlations likely will constrain adaptation to climate change (Etterson & Shaw, 2001). MacLachlan et al. (2021)
suggested that shared variation between adaptive traits in *P. contorta* will result in correlated changes. In Scots pine with largely independent clines, adaptation can be flexible under changing conditions (e.g., de La Mata et al., 2017). Furthermore, high heritabilities in adaptive traits, combined with southern gene flow, will contribute to the adaptive capacity. Similar amounts of genetic variation likely exist in other parts of the clines, too, with a potential exception of range edges (Eckert et al., 2008; Kujala & Knürr et al., 2017).

This analysis concentrated on the capacity for genetic adaptation. To predict future responses, it will be necessary to also consider phenotypic plasticity (e.g., Benito Garzón et al., 2011; Hämälä et al., 2022), effects of demography, such as high survival of old trees limiting establishment possibilities (Kuparinen et al., 2010; Savolainen, 2004), and competition and other interactions with other species. Here our main findings come from a single population, and more empirical data should be obtained from other populations to verify the genetic independence of the seedling traits in other parts of the cline. However, this highly powered analysis of one large population already shows that correlated selection alone has not generated the parallel clines.

Current clines show that historically, genetic adaptation has been a very important component of responses to environmental changes (Davis & Shaw, 2001). Increasing information on the genetics of local adaptation in natural populations of forest trees will enable mitigation of climate change effects in artificially regenerated forests, too, via breeding for new conditions and determining use of regeneration material for changing climate (Aitken et al., 2008; Berlin et al., 2016; Cortés et al., 2020; MacLachlan et al., 2017).

## Supplementary material

Supplementary material is available online at *Evolution Letters*.

## Data and code availability

Data and R commands have been deposited to DRYAD (https://doi.org/10.5061/dryad.vmcvdnd09).

## Author contributions

S.T.K. was the main responsible author, aided by contributions from all authors. S.T.K. and H.K. performed the final quantitative and selection analyses. S.T.K. was responsible for interpretation of the results with O.S. and K.K.. K.A. supervised the greenhouse experiment and made initial quantitative analyses on the greenhouse data. T.A.K. collected the field data and participated in initial quantitative and selection analyses. J.H. examined the spatial autocorrelation. O.S. and K.K. conceived the original idea and supervised the overall project.

## Funding

This work was supported by Academy of Finland funds 309978 (for S.T.K.), 307582 (for O.S.), 307581 (for K.K.), and 351780 (for K.K.), and European Union funded projects Noveltree (for University of Oulu and Natural Resources Institute Finland, FP7 211868) and B4EST (for Natural Resources Institute Finland, H2020 773383).

*Conflict of interest:* The authors declare no conflicts of interest.

## Acknowledgments

We are thankful for Natural Resources Institute Finland staff in Loppi and Punkaharju for carrying out the greenhouse experiment and setting up and maintaining the field experiment. We acknowledge Satu Kumpula for helping with the field measurements, and Timo Knürr for advice with the quantitative data analyses.

## References

Abrahamsson, S., Nilsson, J., -E., Wu, H., García Gil, M. R., & Andersson, B. (2012). Inheritance of height growth and autumn cold hardness based on two generations of full-sib and half-sib families of *Pinus sylvestris*. *Scandinavian Journal of Forest Research*, 27(5), 405-413. https://doi.org/10.1080/02827581.2012.663403

Aho, M. -L. (1994). Autumn frost hardening of one-year-old *Pinus sylvestris* (L) seedlings: Effect of origin and parent trees. *Scandinavian Journal of Forest Research*, 9(1-4), 17-24. https://doi.org/10.1080/02827589409382808

Aitken, S. N., Yeaman, S., Holliday, J. A., Wang, T., & Curtis-McLane, S. (2008). Adaptation, migration or extirpation: Climate change outcomes for tree populations. *Evolutionary Applications*, 1(1), 95-111. https://doi.org/10.1111/j.1752-4571.2007.00013.x

Alberto, F. J., Aitken, S. N., Alía, R., González-Martínez, S. C., Hänninen, H., Kremer, A., Lefèvre, F., Lenormand, T., Yeaman, S., Whetten, R., & Savolainen, O. (2013). Potential for evolutionary responses to climate change—Evidence from tree populations. *Global Change Biology*, 19(6), 1645-1661. https://doi.org/10.1111/gcb.12181

Arnold, S. J., Bürger, R., Hohenlohe, P. A., Ajie, B. C., & Jones, A. G. (2008). Understanding the evolution and stability of the G-matrix. *Evolution*, 62(10), 2451-2461. https://doi.org/10.1111/j.1558-5646.2008.00472.x

Bärring, L., Berlin, M., & Andersson Gull, B. (2017). Tailored climate indices for climate-proofing operational forestry applications in Sweden and Finland. *International Journal of Climatology*, 37(1), 123-142. https://doi.org/10.1002/joc.4691

Barton, N. H. (1999). Clines in polygenic traits. *Genetical Research*, 74(3), 223-236. https://doi.org/10.1017/s001667239900422x

Benito Garzón, M., Alía, R., Robson, T. M., & Zavala, M. A. (2011). Intra-specific variability and plasticity influence potential tree species distributions under climate change. *Global Ecology and Biogeography*, 20(5), 766-778. https://doi.org/10.1111/j.1466-8238.2010.00646.x

Berlin, M., Persson, T., Jansson, G., Haapanen, M., Ruotsalainen, S., Bärring, L., & Andersson Gull, B. (2016). Scots pine transfer effect models for growth and survival in Sweden and Finland. *Silva Fennica*, 50(3), 1-21. https://doi.org/10.14214/sf.1562

Blows, M. W., & Brooks, R. (2003). Measuring nonlinear selection. *The American Naturalist*, 162(6), 815-820. https://doi.org/10.1086/378905

Bradshaw, W. E., & Holzapfel, C. M. (2001). Genetic shift in photo-periodic response correlated with global warming. *Proceedings of the National Academy of Sciences*, 98(25), 14509-14511. https://doi.org/10.1073/pnas.241391498

Brooks, R., Hunt, J., Blows, M. W., Smith, M. J., Bussiére, L. F., & Jennions, M. D. (2005). Experimental evidence for multivariate stabilizing sexual selection. *Evolution*, 59(4), 871-880. https://doi.org/10.1111/j.0014-3820.2005.tb01760.x

Caruso, C. M., Maherali, H., & Martin, R. A. (2020). A meta-analysis of natural selection on plant functional traits. *International Journal of Plant Sciences*, 181(1), 44-55. https://doi.org/10.1086/706199
suggested that shared variation between adaptive traits in *P. contorta* will result in correlated changes. In Scots pine with largely independent clines, adaptation can be flexible under changing conditions (e.g., de La Mata et al., 2017). Furthermore, high heritabilities in adaptive traits, combined with southern gene flow, will contribute to the adaptive capacity. Similar amounts of genetic variation likely exist in other parts of the clines, too, with a potential exception of range edges (Eckert et al., 2008; Kujala & Knürr et al., 2017).

This analysis concentrated on the capacity for genetic adaptation. To predict future responses, it will be necessary to also consider phenotypic plasticity (e.g., Benito Garzón et al., 2011; Hämälä et al., 2022), effects of demography, such as high survival of old trees limiting establishment possibilities (Kuparinen et al., 2010; Savolainen, 2004), and competition and other interactions with other species. Here our main findings come from a single population, and more empirical data should be obtained from other populations to verify the genetic independence of the seedling traits in other parts of the cline. However, this highly powered analysis of one large population already shows that correlated selection alone has not generated the parallel clines.

Current clines show that historically, genetic adaptation has been a very important component of responses to environmental changes (Davis & Shaw, 2001). Increasing information on the genetics of local adaptation in natural populations of forest trees will enable mitigation of climate change effects in artificially regenerated forests, too, via breeding for new conditions and determining use of regeneration material for changing climate (Aitken et al., 2008; Berlin et al., 2016; Cortés et al., 2020; MacLachlan et al., 2017).

## Supplementary material

Supplementary material is available online at *Evolution Letters*.

## Data and code availability

Data and R commands have been deposited to DRYAD (https://doi.org/10.5061/dryad.vmcvdnd09).

## Author contributions

S.T.K. was the main responsible author, aided by contributions from all authors. S.T.K. and H.K. performed the final quantitative and selection analyses. S.T.K. was responsible for interpretation of the results with O.S. and K.K.. K.A. supervised the greenhouse experiment and made initial quantitative analyses on the greenhouse data. T.A.K. collected the field data and participated in initial quantitative and selection analyses. J.H. examined the spatial autocorrelation. O.S. and K.K. conceived the original idea and supervised the overall project.

## Funding

This work was supported by Academy of Finland funds 309978 (for S.T.K.), 307582 (for O.S.), 307581 (for K.K.), and 351780 (for K.K.), and European Union funded projects Noveltree (for University of Oulu and Natural Resources Institute Finland, FP7 211868) and B4EST (for Natural Resources Institute Finland, H2020 773383).

*Conflict of interest:* The authors declare no conflicts of interest.

## Acknowledgments

We are thankful for Natural Resources Institute Finland staff in Loppi and Punkaharju for carrying out the greenhouse experiment and setting up and maintaining the field experiment. We acknowledge Satu Kumpula for helping with the field measurements, and Timo Knürr for advice with the quantitative data analyses.

## References

Abrahamsson, S., Nilsson, J., -E., Wu, H., García Gil, M. R., & Andersson, B. (2012). Inheritance of height growth and autumn cold hardness based on two generations of full-sib and half-sib families of *Pinus sylvestris*. *Scandinavian Journal of Forest Research*, 27(5), 405-413. https://doi.org/10.1080/02827581.2012.663403

Aho, M. -L. (1994). Autumn frost hardening of one-year-old *Pinus sylvestris* (L) seedlings: Effect of origin and parent trees. *Scandinavian Journal of Forest Research*, 9(1-4), 17-24. https://doi.org/10.1080/02827589409382808

Aitken, S. N., Yeaman, S., Holliday, J. A., Wang, T., & Curtis-McLane, S. (2008). Adaptation, migration or extirpation: Climate change outcomes for tree populations. *Evolutionary Applications*, 1(1), 95-111. https://doi.org/10.1111/j.1752-4571.2007.00013.x

Alberto, F. J., Aitken, S. N., Alía, R., González-Martínez, S. C., Hänninen, H., Kremer, A., Lefèvre, F., Lenormand, T., Yeaman, S., Whetten, R., & Savolainen, O. (2013). Potential for evolutionary responses to climate change—Evidence from tree populations. *Global Change Biology*, 19(6), 1645-1661. https://doi.org/10.1111/gcb.12181

Arnold, S. J., Bürger, R., Hohenlohe, P. A., Ajie, B. C., & Jones, A. G. (2008). Understanding the evolution and stability of the G-matrix. *Evolution*, 62(10), 2451-2461. https://doi.org/10.1111/j.1558-5646.2008.00472.x

Bärring, L., Berlin, M., & Andersson Gull, B. (2017). Tailored climate indices for climate-proofing operational forestry applications in Sweden and Finland. *International Journal of Climatology*, 37(1), 123-142. https://doi.org/10.1002/joc.4691

Barton, N. H. (1999). Clines in polygenic traits. *Genetical Research*, 74(3), 223-236. https://doi.org/10.1017/s001667239900422x

Benito Garzón, M., Alía, R., Robson, T. M., & Zavala, M. A. (2011). Intra-specific variability and plasticity influence potential tree species distributions under climate change. *Global Ecology and Biogeography*, 20(5), 766-778. https://doi.org/10.1111/j.1466-8238.2010.00646.x

Berlin, M., Persson, T., Jansson, G., Haapanen, M., Ruotsalainen, S., Bärring, L., & Andersson Gull, B. (2016). Scots pine transfer effect models for growth and survival in Sweden and Finland. *Silva Fennica*, 50(3), 1-21. https://doi.org/10.14214/sf.1562

Blows, M. W., & Brooks, R. (2003). Measuring nonlinear selection. *The American Naturalist*, 162(6), 815-820. https://doi.org/10.1086/378905

Bradshaw, W. E., & Holzapfel, C. M. (2001). Genetic shift in photo-periodic response correlated with global warming. *Proceedings of the National Academy of Sciences*, 98(25), 14509-14511. https://doi.org/10.1073/pnas.241391498

Brooks, R., Hunt, J., Blows, M. W., Smith, M. J., Bussiére, L. F., & Jennions, M. D. (2005). Experimental evidence for multivariate stabilizing sexual selection. *Evolution*, 59(4), 871-880. https://doi.org/10.1111/j.0014-3820.2005.tb01760.x

Caruso, C. M., Maherali, H., & Martin, R. A. (2020). A meta-analysis of natural selection on plant functional traits. *International Journal of Plant Sciences*, 181(1), 44-55. https://doi.org/10.1086/706199
Chen, J., Tsuda, Y., Stocks, M., Källman, T., Xu, N., Kärkkäinen, K., Huotari, T., Semerikov, V. L., Vendramin, G. G., & Lascoux, M. (2014). Clinal variation at phenology-related genes in spruce: Parallel evolution in FTL2 and Gigantea? *Genetics*, 197(3), 1025–1038. https://doi.org/10.1534/genetics.114.163063

Cortés, A. J., Restrepo-Montoya, M., & Bedoya-Canas, L. E. (2020). Modern strategies to assess and breed forest tree adaptation to changing climate. *Frontiers in Plant Science*, 11, 1606. https://doi.org/10.3389/fpls.2020.583323

Davis, M. B., & Shaw, R. G. (2001). Range shifts and adaptive responses to Quaternary climate change. *Science*, 292(5517), 673–679. https://doi.org/10.1126/science.292.5517.673

de La Mata, R., Hood, S., Sala, A., & Baldwin, I. T. (2017). Insect out-break shifts the direction of selection from fast to slow growth rates in the long-lived conifer *Pinus ponderosa*. *PNAS*, 114(28), 7391–7396. https://doi.org/10.2737/RDS-2017-0026

Dickerson, G. E. (1969). Techniques for research in quantitative animal genetics. In: Techniques and procedures in animal science research. *Am Soc Anim Sci*, Albany, N.Y. 36–79.

Eckert, C. G., Samis, K. E., & Lougheed, S. C. (2008). Genetic variation across species' geographical ranges: The central–marginal hypothesis and beyond. *Molecular Ecology*, 17(5), 1170–1188. https://doi.org/10.1111/j.1365-294X.2007.03659.x

Egbäck, S., Liziniewicz, M., Högberg, K. -A., Ekö, P. -M., & Nilsson, U. (2012). Influence of progeny and initial stand density on growth and quality traits of 21 year old half-sib Scots pine (Pinus sylvestris L). *Forest Ecology and Management*, 286, 1–7. https://doi.org/10.1016/j.foreco.2012.09.003

Eiche, V. (1966). Cold damage and plant mortality in experimental provenance plantations with Scots pine in northern Sweden. https://res.slu.se/id/publ/125307

Endler, J. A. (1977). Geographic variation, speciation, and clines (Issue 10). Princeton University Press.

Etterson, J. R., & Shaw, R. G. (2001). Constraint to adaptive evolution in response to global warming. *Science*, 294(5540), 151–154. https://doi.org/10.1126/science.1063656

EVOLTREE. (n.d.). EVOLTREE—Evolution of Trees as Drivers of Terrestrial Biodiversity. https://www.evoltree.eu/resources/intensive-study-sites/sites/site/punkaharju. Date accessed June 15, 2022.

Felsenstein, J. (1977). Multivariate normal genetic models with a finite number of loci. In *Proceedings of the International Conference on Quantitative Genetics* (pp. 227–246).

Fox, J., & Weisberg, S. (2019). An R companion to applied regression (3rd ed.). Sage. https://socialsciences.mcmaster.ca/jfox/Books/Companion/.

Geber, M. A., & Griffen, L. R. (2003). Inheritance and natural selection on functional traits. *International Journal of Plant Sciences*, 164(S3), S21–S42. https://doi.org/10.1086/368233

Geyer, C. J., Wagenius, S., & Shaw, R. G. (2007). Aster models for life history analysis. *Biometrika*, 94(2), 415–426. https://doi.org/10.1093/biomet/asm030

Grant, P.R., & Grant, B. R. (1995). Predicting microevolutionary responses to directional selection on heritable variation. *Evolution*, 49(2), 241–251. https://doi.org/10.1111/j.1558-5646.1995.tb02236.x

Guillaume, F. (2011). Migration-induced phenotypic divergence: The migration–selection balance of correlated traits. *Evolution*, 65(6), 1723–1738. https://doi.org/10.1111/j.1558-5646.2011.01248.x

Haig, I. T., Weidman, R. H., & Davis, K. P. (1941). *Natural regeneration in the western white pine type*. Technical Bulletin. USDA, Washington, D.C.

Hämälä, T., Ning, W., Kuittinen, H., Aryamanesh, N., & Savolainen, O. (2022). Environmental response in gene expression and DNA methylation reveals factors influencing the adaptive potential

of *Arabidopsis* *lyrata*. *ELife*, 11, e83115. https://doi.org/10.7554/ eLife.83115

Hänninen, H., Beuker, E., Johnsen, Ø., Leinonen, I., Murray, M., Sheppard, L., & Skrøppa, T. (2001). Impacts of climate change on cold hardiness of conifers. In *Conifer cold hardiness* (pp. 305–333). Springer. https://doi.org/10.1007/978-94-015-9650-3_12

Harju, A., & Muona, O. (1989). Background pollination in *Pinus sylvestris* seed orchards. *Scandinavian Journal of Forest Research*, 4(1-4), 513–520. https://doi.org/10.1080/02827588909382584

Howe, G. T., Aitken, S. N., Neale, D. B., Jermstad, K. D., Wheeler, N. C., & Chen, T. H. H. (2003). From genotype to phenotype: Unraveling the complexities of cold adaptation in forest trees. *Canadian Journal of Botany*, 81(12), 1247–1266. https://doi.org/10.1139/b03-141

Hurme, P., Repo, T., Savolainen, O., & Pääkkönen, T. (1997). Climatic adaptation of bud set and frost hardiness in Scots pine (*Pinus sylvestris*). *Canadian Journal of Forest Research*, 27(5), 716–723. https://doi.org/10.1139/x97-052

Hurme, P., Sillanpää, M. J., Arjas, E., Repo, T., & Savolainen, O. (2000). Genetic basis of climatic adaptation in Scots pine by Bayesian quantitative trait locus analysis. *Genetics*, 156(3), 1309–1322. https://doi.org/10.1093/genetics/156.3.1309

Jäghagen, K., & Albrektson, A. (1996). Induced competition among Scots pine seedlings and its effect on future timber quality. *New Forests*, 12(2), 163–174. https://doi.org/10.1007/bf00036627

Johansson, K., Hajek, J., Sjölin, O., & Normark, E. (2015). Early performance of Pinus sylvestris and *Picea abies*—A comparison between seedling size, species, and geographic location of the planting site. *Scandinavian Journal of Forest Research*, 30(5), 388–400. https://doi.org/10.1080/02827581.2014.987808

Kärkkäinen, K., & Savolainen, O. (1993). The degree of early inbreeding depression determines the selfing rate at the seed stage: Model and results from *Pinus sylvestris* (Scots pine). *Heredity*, 71(2), 160–166. https://doi.org/10.1038/hdy.1993.120

Kawlecki, T. J., & Ebert, D. (2004). Conceptual issues in local adaptation. *Ecology Letters*, 7(12), 1225–1241. https://doi.org/10.1111/j.1461-0248.2004.00684.x

Kelly, C. A. (1992). Spatial and temporal variation in selection on correlated life-history traits and plant size in *Chamaecrista fasciculata*. *Evolution*, 46(6), 1658–1673. https://doi.org/10.1111/j.1558-5646.1992.tb01160.x

Kingsolver, J. G., Hoekstra, H. E., Hoekstra, J. M., Berrigan, D., Vignieri, S. N., Hill, C. E., Hoang, A., Gibert, P., & Beerli, P. (2001). The strength of phenotypic selection in natural populations. *The American Naturalist*, 157(3), 245–261. https://doi.org/10.1086/319193

Kingsolver, J. G., & Pfennig, D. W. (2007). Patterns and power of phenotypic selection in nature. *Bioscience*, 57(7), 561–572. https://doi.org/10.1641/b570706

Koelewijn, H. P., Koski, V., & Savolainen, O. (1999). Magnitude and timing of inbreeding depression in Scots pine (*Pinus sylvestris* L). *Evolution*, 53(3), 758–768. https://doi.org/10.1111/j.1558-5646.1999.tb05370.x

Korkmaz, S., Göksülük, D., & Zararsiz, G. (2014). MVN: An R package for assessing multivariate normality. *R Journal*, 6(2), 151–162. https://doi.org/10.32614/RJ-2014-031

Koski, V., & Sievänen, R. (1985). Timing of growth cessation in relation to the variations in the growing season. In: Peter M. A. Tigerstedt, Pasi Puttonen, & Veikko Koski (Eds.), *Crop physiology of forest trees*. Helsinki University Press, 167–193.

Kujala, S. T., Knürr, T., Kärkkäinen, K., Neale, D. B., Sillanpää, M. J., & Savolainen, O. (2017). Genetic heterogeneity underlying variation in a locally adaptive clinal trait in *Pinus sylvestris* revealed by a Bayesian multipopulation analysis. *Heredity*, 118(5), 413–423. https://doi.org/10.1038/hdy.2016.115
Chen, J., Tsuda, Y., Stocks, M., Källman, T., Xu, N., Kärkkäinen, K., Huotari, T., Semerikov, V. L., Vendramin, G. G., & Lascoux, M. (2014). Clinal variation at phenology-related genes in spruce: Parallel evolution in FTL2 and Gigantea? *Genetics*, 197(3), 1025–1038. https://doi.org/10.1534/genetics.114.163063

Cortés, A. J., Restrepo-Montoya, M., & Bedoya-Canas, L. E. (2020). Modern strategies to assess and breed forest tree adaptation to changing climate. *Frontiers in Plant Science*, 11, 1606. https://doi.org/10.3389/fpls.2020.583323

Davis, M. B., & Shaw, R. G. (2001). Range shifts and adaptive responses to Quaternary climate change. *Science*, 292(5517), 673–679. https://doi.org/10.1126/science.292.5517.673

de La Mata, R., Hood, S., Sala, A., & Baldwin, I. T. (2017). Insect out-break shifts the direction of selection from fast to slow growth rates in the long-lived conifer *Pinus ponderosa*. *PNAS*, 114(28), 7391–7396. https://doi.org/10.2737/RDS-2017-0026

Dickerson, G. E. (1969). Techniques for research in quantitative animal genetics. In: Techniques and procedures in animal science research. *Am Soc Anim Sci*, Albany, N.Y. 36–79.

Eckert, C. G., Samis, K. E., & Lougheed, S. C. (2008). Genetic variation across species' geographical ranges: The central–marginal hypothesis and beyond. *Molecular Ecology*, 17(5), 1170–1188. https://doi.org/10.1111/j.1365-294X.2007.03659.x

Egbäck, S., Liziniewicz, M., Högberg, K. -A., Ekö, P. -M., & Nilsson, U. (2012). Influence of progeny and initial stand density on growth and quality traits of 21 year old half-sib Scots pine (Pinus sylvestris L). *Forest Ecology and Management*, 286, 1–7. https://doi.org/10.1016/j.foreco.2012.09.003

Eiche, V. (1966). Cold damage and plant mortality in experimental provenance plantations with Scots pine in northern Sweden. https://res.slu.se/id/publ/125307

Endler, J. A. (1977). Geographic variation, speciation, and clines (Issue 10). Princeton University Press.

Etterson, J. R., & Shaw, R. G. (2001). Constraint to adaptive evolution in response to global warming. *Science*, 294(5540), 151–154. https://doi.org/10.1126/science.1063656

EVOLTREE. (n.d.). EVOLTREE—Evolution of Trees as Drivers of Terrestrial Biodiversity. https://www.evoltree.eu/resources/intensive-study-sites/sites/site/punkaharju. Date accessed June 15, 2022.

Felsenstein, J. (1977). Multivariate normal genetic models with a finite number of loci. In *Proceedings of the International Conference on Quantitative Genetics* (pp. 227–246).

Fox, J., & Weisberg, S. (2019). An R companion to applied regression (3rd ed.). Sage. https://socialsciences.mcmaster.ca/jfox/Books/Companion/.

Geber, M. A., & Griffen, L. R. (2003). Inheritance and natural selection on functional traits. *International Journal of Plant Sciences*, 164(S3), S21–S42. https://doi.org/10.1086/368233

Geyer, C. J., Wagenius, S., & Shaw, R. G. (2007). Aster models for life history analysis. *Biometrika*, 94(2), 415–426. https://doi.org/10.1093/biomet/asm030

Grant, P.R., & Grant, B. R. (1995). Predicting microevolutionary responses to directional selection on heritable variation. *Evolution*, 49(2), 241–251. https://doi.org/10.1111/j.1558-5646.1995.tb02236.x

Guillaume, F. (2011). Migration-induced phenotypic divergence: The migration–selection balance of correlated traits. *Evolution*, 65(6), 1723–1738. https://doi.org/10.1111/j.1558-5646.2011.01248.x

Haig, I. T., Weidman, R. H., & Davis, K. P. (1941). *Natural regeneration in the western white pine type*. Technical Bulletin. USDA, Washington, D.C.

Hämälä, T., Ning, W., Kuittinen, H., Aryamanesh, N., & Savolainen, O. (2022). Environmental response in gene expression and DNA methylation reveals factors influencing the adaptive potential

of *Arabidopsis* *lyrata*. *ELife*, 11, e83115. https://doi.org/10.7554/efLife.83115

Hänninen, H., Beuker, E., Johnsen, Ø., Leinonen, I., Murray, M., Sheppard, L., & Skrøppa, T. (2001). Impacts of climate change on cold hardiness of conifers. In *Conifer cold hardiness* (pp. 305–333). Springer. https://doi.org/10.1007/978-94-015-9650-3_12

Harju, A., & Muona, O. (1989). Background pollination in *Pinus sylvestris* seed orchards. *Scandinavian Journal of Forest Research*, 4(1-4), 513–520. https://doi.org/10.1080/02827588909382584

Howe, G. T., Aitken, S. N., Neale, D. B., Jermstad, K. D., Wheeler, N. C., & Chen, T. H. H. (2003). From genotype to phenotype: Unraveling the complexities of cold adaptation in forest trees. *Canadian Journal of Botany*, 81(12), 1247–1266. https://doi.org/10.1139/b03-141

Hurme, P., Repo, T., Savolainen, O., & Pääkkönen, T. (1997). Climatic adaptation of bud set and frost hardiness in Scots pine (*Pinus sylvestris*). *Canadian Journal of Forest Research*, 27(5), 716–723. https://doi.org/10.1139/x97-052

Hurme, P., Sillanpää, M. J., Arjas, E., Repo, T., & Savolainen, O. (2000). Genetic basis of climatic adaptation in Scots pine by Bayesian quantitative trait locus analysis. *Genetics*, 156(3), 1309–1322. https://doi.org/10.1093/genetics/156.3.1309

Jäghagen, K., & Albrektson, A. (1996). Induced competition among Scots pine seedlings and its effect on future timber quality. *New Forests*, 12(2), 163–174. https://doi.org/10.1007/bf00036627

Johansson, K., Hajek, J., Sjölin, O., & Normark, E. (2015). Early performance of Pinus sylvestris and *Picea abies*—A comparison between seedling size, species, and geographic location of the planting site. *Scandinavian Journal of Forest Research*, 30(5), 388–400. https://doi.org/10.1080/02827581.2014.987808

Kärkkäinen, K., & Savolainen, O. (1993). The degree of early inbreeding depression determines the selfing rate at the seed stage: Model and results from *Pinus sylvestris* (Scots pine). *Heredity*, 71(2), 160–166. https://doi.org/10.1038/hdy.1993.120

Kawlecki, T. J., & Ebert, D. (2004). Conceptual issues in local adaptation. *Ecology Letters*, 7(12), 1225–1241. https://doi.org/10.1111/j.1461-0248.2004.00684.x

Kelly, C. A. (1992). Spatial and temporal variation in selection on correlated life-history traits and plant size in *Chamaecrista fasciculata*. *Evolution*, 46(6), 1658–1673. https://doi.org/10.1111/j.1558-5646.1992.tb01160.x

Kingsolver, J. G., Hoekstra, H. E., Hoekstra, J. M., Berrigan, D., Vignieri, S. N., Hill, C. E., Hoang, A., Gibert, P., & Beerli, P. (2001). The strength of phenotypic selection in natural populations. *The American Naturalist*, 157(3), 245–261. https://doi.org/10.1086/319193

Kingsolver, J. G., & Pfennig, D. W. (2007). Patterns and power of phenotypic selection in nature. *Bioscience*, 57(7), 561–572. https://doi.org/10.1641/b570706

Koelewijn, H. P., Koski, V., & Savolainen, O. (1999). Magnitude and timing of inbreeding depression in Scots pine (*Pinus sylvestris* L). *Evolution*, 53(3), 758–768. https://doi.org/10.1111/j.1558-5646.1999.tb05370.x

Korkmaz, S., Göksülük, D., & Zararsiz, G. (2014). MVN: An R package for assessing multivariate normality. *R Journal*, 6(2), 151–162. https://doi.org/10.32614/RJ-2014-031

Koski, V., & Sievänen, R. (1985). Timing of growth cessation in relation to the variations in the growing season. In: Peter M. A. Tigerstedt, Pasi Puttonen, & Veikko Koski (Eds.), *Crop physiology of forest trees*. Helsinki University Press, 167–193.

Kujala, S. T., Knürr, T., Kärkkäinen, K., Neale, D. B., Sillanpää, M. J., & Savolainen, O. (2017). Genetic heterogeneity underlying variation in a locally adaptive clinal trait in *Pinus sylvestris* revealed by a Bayesian multipopulation analysis. *Heredity*, 118(5), 413–423. https://doi.org/10.1038/hdy.2016.115
Kuparinen, A., Savolainen, O., & Schurr, F. M. (2010). Increased mortality can promote evolutionary adaptation of forest trees to climate change. *Forest Ecology and Management*, 259(5), 1003-1008.
https://doi.org/10.1016/j.foreco.2009.12.006

Lande, R. (1979). Quantitative genetic analysis of multivariate evolution, applied to brain: Body size allometry. *Evolution*, 33(1Part2), 402-416. https://doi.org/10.1111/j.1558-5646.1979.tb04694.x

Lande, R., & Arnold, S. J. (1983). The measurement of selection on correlated characters. *Evolution*, 37(6), 1210-1226. https://doi.org/10.1111/j.1558-5646.1983.tb00236.x

Le Corre, V., & Kremer, A. (2003). Genetic variability at neutral markers, quantitative trait loci and trait in a subdivided population under selection. *Genetics*, 164(3), 1205-1219. https://doi.org/10.1093/genetics/164.3.1205

Leck, M. A., Parker, V. T., Simpson, R. L., & Simpson, R. S. (2008). *Seedling ecology and evolution*. Cambridge University Press.

Linderholm, H. W., Walther, A., & Chen, D. (2008). Twentieth-century trends in the thermal growing season in the Greater Baltic Area. *Climatic Change*, 87(3-4), 405-419. https://doi.org/10.1007/s10584-007-9327-3

Liu, Y., & El-Kassaby, Y. A. (2019). Phenotypic plasticity of natural Populus trichocarpa populations in response to temporally environmental change in a common garden. BMC Evolutionary Biology, 19(1), 1-15. https://doi.org/10.1186/s12862-019-1553-6

Lönnroth, E. (1925). Untersuchungen über die innere Struktur und Entwicklung gleichaltriger naturnormaler Kiefernbestände: basiert auf Material aus der Südhälfte Finnlands. *Acta Forestalia Fennica*, 30, 269. http://hdl.handle.net/10138/17704

Luoranen, J., Saksa, T., & Lappi, J. (2018). Seedling, planting site and weather factors affecting the success of autumn plantings in Norway spruce and Scots pine seedlings. *Forest Ecology and Management*, 419-420, 79-90. https://doi.org/10.1016/j.foreco.2018.03.040

MacLachlan, I. R., McDonald, T. K., Lind, B. M., Rieseberg, L. H., Yeaman, S., & Aitken, S. N. (2021). Genome-wide shifts in climate-related variation underpin responses to selective breeding in a widespread conifer. *Proceedings of the National Academy of Sciences*, 118(10), e2016900118. https://doi.org/10.1073/pnas.2016900118

MacLachlan, I. R., Wang, T., Hamann, A., Smets, P., & Aitken, S. N. (2017). Selective breeding of lodgepole pine increases growth and maintains climatic adaptation. *Forest Ecology and Management*, 391, 404-416. https://doi.org/10.1016/j.foreco.2017.02.008

Merilä, J., & Björklund, M. (1999). Population divergence and morphometric integration in the greenfinch (Carduelis chloeris): Evolution against the trajectory of least resistance? *Journal of Evolutionary Biology*, 12(1), 103-112. https://doi.org/10.1046/j.1420-9101.1999.00013.x

Mikola, J. (1980). The effect of seed size and duration of growth on the height of Scots pine (Pinus sylvestris L.) and Norway spruce (Picea abies (L.) Karst.) provenances and progenies at the nursery stage. *Silva Fennica*, 14, 84-94.

Mikola, J. (1982). Bud-set phenology as an indicator of climatic adaptation of Scots pine in Finland [Pinus sylvestris]. *Silva Fennica*, 16, 178-184.

Mimura, M., & Aitken, S. N. (2007). Adaptive gradients and isolation-by-distance with postglacial migration in *Picea sitchensis*. *Heredity*, 99(2), 224-232. https://doi.org/10.1038/sj.hdy.6800987

Mirov, N. T. (1967). *The genus Pinus*. Ronald Press Company.

Muona, O., & Harju, A. (1989). Effective population sizes, genetic variability, and mating system in natural stands and seed orchards of Pinus sylvestris. *Silvae Genetica*, 38(5), 221-228.

Notivol, E., García-Gil, M. R., Alia, R., & Savolainen, O. (2007). Genetic variation of growth rhythm traits in the limits of a latitudinal cline in Scots pine. *Canadian Journal of Forest Research*, 37(3), 540-551. https://doi.org/10.1139/x06-243

Oleksyn, J., Tjoelker, M. G., & Reich, P. B. (1992). Growth and biomass partitioning of populations of European Pinus sylvestris L under simulated 50 and 60 N daylengths: Evidence for photo-periodic ecotypes. *New Phytologist*, 120(4), 561-574. https://doi.org/10.1111/j.1469-8137.1992.tb01806.x

Oleksyn, J., Tjoelker, M. G., & Reich, P. B. (1998). Adaptation to changing environment in Scots pine populations across a latitudinal gradient. *Silva Fennica*, 32(2), 129-140. https://doi.org/10.14214/sf.691

Olsen, J. E. (2010). Light and temperature sensing and signaling in induction of bud dormancy in woody plants. *Plant Molecular Biology*, 73(1-2), 37-47. https://doi.org/10.1007/s11103-010-9620-9

Phillips, P. C., & Arnold, S. J. (1989). Visualizing multivariate selection. *Evolution*, 43(6), 1209-1222. https://doi.org/10.1111/j.1558-5646.1989.tb02569.x

Pinheiro, J., Bates, D., DebRoy, S., & Sarkar, D.; R Core Team (2020). *nlme: Linear and nonlinear mixed effects models*. https://CRAN.R-project.org/package=nlme.

Pyhäjärvi, T., Kujala, S. T., & Savolainen, O. (2019). 275 Years of forestry meets genomics in Pinus sylvestris. *Evolutionary Applications*, 13(1), 11-30. https://doi.org/10.1111/eva.12809

R Core Team. (2022). R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. https://www.R-project.org/

Rausher, M. D. (1992). The measurement of selection on quantitative traits: Biases due to environmental covariances between traits and fitness. *Evolution*, 46(3), 616-626. https://doi.org/10.1111/j.1558-5646.1992.tb02070.x

Roff, D. A., & Fairbairn, D. J. (2012). A test of the hypothesis that correlational selection generates genetic correlations. *Evolution*, 66(9), 2953-2960. https://doi.org/10.1111/j.1558-5646.2012.01656.x

Saikkonen, K., Taulavuori, K., Hyvönen, T., Gundel, P. E., Hamilton, C. E., Vänninen, I., Nissinen, A., & Helander, M. (2012). Climate change-driven species' range shifts filtered by photoperiodism. *Nature Climate Change*, 2(4), 239-242. https://doi.org/10.1038/nclimate1430

Sarvas, R. (1962). Investigations on the flowering and seed crop of Pinus sylvestris. *Communicationes Instituti Forestalis Fenniae*, 53(4), 1-198.

Savolainen, O., Bokma, F., García-Gil, R., Komulainen, P., & Repo, T. (2004). Genetic variation in cessation of growth and frost hardiness and consequences for adaptation of Pinus sylvestris to climatic changes. *Forest Ecology and Management*, 197(1-3), 79-89. https://doi.org/10.1016/j.foreco.2004.05.006

Savolainen, O., Kujala, S. T., Sokol, C., Pyhäjärvi, T., Avia, K., Knürr, T., Kärkkäinen, K., & Hicks, S. (2011). Adaptive potential of northernmost tree populations to climate change, with emphasis on Scots pine (Pinus sylvestris L). *The Journal of Heredity*, 102(5), 526-536. https://doi.org/10.1093/jhered/esr056

Savolainen, O., Pyhäjärvi, T., & Knürr, T. (2007). Gene flow and local adaptation in trees. *Annual Review of Ecology, Evolution, and Systematics*, 38, 595-619. https://doi.org/10.1146/annurev.ecolsys.38.091206.095646

Siepielski, A. M., DiBattista, J. D., & Carlson, S. M. (2009). It's about time: The temporal dynamics of phenotypic selection in the wild. *Ecology Letters*, 12(11), 1261-1276. https://doi.org/10.1111/j.1461-0248.2009.01381.x
Kuparinen, A., Savolainen, O., & Schurr, F. M. (2010). Increased mortality can promote evolutionary adaptation of forest trees to climate change. *Forest Ecology and Management*, 259(5), 1003-1008.
https://doi.org/10.1016/j.foreco.2009.12.006

Lande, R. (1979). Quantitative genetic analysis of multivariate evolution, applied to brain: Body size allometry. *Evolution*, 33(1Part2), 402-416. https://doi.org/10.1111/j.1558-5646.1979.tb04694.x

Lande, R., & Arnold, S. J. (1983). The measurement of selection on correlated characters. *Evolution*, 37(6), 1210-1226. https://doi.org/10.1111/j.1558-5646.1983.tb00236.x

Le Corre, V., & Kremer, A. (2003). Genetic variability at neutral markers, quantitative trait loci and trait in a subdivided population under selection. *Genetics*, 164(3), 1205-1219. https://doi.org/10.1093/genetics/164.3.1205

Leck, M. A., Parker, V. T., Simpson, R. L., & Simpson, R. S. (2008). *Seedling ecology and evolution*. Cambridge University Press.

Linderholm, H. W., Walther, A., & Chen, D. (2008). Twentieth-century trends in the thermal growing season in the Greater Baltic Area. *Climatic Change*, 87(3-4), 405-419. https://doi.org/10.1007/s10584-007-9327-3

Liu, Y., & El-Kassaby, Y. A. (2019). Phenotypic plasticity of natural Populus trichocarpa populations in response to temporally environmental change in a common garden. BMC Evolutionary Biology, 19(1), 1-15. https://doi.org/10.1186/s12862-019-1553-6

Lönnroth, E. (1925). Untersuchungen über die innere Struktur und Entwicklung gleichaltriger naturnormaler Kiefernbestände: basiert auf Material aus der Südhälfte Finnlands. *Acta Forestalia Fennica*, 30, 269. http://hdl.handle.net/10138/17704

Luoranen, J., Saksa, T., & Lappi, J. (2018). Seedling, planting site and weather factors affecting the success of autumn plantings in Norway spruce and Scots pine seedlings. *Forest Ecology and Management*, 419-420, 79-90. https://doi.org/10.1016/j.foreco.2018.03.040

MacLachlan, I. R., McDonald, T. K., Lind, B. M., Rieseberg, L. H., Yeaman, S., & Aitken, S. N. (2021). Genome-wide shifts in climate-related variation underpin responses to selective breeding in a widespread conifer. *Proceedings of the National Academy of Sciences*, 118(10), e2016900118. https://doi.org/10.1073/pnas.2016900118

MacLachlan, I. R., Wang, T., Hamann, A., Smets, P., & Aitken, S. N. (2017). Selective breeding of lodgepole pine increases growth and maintains climatic adaptation. *Forest Ecology and Management*, 391, 404-416. https://doi.org/10.1016/j.foreco.2017.02.008

Merilä, J., & Björklund, M. (1999). Population divergence and morphometric integration in the greenfinch (Carduelis chloeris): Evolution against the trajectory of least resistance? *Journal of Evolutionary Biology*, 12(1), 103-112. https://doi.org/10.1046/j.1420-9101.1999.00013.x

Mikola, J. (1980). The effect of seed size and duration of growth on the height of Scots pine (Pinus sylvestris L.) and Norway spruce (Picea abies (L.) Karst.) provenances and progenies at the nursery stage. *Silva Fennica*, 14, 84-94.

Mikola, J. (1982). Bud-set phenology as an indicator of climatic adaptation of Scots pine in Finland [Pinus sylvestris]. *Silva Fennica*, 16, 178-184.

Mimura, M., & Aitken, S. N. (2007). Adaptive gradients and isolation-by-distance with postglacial migration in *Picea sitchensis*. *Heredity*, 99(2), 224-232. https://doi.org/10.1038/sj.hdy.6800987

Mirov, N. T. (1967). *The genus Pinus*. Ronald Press Company.

Muona, O., & Harju, A. (1989). Effective population sizes, genetic variability, and mating system in natural stands and seed orchards of Pinus sylvestris. *Silvae Genetica*, 38(5), 221-228.

Notivol, E., García-Gil, M. R., Alia, R., & Savolainen, O. (2007). Genetic variation of growth rhythm traits in the limits of a latitudinal cline in Scots pine. *Canadian Journal of Forest Research*, 37(3), 540-551. https://doi.org/10.1139/x06-243

Oleksyn, J., Tjoelker, M. G., & Reich, P. B. (1992). Growth and biomass partitioning of populations of European Pinus sylvestris L under simulated 50 and 60 N daylengths: Evidence for photo-periodic ecotypes. *New Phytologist*, 120(4), 561-574. https://doi.org/10.1111/j.1469-8137.1992.tb01806.x

Oleksyn, J., Tjoelker, M. G., & Reich, P. B. (1998). Adaptation to changing environment in Scots pine populations across a latitudinal gradient. *Silva Fennica*, 32(2), 129-140. https://doi.org/10.14214/sf.691

Olsen, J. E. (2010). Light and temperature sensing and signaling in induction of bud dormancy in woody plants. *Plant Molecular Biology*, 73(1-2), 37-47. https://doi.org/10.1007/s11103-010-9620-9

Phillips, P. C., & Arnold, S. J. (1989). Visualizing multivariate selection. *Evolution*, 43(6), 1209-1222. https://doi.org/10.1111/j.1558-5646.1989.tb02569.x

Pinheiro, J., Bates, D., DebRoy, S., & Sarkar, D.; R Core Team (2020). *nlme: Linear and nonlinear mixed effects models*. https://CRAN.R-project.org/package=nlme.

Pyhäjärvi, T., Kujala, S. T., & Savolainen, O. (2019). 275 Years of forestry meets genomics in Pinus sylvestris. *Evolutionary Applications*, 13(1), 11-30. https://doi.org/10.1111/eva.12809

R Core Team. (2022). R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. https://www.R-project.org/

Rausher, M. D. (1992). The measurement of selection on quantitative traits: Biases due to environmental covariances between traits and fitness. *Evolution*, 46(3), 616-626. https://doi.org/10.1111/j.1558-5646.1992.tb02070.x

Roff, D. A., & Fairbairn, D. J. (2012). A test of the hypothesis that correlational selection generates genetic correlations. *Evolution*, 66(9), 2953-2960. https://doi.org/10.1111/j.1558-5646.2012.01656.x

Saikkonen, K., Taulavuori, K., Hyvönen, T., Gundel, P. E., Hamilton, C. E., Vänninen, I., Nissinen, A., & Helander, M. (2012). Climate change-driven species' range shifts filtered by photoperiodism. *Nature Climate Change*, 2(4), 239-242. https://doi.org/10.1038/nclimate1430

Sarvas, R. (1962). Investigations on the flowering and seed crop of Pinus sylvestris. *Communicationes Instituti Forestalis Fenniae*, 53(4), 1-198.

Savolainen, O., Bokma, F., García-Gil, R., Komulainen, P., & Repo, T. (2004). Genetic variation in cessation of growth and frost hardiness and consequences for adaptation of Pinus sylvestris to climatic changes. *Forest Ecology and Management*, 197(1-3), 79-89. https://doi.org/10.1016/j.foreco.2004.05.006

Savolainen, O., Kujala, S. T., Sokol, C., Pyhäjärvi, T., Avia, K., Knürr, T., Kärkkäinen, K., & Hicks, S. (2011). Adaptive potential of northernmost tree populations to climate change, with emphasis on Scots pine (Pinus sylvestris L). *The Journal of Heredity*, 102(5), 526-536. https://doi.org/10.1093/jhered/esr056

Savolainen, O., Pyhäjärvi, T., & Knürr, T. (2007). Gene flow and local adaptation in trees. *Annual Review of Ecology, Evolution, and Systematics*, 38, 595-619. https://doi.org/10.1146/annurev.ecolsys.38.091206.095646

Siepielski, A. M., DiBattista, J. D., & Carlson, S. M. (2009). It's about time: The temporal dynamics of phenotypic selection in the wild. *Ecology Letters*, 12(11), 1261-1276. https://doi.org/10.1111/j.1461-0248.2009.01381.x
Simms, E. L. (1990). Examining selection on the multivariate phenotype: Plant resistance to herbivores. *Evolution*, 44(5), 1177-1188.
https://doi.org/10.1111/j.1558-5646.1990.tb05224.x

Singh, R. K., Svystun, T., AlDahmash, B., Jönsson, A. M., & Bhalerao, R. P. (2017). Photoperiod- and temperature-mediated control of phenology in trees—A molecular perspective. *New Phytologist*, 213(2), 511–524. https://doi.org/10.1111/nph.14346

Skrøppa, T., Solvin, T. M., & Steffenrem, A. (2023). Diallel crosses in *Picea abies* III variation and inheritance patterns in nursery trials. *Silvae Genetica*, 72, 49–57. https://doi.org/10.2478/sg-2023-0005

Slatkin, M. (1973). Gene flow and selection in a cline. *Genetics*, 75, 733–756. https://doi.org/10.1093/genetics/75.4.733

Slatkin, M. (1978). Spatial patterns in the distributions of polygenic characters. *Journal of Theoretical Biology*, 70, 213–228. https://doi.org/10.1016/0022-5193(78)90348-x

Squillace, A. E. (1974). Average genetic correlations among offspring from open-pollinated forest trees. *Silvae Genetica*, 23(5), 149–156.

Stinchcombe, J. R., Agrawal, A. F., Hohenlohe, P. A., Arnold, S. J., & Blows, M. W. (2008). Estimating nonlinear selection gradients using quadratic regression coefficients: Double or nothing? *Evolution*, 62(9), 2435–2440. https://doi.org/10.1111/j.1558-5646.2008.00449.x

Svensson, E. I., Arnold, S. J., Bürger, R., Csilléry, K., Draghi, J., Henshaw, J. M., Jones, A. G., de Lisle, S., Marques, D. A., McGuigan, K., Simon, M. N., & Runemark, A. (2021). Correlational selection in the age of genomics. *Nature Ecology & Evolution*, 5(5), 562–573. https://doi.org/10.1038/s41559-021-01413-3

Tyrmi, J. S., Vuosku, J., Acosta, J. J., Li, Z., Sterck, L., Cervera, M. T., Savolainen, O., & Pyhäjärvi, T. (2020). Genomics of clinal local adaptation in *Pinus sylvestris* under continuous environmental and spatial genetic setting. *G3 (Bethesda, Md.)*, 10(8), 2683–2696.
https://doi.org/10.1534/g3.120.401285

Vaartaja, O. (1959). Evidence of photoperiodic ecotypes in trees. *Ecological Monographs*, 29(2), 91–111. https://doi.org/10.2307/1942199

Viherä-Aarnio, A., Häkkinen, R., Partanen, J., Luomajoki, A., & Koski, V. (2005). Effects of seed origin and sowing time on timing of height growth cessation of *Betula pendula* seedlings. *Tree Physiology*, 25(1), 101–108. https://doi.org/10.1093/treephys/25.1.101

Walsh, B., & Lynch, M. (2018). *Evolution and selection of quantitative traits*. Oxford University Press.

Warwell, M., & Shaw, R. G. (2019). Phenotypic selection on ponderosa pine seed and seedling traits in the field under three experimentally manipulated drought treatments. *Evolutionary Applications*, 12(2), 159–174. https://doi.org/10.1111/eva.12685

Wheelwright, N. T., Keller, L. F., & Postma, E. (2014). The effect of trait type and strength of selection on heritability and evolvability in an island bird population. *Evolution*, 68(11), 3325–3336. https://doi.org/10.1111/evo.12499

Younginger, B. S., Sirová, D., Cruzan, M. B., & Ballhorn, D. J. (2017). Is biomass a reliable estimate of plant fitness? *Applications in Plant Sciences*, 5(2), 1600094. https://doi.org/10.3732/apps.1600094
Simms, E. L. (1990). Examining selection on the multivariate phenotype: Plant resistance to herbivores. *Evolution*, 44(5), 1177-1188.
https://doi.org/10.1111/j.1558-5646.1990.tb05224.x

Singh, R. K., Svystun, T., AlDahmash, B., Jönsson, A. M., & Bhalerao, R. P. (2017). Photoperiod- and temperature-mediated control of phenology in trees—A molecular perspective. *New Phytologist*, 213(2), 511–524. https://doi.org/10.1111/nph.14346

Skrøppa, T., Solvin, T. M., & Steffenrem, A. (2023). Diallel crosses in *Picea abies* III variation and inheritance patterns in nursery trials. *Silvae Genetica*, 72, 49–57. https://doi.org/10.2478/sg-2023-0005

Slatkin, M. (1973). Gene flow and selection in a cline. *Genetics*, 75, 733–756. https://doi.org/10.1093/genetics/75.4.733

Slatkin, M. (1978). Spatial patterns in the distributions of polygenic characters. *Journal of Theoretical Biology*, 70, 213–228. https://doi.org/10.1016/0022-5193(78)90348-x

Squillace, A. E. (1974). Average genetic correlations among offspring from open-pollinated forest trees. *Silvae Genetica*, 23(5), 149–156.

Stinchcombe, J. R., Agrawal, A. F., Hohenlohe, P. A., Arnold, S. J., & Blows, M. W. (2008). Estimating nonlinear selection gradients using quadratic regression coefficients: Double or nothing? *Evolution*, 62(9), 2435–2440. https://doi.org/10.1111/j.1558-5646.2008.00449.x

Svensson, E. I., Arnold, S. J., Bürger, R., Csilléry, K., Draghi, J., Henshaw, J. M., Jones, A. G., de Lisle, S., Marques, D. A., McGuigan, K., Simon, M. N., & Runemark, A. (2021). Correlational selection in the age of genomics. *Nature Ecology & Evolution*, 5(5), 562–573. https://doi.org/10.1038/s41559-021-01413-3

Tyrmi, J. S., Vuosku, J., Acosta, J. J., Li, Z., Sterck, L., Cervera, M. T., Savolainen, O., & Pyhäjärvi, T. (2020). Genomics of clinal local adaptation in *Pinus sylvestris* under continuous environmental and spatial genetic setting. *G3 (Bethesda, Md.)*, 10(8), 2683–2696.
https://doi.org/10.1534/g3.120.401285

Vaartaja, O. (1959). Evidence of photoperiodic ecotypes in trees. *Ecological Monographs*, 29(2), 91–111. https://doi.org/10.2307/1942199

Viherä-Aarnio, A., Häkkinen, R., Partanen, J., Luomajoki, A., & Koski, V. (2005). Effects of seed origin and sowing time on timing of height growth cessation of *Betula pendula* seedlings. *Tree Physiology*, 25(1), 101–108. https://doi.org/10.1093/treephys/25.1.101

Walsh, B., & Lynch, M. (2018). *Evolution and selection of quantitative traits*. Oxford University Press.

Warwell, M., & Shaw, R. G. (2019). Phenotypic selection on ponderosa pine seed and seedling traits in the field under three experimentally manipulated drought treatments. *Evolutionary Applications*, 12(2), 159–174. https://doi.org/10.1111/eva.12685

Wheelwright, N. T., Keller, L. F., & Postma, E. (2014). The effect of trait type and strength of selection on heritability and evolvability in an island bird population. *Evolution*, 68(11), 3325–3336. https://doi.org/10.1111/evo.12499

Younginger, B. S., Sirová, D., Cruzan, M. B., & Ballhorn, D. J. (2017). Is biomass a reliable estimate of plant fitness? *Applications in Plant Sciences*, 5(2), 1600094. https://doi.org/10.3732/apps.1600094