# Divergent selection predating the Last Glacial Maximum mainly acted on macro-phenotypes in Norway spruce

Mathieu Tiret, Lars Olsson, Thomas Grahn, Bo Karlsson, Pascal Milesi, Martin Lascoux, Sven-Olof Lundqvist, María Rosario García Gil

► To cite this version:

Mathieu Tiret, Lars Olsson, Thomas Grahn, Bo Karlsson, Pascal Milesi, et al. Divergent selection predating the Last Glacial Maximum mainly acted on macro-phenotypes in Norway spruce. Evolutionary Applications, 2022, 16 (1), pp.163-172. 10.1111/eva.13519. hal-03972770

HAL Id: hal-03972770

https://hal.science/hal-03972770

Submitted on 2 Jun 2023

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
# Divergent selection predating the Last Glacial Maximum mainly acted on macro-phenotypes in Norway spruce

Mathieu Tiret<fnref n="1" />,<fnref n="2" />,<fnref n="3" /> ID | Lars Olsson<fnref n="4" /> | Thomas Grahn<fnref n="4" /> | Bo Karlsson<fnref n="5" /> | Pascal Milesi<fnref n="1" /> ID |
Martin Lascoux<fnref n="1" /> ID | Sven-Olof Lundqvist<fnref n="6" /> | Maria Rosario García-Gil<fnref n="2" /> ID

**Correspondence**

Mathieu Tiret, Program in Plant Ecology and Evolution, Department of Ecology and Genetics, EBC and SciLife Lab, Uppsala University, Uppsala, Sweden.

Email: mathieu.tiret@inrae.fr

**Funding information**

B4Est; Sveriges Lantbruksuniversitet; Uppsala Universitet

**Abstract**

The current distribution and population structure of many species were, to a large extent, shaped by cycles of isolation in glacial refugia and subsequent population expansions. Isolation in and postglacial expansion through heterogeneous environments led to either neutral or adaptive divergence. Norway spruce is no exception, and its current distribution is the consequence of a constant interplay between evolutionary and demographic processes. We investigated population differentiation and adaptation of Norway spruce for juvenile growth, diameter of the stem, wood density, and tracheid traits at breast height. Data from 4461 phenotyped and genotyped Norway spruce from 396 half-sib families in two progeny tests were used to test for divergent selection in the framework of $Q_{ST}$ vs. $F_{ST}$. We show that the macroscopic resultant trait (stem diameter), unlike its microscopic components (tracheid dimensions) and juvenile growth, was under divergent selection that predated the Last Glacial Maximum. Altogether, the current variation in these phenotypic traits in Norway spruce is better explained by local adaptation to ancestral environments than to current ones, where populations were partly preadapted, mainly through growth-related traits.

**KEYWORDS**

Norway spruce, population structure, $Q_{ST}$ vs. $F_{ST}$, wood quality

## 1 | INTRODUCTION

During the Quaternary glacial cycles (up to 2.5 million years ago), many European animal and plant species were repeatedly confined to southern refugia. Those were often located in the Iberian, Italian, and Balkan peninsulas, but surviving populations could also be distributed on less discrete spots and, instead, scattered along the southern edge of the ice sheet (Bennett, 1997; Bhagwat et al., 2009; Kremer et al., 2010). Improvement of climatic conditions during post-glacial periods led to demographic expansions accompanied by population admixture (e.g., Petit et al., 2003), new selective pressures

(e.g., Suarez & Tsutsui, 2008), and inter- and/or intraspecific competition (e.g., Savolainen et al., 2007). Although the present population structure is the result of successive cycles of glaciation events that led to a series of demographic contractions and expansions, the current distribution of many species has, at least in part, been shaped by the different ecological processes and selective pressures that acted during the period after the Last Glacial Maximum (e.g., Palmé et al., 2003; Petit et al., 2002).

New abiotic and biotic pressures encountered during the glacial cycles presumably led to new adaptations. In many tree species, current phenotypes partly reflect past selection (Von Gadow et al., 2021),

<footnote n="1">¹Program in Plant Ecology and Evolution, Department of Ecology and Genetics, EBC and SciLife Lab, Uppsala University, Uppsala, Sweden</footnote>

<footnote n="2">²Department of Forest Genetics and Plant Physiology, SLU, Umeå Plant Science Centre (UPSC), Umeå, Sweden</footnote>

<footnote n="3">³IGEPP, INRAE, Institut Agro, Université de Rennes, Domaine de la Motte, Le Rheu, France</footnote>

<footnote n="4">⁴RISE Bioeconomy, Stockholm, Sweden</footnote>

<footnote n="5">⁵Skogforsk, Svalöv, Sweden</footnote>

<footnote n="6">⁶IIC, Stockholm, Sweden</footnote>

<footnote>This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.
© 2022 The Authors. Evolutionary Applications published by John Wiley & Sons Ltd.</footnote>
and selection is known to play an important role in shaping species distribution (Case et al., 2005; Li et al., 2022; Savolainen et al., 2007; Wisz et al., 2013). Among others, wood is one of the main tissues undergoing selection since radial growth of the woody stem provides the mechanical support for woody plants to grow tall and gain a competitive advantage (Falster & Westoby, 2003, 2005; Givnish, 1982; Iwasa et al., 1985). However, adaptation for upright growth is often associated with modified wood properties and evolutionary trade-offs (Chave et al., 2009; Poorter et al., 2006). Wood density is a good example of such a growth-survival trade-off (Kitajima, 1994; Poorter et al., 2010), where low wood density, resulting from wider lumina and, in hardwood species, wider vessels (e.g., Coomes et al., 2008; Enquist et al., 1999; Preston et al., 2006), is associated with fast growth. Conversely, high wood density provides higher biomechanical safety (e.g., Niklas, 1993; Van Gelder et al., 2006) and resistance against pathogens (Augspurger & Kelly, 1984), and is thus associated with higher survival rates. The growth-survival trade-off is a common ecological trade-off observed in both tropical (e.g., Falster & Westoby, 2005; Poorter, 2008) and boreal species (e.g., Pretzsch et al., 2018; Rossi et al., 2015).

Norway spruce (Picea abies (L.) H. Karst.) is a dominant boreal conifer species that thrives nowadays in a vast northern domain spanning from Norway to central Russia, as well as in smaller European domains, such as the Alps and the Carpathians (Li et al., 2022; Nota et al., 2022; Tsuda et al., 2016; Vendramin et al., 2000). Repeated demographic movements associated with climatic cycles and ensuing secondary contacts led to the division of Norway spruce into seven genetic clusters (Chen et al., 2019): Northern Fennoscandia (NFE), Central and South Sweden (CSE), Russia-Baltics (RBA), Northern Poland (NPL), Central Europe (CEU), Alpine (ALP), and Carpathian (ROM) domains. In Southern Sweden, recent hybrids between CSE and ALP clusters can be encountered (CSE-ALP). The timing of the divergence was given by that of the population structure that long predated the Last Glacial Maximum (10–15 million years ago; Chen et al., 2019; Li et al., 2022). The demographic history leading to this population structure went along a history of successive adaptations to biotic and abiotic factors. Following the hypothesis that most adaptations are reflected by cluster-specific allele frequencies (Kremer et al., 2012), the population structure of Norway spruce should partly be associated with signatures of past adaptations. Investigating this association is essential to understand how much of the current phenotypic variation is due to adaptations predating the Last Glacial Maximum.

The comparison of neutral ($F_{ST}$) and quantitative trait ($Q_{ST}$) genetic divergence offers a theoretically well-grounded, accessible, and powerful tool for assessing the role of natural selection in such a population differentiation scenario (Lamy et al., 2012). Under neutrality, it is expected that $F_{ST} = Q_{ST}$ (Whitlock, 2008), even if this framework needs to be updated to overcome some limitations (e.g., De Villemereuil et al., 2020). By implementing a robust way of estimating $F_{ST}$ and $Q_{ST}$, we used this framework to assess whether a phenotype diverges more between the genetic clusters than expected under random genetic drift. Additionally,

in order to account for possible evolutionary trade-offs between wood-related traits, we analyzed several phenotypes: stem diameter and ring widths, wood density, and tracheid traits defining it. Since the nature of selection between trees is expected to change over ontogenetic stages (Miriti, 2006), we studied the development of these phenotypes with cambial age/ring number from increment cores sampled from two progeny trials of the Swedish breeding program. We have conducted the analyses on a large dataset of 4461 20-year-old trees from 396 half-sib families, and 118,145 SNPs from exome sequencing. We show that the current phenotypic variation is partly explained by divergent selection that predated the Last Glacial Maximum, which mainly acted on macro-phenotypes (stem diameter) rather than micro-phenotypes (tracheid dimensions).

## 2 | MATERIALS AND METHODS

### 2.1 | Plant material

The study is based on plant material first presented by Chen et al. (2014), which is here analyzed with special attention to population structure and selection. The material originates from two progeny trials of Norway spruce at two different sites: FF1146 (Höreda, Sweden; 57°60′N 15°04′E) and FF1147 (Erikstorp, Sweden; 56°80′N 13°90′E). Both trials were established by the Forestry Research Institute of Sweden (Skogforsk) in the spring of 1990 from open-pollinated seeds harvested on selected trees from 112 stands in Southern Sweden. FF1146 comprised 1373 half-sib families and FF1147, 1375 half-sib families. In both trials, a randomized incomplete block design was implemented (20 blocks in FF1146 and 23 in FF1147), with single-tree plots and a homogeneous spacing of 1.4 m × 1.4 m. Dead trees and trees of low breeding interest were removed in 2008. After filtering out (i) half-sib families absent in the established Skogforsk clonal archives and (ii) stands with <3 half-sib families, the study included a total of 518 half-sib families with 5644 progeny (2983 from FF1146 and 2661 from FF1147).

### 2.2 | Population structure

Population structure was previously assessed at a fine scale on 1672 individuals using 399,801 unlinked putatively noncoding SNPs (see details in Chen et al., 2019), categorizing families into the following genetic clusters (ordered by decreasing latitude): Central and South Sweden (CSE), Russia-Baltics (RBA), Northern Poland (NPL), Central Europe (CEU), Alps (ALP), Carpathians (ROM), and hybrids between CSE and ALP (CSE-ALP) (for tree geographical origins, see figure 1 in Milesi et al., 2019). The families examined in this study followed the same genetic clustering (Figure 1a,b). Individuals from the ROM genetic cluster were removed, as our sample included only three such families, whereas all other clusters included at least 10 families (Table S1).
<figure><img src="image_2.png" /><figcaption>FIGURE 1 Overall description of phenotypes and genotypes of the sample of Norway spruce. (a) Principal component analysis of the genotypes of the half-sib families. Families are colored according to their genetic clusters (ordered by decreasing latitude in the legend): Central and South Sweden (CSE, dark red), Russia-Baltics (RBA, light green), Northern Poland (NPL, green), Central Europe (CEU, light blue), Alps (ALP, blue), and hybrids between CSE and ALP (CSE-ALP, pink). PC1 and PC2 explained 3.42% and 0.45% of the variance, respectively. (b) Admixture plot of the half-sib families, with K = 3. Individuals are grouped by population and sorted according to their coancestry coefficient. (c) Principal component analysis of the phenotypes: Diameter at breast height (DBH), ring width (RW), wood density (WD), tracheid radial width (TRW), tracheid tangential width (TTW), and tracheid Wall thickness (TWT). Left panel only shows variable axes and 60% confidence ellipses. On the right panel, each datum (dot) is colored according to its cambial age (3rd year in light blue up until 13th year in dark blue). PC1 and PC2 explained 41.2% and 38.3% of the variance, respectively.</figcaption></figure>

## 2.3 | Phenotype measurements

We reused the phenotypic dataset presented by Chen et al. (2014). Briefly, increment cores of 12 mm diameter, extracted from the northern side of the stem, were collected in 2010 and 2011 at breast height (BH, 1.3 m above ground) from six trees per half-sib family per trial. The usual cutting age is 60–70 years or more, while the cores were extracted from 20- or 21-year-old trees: only juvenile stages of wood characteristics were then captured by the cores. The increment cores were analyzed by SilviScan (Evans, 1994, 2006) at Innventia, now RISE (Stockholm, Sweden), assessing high-resolution pith-to-bark radial variation in growth, tracheid, and wood properties. In this study, we focused on the development with cambial age (ring number) of ring width (RW, in mm), and averages per ring of wood density (WD, in kg/m³), tracheid radial width (TRW, in µm), tracheid tangential width (TTW, in µm), tracheid wall thickness (TWT, in µm), diameter at breast height (DBH, in mm) calculated from ring widths, and the age when the trees reached breast height (ABH). Although ABH is related to upright growth, it can be estimated accurately with high-quality increment cores from old trees. The measurements, identification of annual rings, and the importance of ABH are discussed in Lundqvist et al. (2018).

The two innermost rings were removed as the extremely curved structure of the xylem does not allow the same high data quality as for rings further out. In addition, only data up to the 13th ring was kept for analyses since only fast-growing offspring with low ABH are represented at the highest cambial ages. To avoid the effect of the thinning performed, we removed all annual rings formed after 2008. The study then covered six traits (RW, WD, TRW, TTW, TWT, and DBH) for 11 rings per offspring in most of the cases and their ABH.

## 2.4 | Exome capture and SNP filtering

We reused the genotypic information of the plant material available on BioProject PRJNA511374 (Chen et al., 2019) and BioProject PRJNA731384 (Chen et al., 2021). The genotypic dataset is the exome capture (target resequencing) of the mother of each half-sib family, where 120-bp-length probes are designed to capture 26,219 genes of *Picea abies* (Vidalis et al., 2018). DNA extraction protocol of the 518 maternal trees is detailed in Chen et al. (2021), and SNP calling is detailed in Li et al. (2022).

As a complementary filter and as an alternative way to account for minor allele frequency (MAF), we additionally filtered out SNPs with less than five individuals per genotype to avoid a leverage effect. We performed a loose filtering for variant call rate (or missingness) with PLINK v1.9 (Chang et al., 2015), by removing individuals with more than 70% SNP missingness, resulting in 472 half-sib families, and by removing SNPs with more than 70% individual missingness, resulting in 118,145 SNPs (each mother is then genotyped for at least 35,443 SNPs). A total of 70,277 SNPs were successfully mapped on the consensus map of Bernhardsson et al. (2019).

## 2.5 | Statistical model

The traits (DBH, WD, RW, TRW, TTW, and TWT) were normalized by their respective overall means (across all trees and rings 3–13) to allow comparison between traits. Then, they were analyzed at each cambial age separately, resulting in 6 traits over 11 rings, a total of 66 traits in addition to ABH.

To study the family effect underlying the traits, we used the following mixed model:
$$Y_{ijkl} = \mu + s_i + B_{ij} + F_k + \epsilon_{ijkl}, \qquad (1)$$

where “$Y_{ijkl}$” is the trait of the l-th offspring of the k-th family in the j-th block at the i-th site, “$\mu$” is the grand mean, “$s_i$” is the fixed effect of the i-th site, “$B_{ij}$” and “$F_k$” are uncorrelated random effects of the j-th block nested in i-th site (with variance $\sigma_B^2$) and the k-th half-sib family (with variance $\sigma_F^2$), respectively, and “$\epsilon_{ijkl}$” is the residual fitted to a Gaussian distribution of variance $\sigma^2$. Heritability was estimated as $4 \times \sigma_F^2/(\sigma_F^2+\sigma_B^2+\sigma^2)$ as an intraclass correlation estimation of half-sib families (p.166, Falconer, 1989).

To investigate the effect of genetic clusters on traits jointly with the effect of ABH, we used the following linear mixed model:

$$Y_{ijkl} = \mu + c_l + s_i + g_k + (c \times g)_l + (c \times s)_l + B_{ij} + F_k + \epsilon_{ijkl}, \quad (2)$$

where “$c_l$” is the fixed effect of the l-th offspring’s ABH, “$(c \times g)_l$” is the fixed interaction effect between ABH and the genetic cluster of the l-th offspring, “$(c \times s)_l$” is the fixed interaction effect between ABH and the site of the l-th offspring, and the other terms are as defined above in model (1). Model (2) was applied to the 66 traits, as previously mentioned. Since the interaction between ABH and the site effect was always significant, we split the dataset per site, as recommended by Crawley (2007, pp. 323–386). The results are focusing on the Höreda trial, but the conclusions were the same in the Erikstorp trial.

Finally, to disentangle the effect of divergent selection and neutral divergence, we estimated $Q_{ST}$ as defined by Spitze (1993) with the following mixed model:

$$Y_{ijkl} = \mu + c_l + G_k + B_{ij} + F_k + \epsilon_{ijkl}, \qquad (3)$$

Where “$G_k$” is the random effect of the k-th family’s genetic cluster (with variance $\sigma_A^2$), and the other terms are as defined above in model (1). $Q_{ST}$ was estimated as $\sigma_A^2/(\sigma_A^2+2\sigma_F^2)$, and averaged across rings. Following the recommendation of Whitlock (2008), we compared $Q_{ST}$ to the positive part of the distribution of $F_{ST}$ (with VCFTools v0.1.16, option --weir-fst-pop; Danecek et al., 2011) computed on SNPs in chromosomal regions identified as noncoding by Nystedt et al. (2013). This set was not filtered for MAF or missingness, amounting then to 789,786 SNPs. If $Q_{ST}$ of a trait is larger than the 95 percentile of the empiric distribution of $F_{ST}$, the trait was considered to have diverged more than expected under neutrality in a heterogeneous environment. In contrast, if $Q_{ST}$ cannot be distinguished from $F_{ST}$ values, we

would have little evidence that selection is responsible for phenotypic divergence.

All computations were done using R v4.0.4 (R Core Team, 2021), and all models were fitted with the R package *lme4* (function *lmer*; Bates et al., 2015). Statistical significance of the fixed effects was assessed using a type II Wald chi-square test (function *Anova* of the R package *car*), and significance of a slope was assessed using a Student’s t-test (function *t.test* of the R package *stats*). To account for the unbalance caused by consecutive filtering, we only kept families with at least four offspring per site. Applying all the previous filters, we ended up with 396 half-sib families and 4461 offspring per trait. Multiple testing across the 66 traits were accounted for by considering a Bonferroni adjusted cut-off of $p = 7.6 \times 10^{-4}$.

## 3 | RESULTS

### 3.1 | The genetic cluster effect was significant at all ages

At all cambial ages, phenotypic variation of all six phenotypes (DBH, WD, RW, TRW, TTW, and TWT) was associated with the geographical origin of the genetic cluster, as the genetic cluster effect was significant for all rings ($\chi^2 > 11.6$, degree of freedom or df = 5, $p = 0.041$), except for the third ring TWT ($\chi^2 = 10.9$, df = 5, $p = 0.053$). Every trait of the southern populations was consistent with a higher investment in growth than in wood density (higher DBH, RW, TRW, and TTW; lower WD and TWT).

Phenotypic variation exhibited a strong temporal pattern, an ontogenetic shift (Figure 1c). For all rings, phenotypic variation also depended on the geographical origin of the trees. The ontogeny was quantitatively measurable, as it followed the second axis of the PCA on the six phenotypes (41.2% by PC1 and 38.3% by PC2). As expected, the ontogeny was mainly driven by that of DBH, which alone explained 34% of PC2, and on average increased fourfold over the cambial ages investigated, while other phenotypes varied at most by 54% (for further details, see Lundqvist et al., 2018). As cambial age increased, ontogenetic shift was instead driven by other phenotypes, as the coefficient of variation of DBH was the only one decreasing over time (Student’s $t = -6.13$, $p < 0.001$; Figure S1). Overall, the ontogenetic shift resulted from a complex interplay of multiple phenotypes, all with heritabilities higher than 30% (Table 1),

<table><caption>TABLE 1 Genetic variance ($\sigma_F^2$), environmental variance ($\sigma_B^2+\sigma^2$), and heritability ($h^2$) were estimated with model (1)</caption><thead><tr><th></th><th>ABH (year)</th><th>DBH (mm)</th><th>RW (mm)</th><th>WD (kg m<sup>-3</sup>)</th><th>TRW (µm)</th><th>TTW (µm)</th><th>TWT (µm)</th></tr></thead><tbody><tr><td>$\sigma_F^2$</td><td>7.44 &times; 10<sup>-2</sup></td><td>3.03 &times; 10<sup>-2</sup></td><td>5.13 &times; 10<sup>-3</sup></td><td>2.77 &times; 10<sup>-3</sup></td><td>5.99 &times; 10<sup>-4</sup></td><td>4.73 &times; 10<sup>-4</sup></td><td>2.47 &times; 10<sup>-3</sup></td></tr><tr><td>$\sigma_B^2+\sigma^2$</td><td>1.65</td><td>2.15 &times; 10<sup>-1</sup></td><td>5.63 &times; 10<sup>-2</sup></td><td>1.95 &times; 10<sup>-2</sup></td><td>6.04 &times; 10<sup>-3</sup></td><td>4.91 &times; 10<sup>-3</sup></td><td>1.88 &times; 10<sup>-2</sup></td></tr><tr><td>$h^2$ (%)</td><td>17.31</td><td>49.32</td><td>33.47</td><td>49.80</td><td>36.11</td><td>35.20</td><td>46.52</td></tr></tbody></table>

Traits: Calendar age at breast height (ABH), diameter at breast height (DBH), ring width (RW), wood density (WD), tracheid radial width (TRW), tracheid tangential width (TTW), and Tracheid Wall thickness (TWT).
in accordance with previously reported values (Chen et al., 2014, 2021; Hannrup et al., 2004).

## 3.2 | Accounting for early growth rate ruled out environmental causes from the analysis of ontogeny

By definition, ontogenic tempo is given by the ABH. It had a relatively higher environmental variance, and therefore a relatively lower heritability of 17.31% than other traits (Table 1). The lower heritability was not due to a lower genetic variance, in fact, quite the opposite. ABH also depended on the geographical origin of the genetic cluster as the genetic cluster effect was significant ($\chi^2 = 102.8$, df = 1, $p < 10^{-5}$ in model (1)). As mentioned earlier, ABH was latitudinally structured, half-sib families of southern origin having low ABH reflecting higher early growth rates.

Despite evidence supporting the genetic control of ABH, the genetic signals in the model (2) seem to have been absorbed by the genetic cluster and the family effects, and let the feature “ABH” only reflect environmental factors out of a lesser correlation with the response variable: when correcting for multiple testing, in model (2), the effects of ABH and the genetic clusters were significant for more than 70% of the phenotypes (47 of 66 for both), whereas their interaction was significant for <2% of the phenotypes (1 of 66; Table S2). In other words, ABH explained a large part of the phenotypic variation, while being orthogonal to the family or genetic cluster effect – the variable “ABH” was then likely only capturing the environmental part of ABH. In this particular statistical context, we will assume that ABH reflects environmental factors – else than the block effect already accounted for – impacting early growth and orthogonal to genetic factors (e.g., micro-environmental heterogeneity of the field, supposedly orthogonal to genetic factors because of the random block design). Under this assumption, it is possible to disentangle the effect of genetic factors on the phenotype of a “mature” tree from environmental factors giving advantages in early stages of life (i.e., acclimation). We will come back to this interpretation in the Discussion.

## 3.3 | Early favorable environment mainly impacted macro-phenotypes

The effect of acclimation, approximated by ABH, varied mainly with the cambial age, following a consistent pattern for every trait (but RW): acclimation mostly impacted early (3rd and 4th) or later (9th and higher) rings, and least impacted intermediate rings (5th to 8th). More precisely, analysis of variance of model (2) showed that, for early rings, ABH significantly impacted every trait but one trait. For intermediate rings, ABH significantly impacted only 6 of 20 traits. For later rings, ABH significantly impacted every trait but four traits (see details in Supplement method A). Only RW followed another pattern, for which ABH impacted mostly intermediate and later rings (5th ring and beyond) as the significance of ABH followed a

concave curve starting with an increase of significance from the 5th ring ($\chi^2 = 40.70$, df = 1, $p < 10^{-5}$) up until the 9th ring ($\chi^2 = 1199$, df = 1, $p < 10^{-5}$), then a decrease until the 13th ring ($\chi^2 = 54.50$, df = 1, $p < 10^{-5}$). This can also be seen in figure 5a of Lundqvist et al. (2018), where similar age curves, shifting sidewise with ABH, have their steepest phase during the intermediate cambial ages (see also Supplement method B).

On later rings, trees with a low ABH, that is, fast growers, ended up with a markedly larger DBH, lower WD and TWT, and slightly higher TRW and TTW. Fast growers had increasingly higher DBH over time, as the effect of ABH decreased from the 3rd ring to the 13th ring (from $5.36 \times 10^{-3}$ to $-2.87 \times 10^{-2}$). On the opposite, fast growers had increasingly lower WD, as the effect of ABH increased from the 3rd ring to the 13th ring (from $4.93 \times 10^{-3}$ to $3.1 \times 10^{-2}$). Fast growers also had significantly larger 9th rings, as the effect of ABH on RW followed a convex curve matching its concave significance, reaching its largest negative value at the 9th ring ($-5.82 \times 10^{-1}$). At all cambial ages, the effect of ABH was smaller for tracheid properties (TRW, TTW, and TWT) than for macro-phenotypes (DBH, WD) (Figure S2). These results show that environmental factors favoring early growth had an impact on all later-stage phenotypes (except for RW), albeit stronger on macro-phenotypes (DBH, WD) than on micro-phenotypes (TWT, TRW, and TTW).

## 3.4 | Macro-phenotypes were also the ones showing signature of past divergent selection

In addition to all phenotypes being, as shown above, significantly associated with geographical location, the significance of this association increased over cambial ages (Student's $t > 4.36$, $p < 0.002$), apart from TRW (Student's $t = -0.764$, $p = 0.464$). The significance of the genetic cluster in model (2) on the 13th ring was large for DBH and RW ($\chi^2 > 131.26$, df = 5, $p < 10^{-5}$), intermediate for WD and TWT ($\chi^2 > 79.51$, df = 5, $p < 10^{-5}$), and low for TRW and TTW ($\chi^2 < 19.39$, df = 5, $p < 0.002$). Presumably, the significance of TWT might be due to its high correlation with the significant WD ($p = 0.88$). The variation in growth-related macro-phenotypes, namely DBH and RW, therefore reflected the population structure more pronouncedly than did tracheid dimensions or wood density (Figure 2a).

Population structure was even more visible on phenotypes when correcting for environmental factors, as the PCA on family effects clearly exhibited the population structure (Figure 2b): genetic clusters were disposed along PC1 (58.9% of variance explained) in decreasing latitudinal order, except for CEU and ALP being reversed. PC1 was mostly driven by DBH, RW, TRW, and WD, while PC2 (17.0% of variance explained) was mostly driven by TTW and TWT. That the traits are structured by both family and genetic cluster can also be seen in the fact that intra-family (or intra-population) variance was significantly increased when the family (or population) structure was shuffled (see Supplement method C). Not only did phenotypic variation reflect the underlying population structure but it was also quantitatively associated with latitude.
<figure><img src="image_3.png" /></figure>

FIGURE 2 Genetic control of wood properties in Norway spruce. (a) ANOVA of genetic clusters in model (2) as a function of cambial age (from 3 to 13), and of traits: DBH (plain circle), RW (plain triangle), WD (plain square), TRW (circle), TTW (triangle), and TWT (square). (b) Principal component analysis of the residuals of the phenotypes according to model (2): Diameter at breast height (DBH), ring width (RW), wood density (WD), tracheid radial width (TRW), tracheid tangential width (TTW), and tracheid Wall thickness (TWT). Left panel only shows variable axes and 10% confidence ellipses (zoomed ×2.5), where clusters were ordered in decreasing latitude from right to left. On the right panel, each data (dot) are colored according to its genetic cluster (ordered by decreasing latitude in the legend): Central and South Sweden (CSE, dark red), Russia-Baltics (RBA, light green), Northern Poland (NPL, green), Central Europe (CEU, light blue), Alps (ALP, blue), and hybrids between CSE and ALP (CSE-ALP, pink). PC1 and PC2 explained 57.5% and 18.2% of the variance, respectively. (c) Distribution of pairwise $F_{ST}$ across the genome, averaged among pairs of populations. The vertical solid black line is the average $F_{ST}$ over the genome (0.0488), and vertical dotted lines are the $Q_{ST}$ estimated for different traits (averaged across rings): Diameter at breast height (DBH, red, $Q_{ST} = 0.281$), ring width (RW, green, $Q_{ST} = 0.302$), wood density (WD, brown, $Q_{ST} = 0.121$), tracheid radial width (TRW, turquoise, $Q_{ST} = 0.0579$), tracheid tangential width (TTW, blue, $Q_{ST} = 0.177$), Tracheid Wall thickness (TWT, purple, $Q_{ST} = 0.112$), and age at breast height (ABH, pink, $Q_{ST} = 0.066$).

The quantitative association between environmental gradients and phenotypic variation was partly due to selection (Figure 2c), at least for macro-phenotypes. Selection, as assessed by $Q_{ST}$ from model (3), was only significant for DBH ($p = 0.0115$), RW ($p = 0.0086$), and marginally so for TTW ($p = 0.0437$). Despite the $Q_{ST}$ of TTW being marginally significant, what is striking is that TRW, which is a component of DBH and RW, is not ($p = 0.2680$). Comparing $Q_{ST}$ and $F_{ST}$ demonstrate that selection, if any, would primarily act upon growth-related macroscopic traits. This strong association between population structure and signatures of selection made impossible the inference of genetic variants underlying phenotypic variation (see details in Supplement method D; Figures S3-S4). Overall, phenotypic divergence in growth-related macroscopic traits (DBH and RW) is partly due to divergent selection on top of neutral divergence.

4 | DISCUSSION

The Swedish breeding program is a mixture of local and recently introduced trees from other European countries. As such, although incomplete, the progeny trials reflect the population structure due

to glaciation events and constitute an important source of information for evolutionary biology (Chen et al., 2019; Li et al., 2022; Milesi et al., 2019). The breeding conditions constituted a unique opportunity to accurately correct phenotypes for environmental factors, by ruling out recruitment selection and considerably lowering interspecific competition. Furthermore, by correcting for population structure and using common garden experiments, we can rule out neutral and plastic processes so that the observed phenotypic differences can be interpreted as being the outcome of adaptation (Savolainen et al., 2013; de Villemereuil et al. et al., 2020).

In our common garden experiments, growth exhibited a high level of population differentiation, as expected (Kremer et al., 2012; Savolainen et al., 2007). Population differentiation followed a latitudinal pattern, as previously reported (Milesi et al., 2019; Montwé et al., 2016; Rossi et al., 2015). The growth rate was generally differentiated among individual trees at early stages, supporting the fact that young ages are the most dynamic phases (Lundqvist et al., 2018). Differences in apparent growth rate were mainly due to a longer period of growth rather than a faster daily growth (see discussion in Milesi et al., 2019), and despite the low heritability of growth rate, its high genetic and environmental variances suggest that it could be related to fitness (Houle, 1992; Caballero, 2020, p.235). This is
also supported by the fact that the effect of ABH (here, only capturing the micro-environmental effect) was strongest on phenotypes on which signature of selection was also the strongest. Early advantages, possibly acquired by chance, seem to be re-invested in macro-phenotypes, giving thus competitive advantages.

## 4.1 | The quantitative association between phenotypes and population structure predates the last glacial maximum

Both genetic and phenotypic variation had a clear latitudinal component, as previously reported in Norway spruce (e.g., Milesi et al., 2019) and in other species (e.g., Csilléry et al., 2020; Sang et al., 2019). Since local adaptation matches the environment (Savolainen et al., 2013) and postglacial recolonization routes followed environmental gradients (Gaggiotti et al., 2009; Li et al., 2022), population structure of Norway spruce is, as a result, strongly confounded with local adaptation. Population structure is thus a good predictor of selection in Norway spruce, as previously reported (Collignon et al., 2002; Kremer et al., 2010; Savolainen et al., 2007). However, local adaptation is not necessarily limited to recent population history. Some of the adaptations displayed by Norway spruce were associated with ancestral environments and some with current environments. If we consider that adaptation can be dated by the population structure it is associated with, knowing that the population structure of Norway spruce largely predates the Last Glacial Maximum (10–15 million years ago; Chen et al., 2019; Li et al., 2022), the observed phenotypic divergence is partly due to adaptations also predating the Last Glacial Maximum. We posit that the quantitative association among phenotype, latitude, and population structure would not be observed if there were no local adaptation to ancestral environments. The last postglacial expansion of Norway spruce was thus mainly shaped by adaptations *formerly acquired* – *not newly acquired*. In other words, the phenotypic variation within the refugia explains in part the latitudinal variation observed nowadays at large geographical scale.

Because of the association between phenotypic variation and demographic processes, assessing the genetic control underlying adaptation of Norway spruce is challenging: not correcting for population structure in GWASs will lead to a high false-positive rate (Barton et al., 2019; Lawson et al., 2019), and correcting for population structure may remove true association, that is, lead to a high false-negative rate. A comprehensive assessment of the genetic variants underlying phenotypic variation is then not possible when there is such a strong latitudinal component. Increasing the number of samples or the number of SNPs remains inconclusive (Chen et al., 2021). A possible workaround suggested by Milesi et al. (2019) is to consider phenotypes that are less latitudinally associated, such as the longitudinally associated bud burst, and correct other phenotypes by the association with population structure thereby inferred. Alternatively, we can consider several species jointly as suggested by Alberto et al. (2013).

## 4.2 | Divergent selection mostly affected resultant traits, not components

The strongest genetic divergence was observed for growth-related traits, and secondly for wood density, suggesting that most of the adaptation is associated with radial or upright growth. Divergent selection acting at early life stages most likely led to the accumulation of competitive advantages, causing then an observable ontogenetic shift of the effect of population structure. The average value of $F_{ST}$ was slightly lower than previously reported (0.0488 in our study, and ~0.116–0.145 in table 1 of Chen et al., 2019). However, overall, and even accounting for the challenging aspect of interpreting $Q_{ST}$ vs. $F_{ST}$ comparisons (e.g., Edelaar & Björklund, 2011; Whitlock, 2008), correcting the estimates of $F_{ST}$ would not change the conclusion of our study: if there is selection, it would affect growth-related macroscopic traits more strongly than wood density or its microscopic components. Consistently, the PCA on family effects showed that the latitudinal component – and so, population structure – was mostly explained by stem diameter, ring width, and wood density. Neutral divergence could not be excluded for early growth rate (ABH) in spite of its clear relation with growth. Either this trait is neutral, which is highly unlikely as early-stage selection is known to be strong (e.g., Saleh et al., 2022), or, more likely, selection on ABH is invariant among genetic clusters. In the latter case, selection of early growth rate would be strong regardless of the heterogeneity of the locations, where variability would not come from divergence but from, for instance, trade-offs with other life-history traits (e.g., wood density).

The difference between $Q_{ST}$ and $F_{ST}$ is due to the covariance between effect sizes, and the covariance between allelic frequencies (reviewed in Le Corre & Kremer, 2012). Altogether, the main pattern of our analyses is that adaptation is more tightly associated with traits that are composite or resultant traits rather than with their components themselves, in agreement with Caballero (2020, p. 253). Indeed, components can have random variations that are invisible to selection – therefore neutral, as long as the components contribute to a resultant trait that is fit. In other words, components have a higher degree of freedom, and therefore are noisier. Since genetic polymorphism depicted by SNPs is partly linked to selection, genomic selection on resultant traits – such as stem diameter – should be much more efficient than on components – such as treasury traits. An alternative explanation is simply that components were not associated with divergent selection predating the Last Glacial Maximum, and might be under recent selection – not detectable with our approach.

## 5 | CONCLUSION

As expected, phenotypic variation reflected the population structure of Norway spruce. The fact that the phenotypes were latitudinally ordered suggests that the origin of population was strongly confounded with adaptation. Teasing apart the demographic history and
the adaptation history is challenging, more so that it has been demonstrated that population history had a strong impact on phenotypic divergence in Norway spruce (Lagercrantz & Ryman, 1990; Milesi et al., 2019). Divergence of phenotypic traits might in part be the results of past adaptations, which can be dated by the population structure it is associated with. In accordance with previous studies (Martínez-Sancho et al., 2021; Milesi et al., 2019; Savolainen et al., 2013), our results support the hypothesis that the current distribution of Norway spruce reflects in part adaptations to ancestral environments, as some signatures of selection were associated with population structure that predated the Last Glacial Maximum. The divergent selection affected more growth-related macro-phenotypes than micro-phenotypes at tracheid level – presumably invisible to selection.

## ACKNOWLEDGMENTS

The authors acknowledge the Swedish research program Bio4Energy and Skogforsk for access to phenotypic data, and the Swedish Foundation for Strategic Research (SSF), grant number RBP14-0040 for exome capture data. This study was funded by the following sources: Uppsala University, SLU Umeå, and the European Union's Horizon 2020 B4EST for basic functioning and postdoctoral grant for MT. The authors would also like to thank Sophie Karrenberg and Zhi-Qiang Chen for their useful feedback. The computations and data handling were enabled by resources provided by the Swedish National Infrastructure for Computing (SNIC 2017-7-296) at UppMax partially funded by the Swedish Research Council through grant agreement no. 2018-05973.

## CONFLICT OF INTEREST

The authors declare that there is no conflict of interest.

## DATA AVAILABILITY STATEMENT

The exome capture raw reads and the RNA-seq data have been deposited in NCBI's sequence read archive (SRA) under accession number PRJNA511374 (Chen et al., 2019) and PRJNA731384 (Chen et al., 2021). For question on availability of the raw phenotypic data, contact RGG.

## ORCID

Mathieu Tiret ID https://orcid.org/0000-0001-5389-4886

Pascal Milesi ID https://orcid.org/0000-0001-8580-4291

Martin Lascoux ID https://orcid.org/0000-0003-1699-9042

Maria Rosario García-Gil ID https://orcid.org/0000-0002-6834-6708

## REFERENCES

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.

Augspurger, C. K., & Kelly, C. K. (1984). Pathogen mortality of tropical tree seedlings: Experimental studies of the effects of dispersal distance, seedling density, and light conditions. *Oecologica, 61*, 211–217.

Barton, N., Hermisson, J., & Nordborg, M. (2019). Population genetics: Why structure matters. *eLife*, 8, e45380.

Bates, D., Mächler, M., Bolker, B., & Walker, S. (2015). Fitting linear mixed-effects models using lme4. *Journal of Statistical Software*, 67(1), 1–48.

Bennett, K. D. (1997). *Evolution and ecology: The pace of life*. Cambridge University Press.

Bernhardsson, C., Vidalis, A., Wang, X., Scofield, D. G., Schiffthaler, B., Baison, J., Street, N. R., García-Gil, M. R., & Ingarsson, P. K. (2019). An ultra-dense haploid genetic map for evaluating the highly fragmented genome assembly of Norway spruce (*Picea abies*). *G3 Genes* | *Genomes* | *Genetics*, 9, 1623–1632.

Binney, H. A., Willis, K. J., Edwards, M. E., Bhagwat, S. A., Anderson, P. M., Andreev, A. A., Blaauw, M., Damblon, F., Haesaerts, P., Kienast, F., Kremenetski, K. V., Krivonogov, S. K., Lozhkin, A. V., MacDonald, G. M., Novenko, E. Y., Oksanen, P., Sapelko, T. V., Väliranta, M., & Vazhenina, L. (2009). The distribution of late-quaternary woody taxa in northern Eurasia: Evidence from a new macrofossil database. *Quaternary Science Reviews*, 28, 2445–2464.

Caballero, A. (2020). *Quantitative genetics* (p. 235). Cambridge University Press.

Case, T. J., Holt, R. D., McPeek, M. A., & Keitt, T. H. (2005). The community context of species' borders: Ecological and evolutionary perspectives. *Oikos*, 108, 28–46.

Chang, C. C., Chow, C. C., Tellier, L. C. A. M., Vattikuti, S., Purcell, S. M., & Lee, J. J. (2015). Second-generation PLINK: Rising to the challenge of larger and richer datasets. *GigaScience*, 4(1), 7.

Chave, J., Coomes, D., Jansen, S., Lewis, S. L., Swenson, N. G., & Zanne, A. E. (2009). Towards a worldwide wood economics spectrum. *Ecology Letters*, 12(4), 351–366.

Chen, J., Li, L., Milesi, P., Jansson, G., Berlin, M., Karlsson, B., Aleksic, J., Vendramin, G. G., & Lascoux, M. (2019). Genomic data provide new insights on the demographic history and the extent of recent material transfers in Norway spruce. *Evolutionary Applications*, 12, 1539–1551.

Chen, Z. Q., García-Gil, M. R., Karlsson, B., Lundqvist, S.-O., Olsson, L., & Wu, H. X. (2014). Inheritance of growth and solid wood quality traits in a large Norway spruce population tested at two locations in southern Sweden. *Tree Genetics & Genomes*, 10(5), 1291–1303.

Chen, Z.-Q., Zan, Y., Milesi, P., Zhou, L., Chen, J., Li, L., Cui, B. B., Niu, S., Westin, J., Karlsson, B., García-Gil, M. R., Lascoux, M., & Wu, H. X. (2021). Leveraging breeding programs and genomic data in Norway spruce (*Picea abies* L. karst) for GWAS analysis. *Genome Biology*, 22, 179.

Collignon, A. M., Sype, H. V. D., & Favre, J. M. (2002). Geographical variation in random amplified polymorphic DNA and quantitative traits in Norway spruce. *Canadian Journal of Forest Research*, 32(2), 266–282.

Coomes, D. A., Heathcote, S., Godfrey, E. R., Shepherd, J. J., & Sack, L. (2008). Scaling of xylem vessels and veins within the leaves of oak species. *Biology Letters*, 4(3), 302–306.

Crawley, M. J. (2007). Statistical modelling. In *The R book* (pp. 323–386). John Wiley & Sons, Ltd.

Csilléry, K., Ovaskainen, O., Sperisen, C., Buchmann, N., Widmer, A., & Gugerli, F. (2020). Adaptation to local climate in multi-trait space: Evidence from silver fir (*Abies alba* mill.) populations across a heterogeneous environment. *Heredity*, 124(1), 77–92.

Danecek, P., Auton, A., Abecasis, G., Albers, C. A., Banks, E., DePristo, M., Handsaker, R. E., Lunter, G., Marth, G. T., Sherry, S. T., McVean, G., Durbin, R., & 1000 Genomes Project Analysis Group. (2011). The variant call format and VCFtools. *Bioinformatics*, 27(15), 2156–2158.

De Villemereuil, P., Gaggiotti, O. E., & Gouder, J. (2020). Common garden experiments to study local adaptation need to account for population structure. *Journal of Ecology*, 100, 1–5.

Edelaar, P. I. M., & Björklund, M. (2011). If FST does not measure neutral genetic differentiation, then comparing it with QST is misleading. Or is it? *Molecular Ecology*, 20(9), 1805–1812.
Enquist, B. J., West, G. B., Charnov, E. L., & Brown, J. H. (1999). Allometric scaling of production and life-history variation in vascular plants. *Nature*, **401**, 907-911.

Evans, R. (1994). Rapid measurement of the transverse dimensions of tracheids in radial wood sections from Pinus Radiata. *Holzforschung*, **48**, 168-172.

Evans, R. (2006). Wood stiffness by X-ray diffractometry. In D. D. Stokke & H. L. Groom (Eds.), *Characterization of the cellulosic cell wall* (pp. 138-146). Wiley.

Falconer, D. S. (1989). *Introduction to quantitative genetics* (3rd ed., pp. 116). Pearson Education India.

Falster, D. S., & Westoby, M. (2003). Plant height and evolutionary games. *Trends in Ecology & Evolution*, **18**(7), 337-343.

Falster, D. S., & Westoby, M. (2005). Alternative height strategies among 45 dicot rain forest species from tropical Queensland. *Australian Journal of Ecology*, **93**(3), 521-535.

Gaggiotti, O. E., Bekkevold, D., Jørgensen, H. B. H., Foll, M., Carvalho, G. R., Andre, C., & Ruzzante, D. E. (2009). Disentangling the effects of evolutionary, demographic, and environmental factors influencing genetic structure of natural populations: Atlantic herring as a case study. *Evolution*, **63**, 2939-2951.

Givnish, T. J. (1982). On the adaptive significance of leaf height in forest herbs. *The American Naturalist*, **120**(3), 353-381.

Hannrup, B., Cahalan, C., Chantre, G., Grabner, M., Karlsson, B., Bayon, I. L., Jones, G. L., Müller, U., Pereira, H., Rodrigues, J. C., Rosner, S., Wilhelmsson, L., & Wimmer, R. (2004). Genetic parameters of growth and wood quality traits in *Picea abies*. *Scandinavian Journal of Forest Research*, **19**(1), 14-29.

Houle, D. (1992). Comparing evolvability and variability of quantitative traits. *Genetics*, **130**(1), 195-204.

Iwasa, Y., Cohen, D., & Léon, J. A. (1985). Tree height and crown shape, as results of competitive games. *Journal of Theoretical Biology*, **112**(2), 279-297.

Kitajima, K. (1994). Relative importance of photosynthetic traits and allocation patterns as correlates of seedling shade tolerance of 13 tropical trees. *Oecologia*, **98**, 419-428.

Kremer, A., Le Corre, V., Petit, R. J., & Ducousso, A. (2010). Historical and contemporary dynamics of adaptive differentiation in European oaks. In *Molecular approaches in natural resource conservation* (pp. 101-122). Cambridge University Press.

Kremer, A., Ronce, O., Robledo-Arnuncio, J. J., Guillaume, F., Bohrer, G., Nathan, R., Bridle, J. R., Gomulkiewicz, R., Klein, E. K., Ritland, K., Kuparinen, A., Gerber, S., & Schueler, S. (2012). Long-distance gene flow and adaptation of forest trees to rapid climate change. *Ecology Letters*, **15**(4), 378-392.

Lagercrantz, U., & Ryman, N. (1990). Genetic structure of Norway spruce (*Picea abies*): Concordance of morphological and allozymic variation. *Evolution*, **44**(1), 38-53.

Lamy, J. B., Plomion, C., Kremer, A., & Delzon, S. (2012). QST<FST as a signature of canalization. *Molecular Ecology*, **21**(23), 5646-5655.

Lawson, D. J., Davies, N. M., Haworth, S., Ashraf, B., Howe, L., Crawford, A., Hemani, G., Smith, G. D., & Timpson, N. J. (2019). Is population structure in the genetic biobank era irrelevant, a challenge, or an opportunity? *Human Genetics*, **139**, 23-41.

Le Corre, V., & Kremer, A. (2012). The genetic differentiation at quantitative trait loci under local adaptation. *Molecular Ecology*, **21**(7), 1548-1566.

Li, L., Milesi, P., Tiret, M., Chen, J., Sendrowski, J., Baison, J., Chen, Z. Q., Zhou, L., Karlsson, B., Berlin, M., Westin, J., Garcia-Gil, M. R., Wu, H. X., & Lascoux, M. (2022). Teasing apart the joint effect of demography and natural selection in the birth of a contact zone. *New Phytologist*, **236**, 1976-1987.

Lundqvist, S.-O., Seifert, S., Grahn, T., Olsson, L., García-Gil, M. R., Karlsson, B., & Seifert, T. (2018). Age and weather effects on between and within ring variations of number, width and coarseness of tracheids and radial growth of young Norway spruce. *European Journal of Forest Research*, **137**, 719-743.

Martínez-Sancho, E., Rellstab, C., Guillaume, F., Bigler, C., Fonti, P., Wohlgemuth, T., & Vitasse, Y. (2021). Post-glacial re-colonization and natural selection have shaped growth responses of silver fir across Europe. *Science of the Total Environment*, **779**, 146393.

Milesi, P., Berlin, M., Chen, J., Orsucci, M., Li, L., Jansson, G., Karlsson, B., & Lascoux, M. (2019). Assessing the potential for assisted gene flow using past introduction of Norway spruce in southern Sweden: Local adaptation and genetic basis of quantitative traits in trees. *Evolutionary Applications*, **12**(10), 1946-1959.

Miriti, M. N. (2006). Ontogenetic shift from facilitation to competition in a desert shrub. *Journal of Ecology*, **94**(5), 973-979.

Montwé, D., Isaac-Renton, M., Hamann, A., & Spiecker, H. (2016). Drought tolerance and growth in populations of a wide-ranging tree species indicate climate change risks for the boreal north. *Global Change Biology*, **22**(2), 806-815.

Niklas, K. J. (1993). Influence of tissue density-specific mechanical properties on the scaling of plant height. *Annals of Botany*, **72**(2), 173-179.

Nota, K., Klaminder, J., Milesi, P., Bindler, R., Nobile, A., van Steijn, T., Bertilsson, S., Svensson, B., Hirota, S. K., Matsuo, A., Gunnarsson, U., Seppä, H., Väliranta, M. M., Wohlfarth, B., Suyama, Y., & Parducci, L. (2022). Norway spruce postglacial recolonization of Fennoscandia. *Nature Communications*, **13**(1), 1-9.

Nystedt, B., Street, N. R., Wetterbom, A., Zuccolo, A., Lin, Y. C., Scofield, D. G., Vezzi, F., Delhomme, N., Giacomello, S., Alexeyenko, A., Vicedomini, R., Sahlin, K., Sherwood, E., Elfstrand, M., Gramzow, L., Holmberg, K., Hällman, J., Keech, O., Klasson, L., ... Jansson, S. (2013). The Norway spruce genome sequence and conifer genome evolution. *Nature*, **497**, 579-584.

Palmé, A. E., Su, Q., Rautenberg, A., Manni, F., & Lascoux, M. (2003). Postglacial recolonization and cpdn variation of silver birch, *Betula pendula*. *Molecular Ecology*, **12**, 201-212.

Petit, R. J., Aguinagalde, I., de Beaulieu, J. L., Bittkau, C., Brewer, S., Cheddadi, R., Ennos, R., Fineschi, S., Grivet, D., Lascoux, M., Mohanty, A., Müller-Starck, G., Demesure-Musch, B., Palmé, A., Martín, J. P., Rendell, S., & Vendramin, G. G. (2003). *Glacial refugia*: Hotspots but not melting pots of genetic diversity. *Science*, **300**(5625), 1563-1565.

Petit, R. J., Brewer, S., Bordács, S., Burg, K., Cheddadi, R., Coart, E., Cottrell, J., Csaikl, U. M., van Dam, B., Deans, J. D., Espinel, S., Fineschi, S., Finkeldey, R., Glaz, I., Goicoechea, P. G., Jensen, J. S., König, A. O., Lowe, A. J., Madsen, S. F., ... Kremer, A. (2002). Identification of refugia and post-glacial colonisation routes of European white oaks based on chloroplast DNA and fossil pollen evidence. *Forest Ecology and Management*, **156**(1-3), 49-74.

Poorter, L. (2008). The relationships of wood-, gas-and water fractions of tree stems to performance and life history variation in tropical trees. *Annals of Botany*, **102**(3), 367-375.

Poorter, L., Bongers, L., & Bongers, F. (2006). Architecture of 54 moist-forest tree species: Traits, trade-offs, and functional groups. *Ecology*, **87**(5), 1289-1301.

Poorter, L., McDonald, I., Alarcón, A., Fichtler, E., Licona, J. C., Peña-Claros, M., Sterck, F., Villegas, Z., & Sass-Klaassen, U. (2010). The importance of wood traits and hydraulic conductance for the performance and life history strategies of 42 rainforest tree species. *New Phytologist*, **185**(2), 481-492.

Preston, K. A., Cornwell, W. K., & DeNoyer, J. L. (2006). Wood density and vessel traits as distinct correlates of ecological strategy in 51 California coast range angiosperms. *New Phytologist*, **170**(4), 807-818.

Pretzsch, H., Biber, P., Schütze, G., Kemmerer, J., & Uhl, E. (2018). Wood density reduced while wood volume growth accelerated in central European forests since 1870. *Forest Ecology and Management*, **429**, 589-616.

R Core Team. (2021). *R: A language and environment for statistical computing*. R Foundation for Statistical Computing. https://www.R-project.org/
Rossi, S., Cairo, E., Krause, C., & Deslauriers, A. (2015). Growth and basic wood properties of black spruce along an alti-latitudinal gradient in Quebec, Canada. *Annals of Forest Science*, 72, 77-87.

Saleh, D., Chen, J., Leplé, J. C., Leroy, T., Truffaut, L., Lalanne, C., Labadie, K., Lesur, I., Bert, D., Lagane, F., Morneau, F., Aury, J. M., Plomion, C., Lascoux, M., & Kremer, A. (2022). Genome-wide evolutionary response of European oaks during the Anthropocene. *Evolution letters*, 6(1), 4-20.

Sang, Z., Sebastian-Azcona, J., Hamann, A., Menzel, A., & Hacke, U. (2019). Adaptive limitations of white spruce populations to drought imply vulnerability to climate change in its western range. *Evolutionary Applications*, 12(9), 1850-1860.

Savolainen, O., Lascoux, M., & Merilä, J. (2013). Ecological genomics of local adaptation. *Nature Reviews Genetics*, 14, 807-820.

Savolainen, O., Pyhjäjärvi, T., & Knürr, T. (2007). Gene flow and local adaptation in trees. *Annual Review of Ecology, Evolution, and Systematics*, 38, 595-619.

Spitze, K. (1993). Population structure in *Daphnia obtusa*: Quantitative genetic and allozymic variation. *Genetics*, 135(2), 367-374.

Suarez, A. V., & Tsutsui, N. D. (2008). The evolutionary consequences of biological invasions. *Molecular Ecology*, 17(1), 351-360.

Tsuda, Y., Chen, J., Stocks, M., Källman, T., Sønstebø, J. H., Parducci, L., Semerikov, V., Sperisen, C., Politov, D., Ronkainen, T., Väliranta, M., Vendramin, G. G., Tollefsrud, M. M., & Lascoux, M. (2016). The extent and meaning of hybridization and introgression between Siberian spruce (Picea obovata) and Norway spruce (Picea abies): Cryptic refugia as stepping stones to the west? *Molecular Ecology*, 25(12), 2773-2789.

Van Gelder, H. A., Poorter, L., & Sterck, F. J. (2006). Wood mechanics, allometry, and life-history variation in a tropical rain forest tree community. *New Phytologist*, 171(2), 367-378.

Vendramin, G. G., Anzidei, M., Madaghiele, A., Sperisen, C., & Bucci, G. (2000). Chloroplast microsatellite analysis reveals the presence of population subdivision in Norway spruce (Picea abies K.). *Genome*, 43(1), 68-78.

Vidalis, A., Scofield, D. G., Neves, L. G., Bernhardsson, C., García-Gil, M. R., & Ingvarsson, P. K. (2018). Design and evaluation of a large sequence-capture probe set and associated SNPs for diploid and haploid samples of Norway spruce (Picea abies). *BioRxiv*, 291716.

Von Gadow, K., Álvarez González, J. G., Zhang, C., Pukkala, T., & Zhao, X. (2021). The shape and growth of Forest trees. In *Sustaining Forest ecosystems* (pp. 159-219). Springer.

Whitlock, M. C. (2008). Evolutionary inference from QST. *Molecular Ecology*, 17(8), 1885-1896.

Wisz, M. S., Pottier, J., Kissling, W. D., Pellissier, L., Lenoir, J., Damgaard, C. F., Dormann, C. F., Forchhammer, M. C., Grytnes, J. A., Guisan, A., Heikkinen, R. K., Høye, T. T., Kühn, I., Luoto, M., Maiorano, L., Nilsson, M. C., Normand, S., Öckinger, E., Schmidt, N. M., ... Svenning, J. C. (2013). The role of biotic interactions in shaping distributions and realised assemblages of species: Implications for species distribution modelling. *Biological Reviews*, 88, 15-30.

## SUPPORTING INFORMATION

Additional supporting information can be found online in the Supporting Information section at the end of this article.

**How to cite this article:** Tiret, M., Olsson, L., Grahn, T., Karlsson, B., Milesi, P., Lascoux, M., Lundqvist, S.-O., & García-Gil, M. R. (2023). Divergent selection predating the Last Glacial Maximum mainly acted on macro-phenotypes in Norway spruce. *Evolutionary Applications*, 16, 163-172. https://doi.org/10.1111/eva.13519