# Reduced within-population quantitative genetic variation is associated with climate harshness in maritime pine

Juliette Archambeau, Marta Benito-Garzón, Marina de Miguel, Benjamin Brachi, Frederic Barraquand, Santiago C. Gonzalez-Martinez

► To cite this version:

Juliette Archambeau, Marta Benito-Garzón, Marina de Miguel, Benjamin Brachi, Frederic Barraquand, et al.. Reduced within-population quantitative genetic variation is associated with climate harshness in maritime pine. Heredity, inPress, 10.1038/s41437-023-00622-9 . hal-04101611

HAL Id: hal-04101611

https://hal.inrae.fr/hal-04101611

Submitted on 20 May 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>
Reduced within-population quantitative genetic
variation is associated with climate harshness in
maritime pine

Juliette Archambeau¹,², Marta Benito Garzón¹, Marina de Miguel¹,³,
Benjamin Brachi¹, Frédéric Barraquand⁴ and Santiago C.
González-Martínez¹

¹INRAE, Univ. Bordeaux, BIOGECO, F-33610 Cestas, France

²UK Centre for Ecology & Hydrology, Bush Estate, Penicuik, United
Kingdom

³EGFV, Univ. Bordeaux, Bordeaux Sciences Agro, INRAE, ISVV, F-33882,
Villenave d'Ornon, France

⁴CNRS, Institute of Mathematics of Bordeaux, F-33400 Talence, France

**Corresponding author:**

Juliette Archambeau
juli.archambeau@orange.fr

**Authors' ORCIDs:**

Juliette Archambeau: https://orcid.org/0000-0003-0155-9538

Marta Benito Garzón: https://orcid.org/0000-0002-3436-123X

Marina de Miguel: https://orcid.org/0000-0001-6398-2660

Benjamin Brachi: https://orcid.org/0000-0001-5988-7150

Frédéric Barraquand: https://orcid.org/0000-0002-4759-0269

Santiago C. González-Martínez: https://orcid.org/0000-0002-4534-3766

**Running title:** Evolutionary forces maintaining genetic variation

**Keywords:** quantitative genetic variation, adaptive potential, forest tree, natural selection, environmental heterogeneity, cold winter temperatures.

**Acknowledgments:** We thank A. Saldaña, F. del Caño, E. Ballesteros and D. Barba (INIA) and the ‘Unité Expérimentale Forêt Pierroton’ (UEFP, INRAE; https://doi.org/10.15454/1.5483264699193726E12) for field assistance (plantation and measurements). Data used in this
research are part of the Spanish Network of Genetic Trials (GENFORED, http://www.genfordes.es). We thank all persons and institutions linked to the establishment and maintenance of field trials used in this study. We are very grateful to Ricardo Alía and Juan Majada who initiated and supervised the establishment of the CLONAPIN network and provided the progeny test height data that we used in the validation analysis. We thank Maurizio Marchi for kindly providing the raster files of the climate variables corresponding to the desired time period and spatial extent and for providing detailed explanations in response to our questions regarding the extraction of climatic data with the Climatic Downscaling Tool. JA was funded by the University of Bordeaux (ministerial grant). This study was funded by the 'Initiative d'Excellence (IdEx) de l'Université de Bordeaux - Chaires d'installation 2015' (EcoGenPin) and the European Union's Horizon 2020 research and innovation programme under grant agreements No 773383 (B4EST) and No 862221 (FORGENIUS).

**Author contributions:** SCG-M designed the experiment and supervised the curation of field data. MdM cleaned and formatted the phenotypic data. SCG-M, JA, FB, MBG and BB conceived the paper methodology. JA and FB built the model equations and codes. JA conducted the data and simulation analyses. All authors interpreted the results. JA led the writing of the manuscript. All authors contributed to the manuscript and gave final approval for publication.

**Data and script availability:** Data are publicly available. SNP data were deposited in the Dryad repository at http://dx.doi.org/10.5061/dryad.8d6k1. Height data have been deposited in GENFORED, the Spanish Network of Genetic Trials (http://www.genfordes.es). The scripts are available in the following Github repository: https://github.com/JulietteArchambeau/H2Pinpin.

**Preregistration:** The objectives and methods of this manuscript have been pre-registered at the Center For Open Science and are available here (anonymized repository): https://osf.io/knx6z/?view_only=41bb7b5cbf7241d0856e8b9e393cc795.

## Abstract

How evolutionary forces interact to maintain genetic variation within populations has been a matter of extensive theoretical debates. While mutation and exogenous gene flow increase genetic variation, stabilizing selection and genetic drift are expected to deplete it. To date, levels of genetic variation observed in natural populations are hard to predict without accounting for other processes, such as balancing selection in heterogeneous environments. We aimed to empirically test three hypotheses: (i) admixed populations have higher quantitative genetic variation due to introgression from other gene pools, (ii) quantitative genetic variation is lower in populations from harsher environments (i.e. experiencing stronger selection), and (iii) quantitative genetic variation is higher in populations from heterogeneous environments. Using growth, phenological and functional trait data from three clonal common gardens and 33 populations (522 clones) of maritime pine (Pinus pinaster Aiton), we estimated the asso-
ciation between the population-specific total genetic variances (i.e. among-clone variances) for these traits and ten population-specific indices related to admixture levels (estimated based on 5,165 SNPs), environmental temporal and spatial heterogeneity and climate harshness. Populations experiencing colder winters showed consistently lower genetic variation for early height growth (a fitness-related trait in forest trees) in the three common gardens. Within-population quantitative genetic variation was not associated with environmental heterogeneity or population admixture for any trait. Our results provide empirical support for the potential role of natural selection in reducing genetic variation for early height growth within populations, which indirectly gives insight into the adaptive potential of populations to changing environments.

# 1 Introduction

Most complex traits show marked heritable variation in natural populations. How evolutionary forces interact to maintain such variation remains a long-standing dilemma in evolutionary biology and quantitative genetics (Johnson and Barton 2005). While mutation and genetic drift have straightforward roles, generating and eliminating variation respectively, the effect of natural selection is more complicated (Walsh and Lynch 2018). Stabilizing selection, i.e. the selection of intermediate phenotypes, is often strong in natural populations (Hereford et al. 2004). This type of selection is expected to deplete genetic variation (Fisher 1930), either directly on the focal trait or indirectly via pleiotropic effects (Johnson and Barton 2005). Theoretical models based on the balance between mutation, drift and stabilizing selection support this idea, but they suggest lower fitness heritability values than those generally observed in empirical studies (Johnson and Barton 2005). Balancing selection encompasses various evolutionary processes that can maintain greater than neutral genetic variation within populations (Mitchell-Olds et al. 2007), such as temporally or spatially fluctuating selection pressures in subpopulations connected by gene flow (Felsenstein 1976). Since Levene's archetypal model (1953), a large corpus of single-locus and polygenic models, most often deterministic, have generally concluded that genetic polymorphisms in spatially heterogeneous environments can only be maintained under restrictive conditions (Spichtig and Kawecki 2004, Byers 2005). In this line, McDonald and Yeaman (2018) showed with stochastic individual-based simulations that substantial within-population genetic variation can be maintained in spatially (and to a much lesser extent, temporally) heterogeneous environments, at intermediate migration rates and regardless of population size.

In empirical studies, the causes that maintain genetic variation in populations have been largely studied based on genetic and discrete-trait polymorphisms (reviewed in Hedrick 1986, 2006), but much less via quantitative genetic analysis of complex traits. Lower quantitative genetic variation in populations evolving under strong selection pressures has been found in a number of studies, e.g. for flowering time and size at reproduction in an invasive plant species (Colautti et al. 2010), for desiccation resistance in *Drosophila birchii* (van Heerwaarden et al. 2009) and for body mass in great tits (Charmantier et al. 2004), but not in others (e.g. Merilä et al. 2004, Stock et al. 2014). Higher genetic variation in populations evolving under
spatially or temporally varying selection pressures has been found in *Drosophila* populations for body mass, survival and the sternopleural chaeta number (Mackay 1981, Huang et al. 2015), but not for wing shape (Yeaman et al. 2010). The lack of general trends in empirical studies can be explained by method-specific pitfalls to accurately estimate quantitative genetic variation, e.g. the genetic and environmental variances are hard to disentangle in the wild, and when estimated in common gardens, their environment-dependent nature does not allow for wide generalization of estimates (Hoffmann and Parsons 1991, Merilä et al. 2001, Charmantier et al. 2004).

Forest trees have specific life-history traits and genomic features making them interesting model species in population and quantitative genetic studies (Petit and Hampe 2006, Savolainen et al. 2007). Compared to crop species, they remain largely undomesticated (Neale and Savolainen 2004). Most forest trees are outcrossing, have high lifetime reproductive output and long generation times. They often display important gene flow among populations through long-distance pollen dispersal (Kremer et al. 2012). They show slow rates of macroevolution (i.e. low nucleotide substitution rates and low speciation rates; Petit and Hampe 2006), generally have large effective population sizes, with distributions often covering a wide range of environmental conditions (Alberto et al. 2013). Extensive work has revealed strong clines at large geographical scales in the population-specific mean values of phenotypic traits (reviewed in Savolainen et al. 2007, Benito Garzón et al. 2019), e.g. phenological traits with latitude or altitude (Alberto et al. 2011, Thibault et al. 2020) or height growth with cold hardiness (Rehfeldt et al. 1999, Leites et al. 2012). Genetic differentiation at microgeographic spatial scales has also been repeatedly observed (reviewed in Linhart and Grant 1996, Jump and Peñuelas 2005, Scotti et al. 2016), suggesting rapid rates of microevolution (Petit and Hampe 2006). Possible explanations include the fact that forest trees have high levels of genetic diversity and that most of their quantitative and neutral genetic variation is within populations (Hamrick 2004). To our knowledge, only two empirical studies investigated the potential causes underlying the maintenance of quantitative trait variation within forest tree populations. Yeaman and Jarvis (2006) showed that 20% of growth genetic variation in lodgepole pine populations was attributable to regional spatial heterogeneity, suggesting an important role of gene flow and varying selection pressures. In the neotropical oak *Quercus oleoides*, Ramírez-Valiente et al. (2019) found lower quantitative genetic variation of functional traits in harsher environments, but not higher quantitative genetic variation in temporally fluctuating environments. They also suggested only a marginal effect of genetic structure and diversity on the maintenance of within-population genetic variation. More empirical studies on the effects of evolutionary forces on quantitative genetic variation across natural forest tree populations are therefore needed to assess their evolutionary potential, anticipate their responses to ongoing global change, and develop sound adaptive management and conservation strategies that guarantee their future persistence.

In this study, we aimed to test three competing, but not mutually exclusive, hypotheses regarding the potential drivers maintaining or reducing quantitative genetic variation within maritime pine populations, namely: i) the most admixed populations have higher quantitative genetic variation due to exogenous gene flow from other gene pools, and this relationship is proportional to the divergence between sink and source gene pools; ii) quantitative genetic
variation is lower in populations that have evolved in harsher environments, as a result of
higher selection pressures in these regions; and iii) quantitative genetic variation is higher
in populations that have evolved in spatially or temporally heterogeneous environments.
Importantly, the last two hypotheses require the action of natural selection, while the first
does not. Therefore, under the last two hypotheses, detecting changes in genetic variances
will be more likely for traits more closely related to fitness, whereas under the first hypoth-
esis, detecting changes in genetic variances will be equally likely across traits. We used
growth (four height measurements), phenological (bud burst and duration of bud burst)
and functional ($\delta^{13}$C and specific leaf area, SLA) trait data from three clonal common gar-
dens, consisting of 522 clones (i.e. genotypes) from 33 populations, spanning all known gene
pools in the species (Jaramillo-Correa et al. 2015) and genotyped for 5,165 SNPs. We calcu-
lated ten population-specific indices related to our three hypotheses: two indexes describing
the level and origin of admixture in the populations (estimated with SNP markers), two
indexes of climate harshness (capturing summer droughts and winter cold temperatures),
four indexes describing the environmental spatial heterogeneity in the forested areas sur-
rounding the populations and two indexes describing the climatic temporal heterogeneity
(i.e. annual variability in summer droughts and winter cold temperatures). We then used
quantitative genetic models in a Bayesian framework to estimate for each trait the linear
association among the population-specific total genetic variances (i.e. hereafter referred to as
the within-population quantitative genetic variation) and each of the ten population-specific
indices described above.

2 Materials & Methods

## 2.1 Maritime pine, a forest tree growing in heterogeneous environments

Maritime pine (*Pinus pinaster* Ait., Pinaceae) is a wind-pollinated, outcrossing and long-lived tree species with large ecological and economical importance in western Europe and North Africa. Maritime pine is largely appreciated for its wood, for stabilizing coastal and fossil dunes and, as a keystone species, for supporting biodiversity (Viñas et al. 2016). The distribution of maritime pine natural populations is scattered and covers a wide range of environmental conditions. Several studies have provided evidence of genetic differentiation for adaptive traits in this species, suggesting local adaptation (e.g. González-Martínez et al. 2002, de Miguel et al. 2022). Maritime pine can grow in widely different climates: the dry climate along the northern coasts of the Mediterranean Basin (from Portugal to western Italy), the mountainous climates of south-eastern Spain and Morocco, the wetter climate of the Atlantic region (from the Spanish Iberian region to the western part of France) and the continental climate of central Spain. Maritime pine can also grow on a wide range of substrates, from sandy and acidic soils to more calcareous ones (Viñas et al. 2016).

Maritime pine presents a strong population genetic structure with six main gene pools,
located in the French Atlantic region, Iberian Atlantic region, central Spain, south-eastern Spain, Corsica and northern Africa (Fig. S10; Alberto et al. 2013, Jaramillo-Correa et al. 2015). These gene pools consist of homogeneous (albeit interconnected by gene flow) genetic clusters, as shown by molecular markers (Jaramillo-Correa et al. 2015), and probably resulted from the expansion of distinct glacial refugia (Bucci et al. 2007).

## 2.2 Phenotypic data

Phenotypic data were collected in three clonal common gardens planted in 2011 (Fig. 1): Asturias (Spain) and Fundão (Portugal) in the Iberian Atlantic region and Bordeaux (France) in the French Atlantic region. As evidenced by the high survival rate at these sites (Table S1), the common gardens are located in environments considered favorable to maritime pine, with mild winters, high annual rainfall and relatively wet summers (Tables S11, S12 and Fig. 1). In each of these common gardens, trees belonging to 522 clones (i.e. genotypes) from 33 populations, covering the six known gene pools in the species (Fig. S10), were planted following a randomized complete block design with 8 blocks. To obtain the clones, trees at least 50 m apart were sampled in natural stands, and one seed per tree was planted in a nursery and vegetatively propagated by cuttings (see Rodríguez-Quilón et al. 2016 for details). Clones were therefore considered unrelated.

One growth trait, height, was measured in all common gardens and at different tree ages (Table S1). Two phenology-related traits, the mean bud burst date over four years and the mean duration of bud burst over three years, were measured in Bordeaux and were averaged over several years to suppress differences across years and approximate a normal distribution of their trait values (Table S1). Bud burst corresponds to the date of brachyblast emergence in accumulated degree-days (with base temperature 0°C) from the first day of the year to account for between-year variability in temperature. The duration of bud burst corresponds to the number of degree-days between the beginning of bud elongation and the total elongation of the needles (see Hurel et al. 2021). Last, two functional traits, δ¹³C and the specific leaf area (SLA) were measured in Fundão (Table S1). For each trait and each population, the number of trees and clones used in the analyses, the raw phenotypic means and variances, and histograms of trait distribution are provided in section 1.1 of the Supplementary Information. Prior to analyses, some traits were log-transformed to get closer to normality (Fig. S1) or mean-centered to help model convergence (Table S1).

## 2.3 SNP genotyping and population admixture

The 522 clones planted in the Asturias common garden were genotyped with the Illumina Infinium assay described in Plomion et al. (2016), resulting in 5,165 high-quality polymorphic SNPs. There were on average only 3.3 missing SNP values per genotype (ranging between 0 and 142). For each clone, the proportions of ancestry from each of the six known gene pools were estimated in Jaramillo-Correa et al. (2015) using the Bayesian approach available in
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Locations of the three common gardens and the 33 populations used in the study. The colors represent the gradients of mean coldest month temperature (MCMT) and summer (May-September) precipitation (SP) over the period 1901-1950 within the maritime pine range. These gradients were obtained by performing a centered and scaled principal component analysis on the MCMT and SP values. The maritime pine distribution combines the EUFORGEN distribution (http://www.euforgen.org/) and 10-km radius areas around the National Forest Inventory plots with this species.</figcaption></figure>

Structure (Pritchard et al. 2000), and were then averaged by population (Table S14). For each population, the main gene pool was defined as the gene pool constituting the higher proportion of population ancestry whereas the other gene pools were considered as a product of (past or present) exogenous gene flow. The accuracy of the assignment of each population to a given gene pool was validated by observing that populations with the same main gene pool ancestry were located close to each other. First, we calculated a population admixture score A, as the proportion of ancestry from foreign gene pools (Table S14). Second, we calculated a population admixture score D that considers both the proportion of foreign ancestries and the divergence between the main and foreign gene pools (Table S14). For that, we weighted the proportions of ancestry from foreign gene pools by the sum of the allele frequency divergence of the main and foreign gene pool from the common ancestral one ($F_k$, which should be numerically similar to $F_{ST}$; Falush et al. 2003). We developed D considering that some gene pools are more divergent than others and thus may bring higher genetic diversity to an admixed population at the same level of introgression. A was highly correlated with D (Pearson correlation coefficient of 0.91; Table S14), and also with a related index that used weights based on pairwise $F_{ST}$ (Table S14). A and D were mean-centered
and divided by their standard deviation prior to analyses.

## 2.4 Climate harshness and temporal and spatial heterogeneity

To describe the winter temperature and drought summer conditions under which the populations have evolved, we extracted the annual values of the mean coldest month temperature (MCMT) and the summer (May-September) precipitation (SP) at the population location for the period 1901-1950 using the bilinear interpolation of the Climate Downscaling Tool (https://www.ibbr.cnr.it/climate-dt/). SP was selected to better understand how trees adapt in areas subjected to frequent summer droughts, which is key to anticipate tree population responses to ongoing climate change, while MCMT was chosen because marked patterns of cold adaptation have been previously reported in maritime pine (e.g. Grivet et al. 2011). The means of the annual values of MCMT and SP over the period considered were used as two indexes of climate harshness while their variances were used as two indexes of temporal heterogeneity.

In addition, to describe the environmental spatial heterogeneity around the population location, we used raster files of climatic, topographic and soil data at 1-km resolution. Climatic raster data correspond to climatic variables from the *ClimateEU* database (Marchi et al. 2020) and averaged over the period 1901-1950. Topographic raster data were generated from NASA's Shuttle Radar Topography Mission (SRTM) and SAGA v 2.3.1 (Conrad et al. 2015), and were used to calculate the topographic ruggedness index (TRI), which quantifies the terrain heterogeneity, i.e. differences in elevation among adjacent cells (Riley et al. 1999). Soil raster data were downloaded from the European Soil Database (Hiederer et al. 2013). All environmental variables used are listed in Table S15.

We extracted raster cell values of the environmental variables within a 20-km radius around each population location, and kept only raster cells that fell within forested areas, to avoid including environmental data from non-suitable areas (e.g. lakes, mountain peaks; section 1.4.2 of the Supplementary Information). We then performed a principal component analysis (PCA) on the set of the selected raster cells and extracted the PC1 and PC2 scores of each cell, accounting for 45.2% and 34.1% of the variance, respectively (Fig. S11). PC1 mainly captured precipitation and maximum temperature conditions in summer, while PC2 primarily represented minimum temperature conditions in winter. To obtain the indexes of spatial environmental heterogeneity, we calculated the variances of the PC1 and PC2 scores in a 20-km and 1.6-km radius around each population location (Figs S12 & S13). The environmental heterogeneity indexes were only very weakly correlated (Pearson correlation coefficients lower than 0.36) with the number of forested cells (i.e. the area considered to calculate the indexes), ensuring that the estimated effects of spatial heterogeneity in further analyses were not due to the area per se (Triantis et al. 2003, Stein et al. 2014).

The population-specific values of the spatial and temporal heterogeneity and climate harshness indexes are shown in Tables S16 and S17. All ten indexes were mean-centered and divided by their standard deviation prior to analyses.
## 2.5 Univariate quantitative genetic models

We constructed ten models for each trait y, that is one model per potential driver of the within-population genetic variance (i.e. the two admixture scores, the six environmental heterogeneity indexes and the two climate harshness indexes). Each model was derived from the following baseline quantitative genetic model:

$$y_{bpcr} = \mathcal{N}(\beta_0 + B_b + P_p + C_{c(p)}, \sigma_r^2) \quad (1)$$

$\sigma_r^2$ is the residual variance. $\beta_0$ is the global intercept with a weakly informative prior; $\beta_0 \sim \mathcal{N}(\bar{y}, 2)$, $\bar{y}$ being the mean of the observed trait values. $B_b$, $P_p$ and $C_{c(p)}$ are the block, population and clone (nested within population) varying intercepts; they are distributed as follows:

$$\begin{gathered} \begin{bmatrix} B_b \\ P_p \end{bmatrix} \sim \mathcal{N} \left( 0, \begin{bmatrix} \sigma_B^2 \\ \sigma_P^2 \end{bmatrix} \right) \\ C_{c(p)} \sim \mathcal{N}(0, \sigma_{C_p}^2) \end{gathered} \quad (2)$$

where $\sigma_B^2$ and $\sigma_P^2$ are the variance among blocks and populations, respectively, and $\sigma_{C_p}^2$ are the population-specific variances among clones (i.e. population-specific total genetic variances), which we refer to herein as the within-population quantitative genetic variation. The total variance $\sigma_{tot}^2$ is partitioned as follows:

$$\begin{align} \sigma_{tot}^2 &= \sigma_r^2 + \sigma_B^2 + \overline{\sigma_{C_p}^2} + \sigma_P^2 \nonumber \\ \sigma_r &= \sigma_{tot} \times \sqrt{(\pi_r)} \nonumber \\ \sigma_B &= \sigma_{tot} \times \sqrt{(\pi_B)} \nonumber \\ \sigma_P &= \sigma_{tot} \times \sqrt{(\pi_P)} \tag{3} \\ \overline{\sigma_{C_p}} &= \sigma_{tot} \times \sqrt{(\pi_C)} \nonumber \\ 1 &= \pi_r + \pi_B + \pi_C + \pi_P \nonumber \\ \sigma_{tot} &\sim \mathcal{S}^*(0, 1, 3) \end{align}$$

where $\overline{\sigma_{C_p}^2}$ is the arithmetic mean of the population-specific among-clone variances ($\sigma_{C_p}^2$) and $\overline{\sigma_{C_p}}$ is the quadratic mean of the standard deviations. $\sigma_{C_p}$ is expressed as follows:

$$\sigma_{C_p} \sim \mathcal{LN} \left( \ln(\overline{\sigma_{C_p}}) - \frac{\sigma_K^2}{2} + \beta_X X_p, \sigma_K^2 \right) \quad (4)$$
$\sigma_K^2$ is the variance of the population-specific among-clone standard deviations $\sigma_{C_p}$, with $\sigma_K \sim \exp(1)$ (see section 2 in the Supplementary Information for more details). $X_p$ is one of the ten potential drivers and $\beta_X$ stands for the linear association between $\sigma_{C_p}$ and the potential driver considered, with $\beta_X \sim N(0, 1)$.

The modeling approach described above has two key advantages. First, estimating the total genetic variances for each trait (i.e. $\sigma_{C_p}^2$ coefficients) and their association with the ten potential drivers (i.e. $\beta_X$ coefficients) in a one-step quantitative genetic model within a Bayesian framework allowed us to propagate the uncertainty in the parameter estimates. Second, we used weakly informative priors, which are a mean to statistically regularize the parameters and thus avoid getting aberrant values (Lemoine 2019). These priors are particularly useful to deal with small sample sizes and therefore we were able to estimate the total genetic variances for each population despite the small number of clones in some of them. We evaluated the model ability to recover the true parameter values under the assumption that the mathematical model represents the biological processes at stake (i.e. how the real data was generated). For each common garden, we generated 100 simulated height samples with the same sample sizes as the real experimental design (see Tables S3, S4 and S5), and based on the mean and variance of the observed height data. We extracted the mean standard error and bias error of the estimates, and the coverage of the 80% and 95% credible intervals, to compare the model estimates with the true parameter values (see details in section 3 of the Supplementary Information).

It can be noted that the estimates of the total genetic variances may be inflated by C effects (environmental effects common to a group of relatives such as clones, Lerner 1958). However, we expect C effects to inflate the total genetic variances uniformly across populations and are therefore unlikely to bias the $\beta_X$ estimates.

Model specification and fit were performed using the Stan probabilistic programming language (Carpenter et al. 2017), based on the no-U-turn sampler algorithm. Models were run with four chains and 2,500 iterations per chain (including 1,250 warm-up samples not used for the inference). All analyses were undertaken in R version 3.6.3 (R Core Team 2020) and scripts are available at https://github.com/JulietteArchambeau/H2Pinpin.

## 2.6 Validation analysis for height and SLA in a family-based progeny test

To validate our results, we used height and SLA measurements from a progeny test near Asturias (thus in a similar environment to the nearby clonal common garden of our study) which shares 23 populations with the CLONAPIN common garden network (progeny test data kindly provided by Dr. Ricardo Alia, CSIC, Madrid). As the progeny test is based on families, we were able to estimate the additive genetic variance within populations (and not only the total genetic variance) for height at 3 and 6 years and SLA at 5 years (see section 10 of the Supplementary Information for more details).
# 3 Results

When comparing the estimated height values to the true values ($\sigma_K = 0.1$ and $\beta_X = 0.1$) over the 100 simulated samples, the mean standard error was around 0.06 for $\sigma_K$ and 0.05 for $\beta_X$, the mean bias error was around 0.01 for $\sigma_K$ and between -0.004 and 0.006 for $\beta_X$, the coverage of the 80% credible interval was higher than 93% for $\sigma_K$ and higher than 76% for $\beta_X$, and the coverage of the 95% credible interval was higher than 99% for $\sigma_K$ and 96% for $\beta_X$ (Table S19, see also Figs S14-S16). These simulations therefore showed that, under the assumption that the statistical model reflects the processes at work, our model displayed a satisfactory accuracy to be used in the following analyses.

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Median and 95% credible intervals of the $\beta_X$ coefficients, which stand for the association among within-population quantitative genetic variation and its potential underlying drivers on the x-axis (see equation 4). A and D are the two admixture scores; mean(MCMT) and mean(SP) are the two indexes of climate harshness and correspond to the mean of the annual values of the mean coldest month temperature (MCMT) and the summer precipitation (SP) over the period 1901-1950, respectively; var(MCMT) and var(SP) are the two indexes of temporal heterogeneity and correspond to the variance of the annual values of MCMT and SP over the period 1901-1950; finally, SH1[20km], SH2[20km], SH1[1.6km] and SH2[1.6km] are the indexes of spatial heterogeneity in a 20-km and 1.6-km radius around the population location calculated as the variance of the first two principal components of a PCA summarizing the environmental variation (see details in section 1.4.2 of the Supplementary Information). Colors stand for the different traits under study and the shapes for the different types of traits, i.e. height (circles), phenology-related traits (triangles) and functional traits (squares).</figcaption></figure>

The proportion of variance explained by the models (i.e. the sum of the among-population, among-clone and among-block variances) and the variance partitioning varied broadly across traits (see section 5 in the Supplementary Information). More specifically, the models explained between 40% and 50% of the variance for phenology-related traits, between 30% and 40% for functional traits, and from 20% for height in Fundão to almost 60% for height in Bordeaux at 85-month old (Fig. S17). Residual variance explained most of the variance for all traits, except for height in Bordeaux at 85-month old, where ~40% of the variance
came from variation among populations, ~40% from residuals and the remaining ~20% from variation among clones (Fig. S29). Variation among populations was higher than variation among clones for height and δ¹³C (Figs. S21, S25, S29, S33 and S49), but not for SLA and phenology-related traits (Figs. S45, S37 and S41).

We found that populations experiencing colder winter conditions (i.e. lower MCMT) consistently displayed lower quantitative genetic variation for early height growth in the three common gardens (Fig. 2). Holding all other parameters constant, a one-standard deviation decrease in the mean temperature of the coldest month over the period considered was associated, on average, with a 42%, 21% and 25% decrease of σ<sub>Cp</sub> for height in Fundão, Bordeaux at 25-month old and Asturias, respectively (see details of the calculation in the section 4 of the Supplementary Information). Interestingly, this association was not found in Bordeaux when the trees were older (i.e. at 85-month old; Fig. 2). The negative association between colder winter conditions and quantitative genetic variation was also detected in a model combining height data from the three common gardens (except Bordeaux at 85-month old; Fig. S52), and in which the genotype-by-environment interaction explained only 3.7% of the phenotypic variance (Fig. S54). Unexpectedly, populations experiencing colder winter conditions also showed higher genetic variation for SLA (Fig. 2). The detected associations were unlikely to have been biased by unequal sample sizes across clones and populations, as detailed in section 11 of the Supplementary Information.

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Validation step using height and SLA measurements from a family-based progeny test near Asturias (data kindly provided by Dr. Ricardo Alía, CSIC, Madrid). Median and 95% credible intervals of the $β_X$ posterior distributions are shown. In the validation analysis, $β_X$ coefficients stand for the association among within-population additive genetic variation and its potential underlying drivers on the x-axis. A description of the drivers can be found in the legend of Fig. 2.</figcaption></figure>

Importantly, in the validation analysis using a family-based progeny test, we also found that populations experiencing colder winter conditions have lower additive genetic variation for height at 3-year old, but not at 6-year old, and we did not find any association with the other potential drivers (Fig. 3). In addition, the negative association between cold
winter temperatures and the within-population quantitative genetic variation for SLA was
not replicated in the validation analysis (Fig. 3). Therefore, we consider our results to
provide only robust support for a decrease in genetic variation for early height growth under
colder winter conditions. The other detected associations for early height growth were not
consistently replicated across the three common gardens or in the validation analysis (Figs.
2 and 3), and thus we do not interpret them in the discussion.

# 4 Discussion

How quantitative genetic variation is maintained within populations remains a long-standing
open question that has been extensively explored in theoretical work but still lacks sufficient
empirical evidence (Johnson and Barton 2005). We found that genetic variation for early
height growth in maritime pine is lower in populations exposed to colder winters, which
is consistent with some theoretical models predicting that quantitative genetic variation in
fitness-related traits is lower in populations under strong selection (Fisher 1930). However,
we did not find this pattern for later growth, phenology-related traits or functional traits,
and we propose below different non-exclusive explanations for this. Moreover, we did not
find higher genetic variation in populations located in temporally or spatially variable envi-
ronments for any trait, which goes against predictions from theoretical or simulation work
(McDonald and Yeaman 2018, Walsh and Lynch 2018) and an empirical study in lodgepole
pine (Yeaman and Jarvis 2006). Admixed populations did not show higher quantitative
genetic variation, suggesting that the observed patterns are not confounded by historical
gene flow among distinct gene pools. Empirically-based detection of the potential footprints
of natural selection on within-population quantitative genetic variation is much needed to
understand how populations are adapted to their current environments and will evolve under
changing conditions.

## 4.1 Lower genetic variation for early height growth under colder winters

Genetic variation for early height growth was lower in populations experiencing colder winters (Figs. 2 and S52). This result is in line with the hypothesis that strong stabilizing selection in harsh environments depletes quantitative genetic variation within populations (Fisher 1930) and echoes similar results in another forest tree, *Q. oleoides*. For this Mesoamerican white oak species, Ramírez-Valiente et al. (2019) found lower genetic variation averaged over functional and growth traits in populations experiencing low precipitation and high temperatures during the dry season. The importance of cold temperatures as a driver of height genetic variation in maritime pine is supported by the association between candidate-gene allele frequency and temperature gradients in other studies (Grivet et al. 2011, Jaramillo-Correa et al. 2015), suggesting a major role of minimum temperatures in the species adaptive evolution. Lower quantitative genetic variation in environments with cold winters may therefore reflect
adaptation to local conditions, but it may also alter the ability of populations to respond to the new selection pressures induced by climate change (Lande and Shannon 1996, Storfer 1996) and decrease the fitness of individuals through the expression of inbreeding depression (but see Carley et al. 2022).

Mean coldest month temperatures were highly correlated with altitude in our study (Pearson's correlation of 0.79), and adaptation patterns along altitudinal gradients are common in forest trees (e.g. Kurt et al. 2012). Therefore, we cannot exclude that the association between genetic variation for early height growth and coldest month temperatures is triggered by more complex environmental factors typical of high altitude conditions (e.g. reduced vapor pressure deficit, higher maximum solar radiation; Körner 1995).

Differences in genetic variation for early height growth among populations were unlikely to originate from the expression of hidden genetic variation in novel environments (i.e. cryptic genetic variation caused by genotype-by-environment interactions; Gibson and Dworkin 2004, Schlichting 2008). First, the lower genetic variation for early height growth in populations experiencing colder winters was consistent across the three common gardens (i.e. independent of their environmental conditions), which was confirmed by the low genotype-by-environment interaction estimated in a complementary analysis (see section 8 of the Supplementary Information; Archambeau et al. 2022). Second, populations with higher quantitative genetic variation for early height growth were not those with larger differences in winter temperatures between their locations of origin and the common garden (i.e. higher climatic transfer distances; Fig. S50). In contrast, the expression of cryptic genetic variation may explain why populations with larger climatic transfer distances for winter temperatures showed higher quantitative genetic variation for SLA in Fundão (see section 6 of the Supplementary Information). From an evolutionary perspective, this pattern requires attention as cryptic genetic variation can enhance the evolutionary potential of populations by generating new phenotypes that are then available for natural selection (Ledón-Rettig et al. 2014), and may ultimately get fixed in the population, a process called 'genetic assimilation' (Waddington 1953).

Most importantly, the validation analysis in a family-based progeny test provided further evidence that additive within-population genetic variation for height was lower in populations experiencing colder winters for young trees but not for older trees (Fig. 3). This supports the robustness of our study and suggests that our results were unlikely to be biased by considering the total variance instead of the additive one, which was somehow expected as two previous studies in maritime pine found low non-additive effects for growth (Gaspar et al. 2013), and height and diameter (Lepoittevin et al. 2011).

## 4.2 Potential effects of demographic processes

Historical demographic processes, by generating population structure or changes in effective population size independently from natural selection, can also shape patterns of genetic variation (Lawton-Rauh 2008). Theory predicts that small and isolated populations are
likely to lose genetic variation more rapidly than large populations, due to genetic drift and limited gene flow from other populations. In this line, empirical studies have provided evidence that genetic diversity (measured with DNA markers) increases with population sizes (Frankham 1996, Leimu et al. 2006), although the magnitude of differences in population genetic diversity is marginal relative to the magnitude of changes in population size, a paradox which remains unresolved (Buffalo 2021). The relationship between quantitative genetic variation and population size is less clear. In their review, Willi et al. (2006) found that reduced quantitative genetic variation in small populations was supported in laboratory experiments in which populations evolve at different sizes, whereas most field studies do not support this prediction (see also Wood et al. 2015). In addition, individual-based simulations revealed that the additive genetic variance was less altered by a 10-fold change in population size than a reduction in migration (McDonald and Yeaman 2018).

Estimating effective population sizes of large and widespread forest tree populations is par-
ticularly challenging due to extensive sampling requirements (i.e. sufficient spatial coverage,
sample sizes, and genotyping intensity), and thus current estimates are probably unreliable
(Santos-del-Blanco et al. 2022). Therefore, although population sizes are normally large in
maritime pine, we can not rule out the hypothesis that low effective population sizes in harsh
environments, due to stochastic demographic events or recent colonization history, may have
caused the lower genetic variation observed for early height growth in our study. However,
several observations makes this hypothesis unlikely. First, reduced genetic variation was only
observed for early height growth and not for the other traits, which is more likely to occur un-
der the hypothesis that the underlying process is natural selection, rather than demographic
factors (under this second hypothesis, all traits would have the same probability of experi-
encing a decrease in variance). Second, we did not find any association between quantitative
genetic variation and genetic diversity estimated with molecular markers (i.e. expected het-
erozygosity) for any trait (see section 7 of the Supplementary Information; Rodríguez-Quilón
et al. 2015). Third, colonization history is unlikely to explain the reduced genetic variation
for early height growth in populations experiencing the coldest winter conditions, as these
populations are located in central and southeastern Spain, which are probably among the
first regions recolonized after the last glaciation (Bucci et al. 2007).

Populations in harsh environments may also be more isolated, and therefore less likely to benefit from increased quantitative genetic variance as a result a gene flow from nearby populations (as shown for instance in Reid and Arcese 2020). In our study, historical gene flow among gene pools most likely had a negligible impact on the expression of quantitative genetic variation in maritime pine populations, as the most admixed populations did not show higher quantitative genetic variation for any trait (even when the divergence between source and sink gene pools was accounted for). However, to our knowledge, evaluating the effect of contemporary gene flow on quantitative genetic variation using quantitative genetic models (such as the one used in our study) was until recently only possible using pedigree data (but see Aase et al. 2022), and therefore we were not able to rule out the hypothesis that recent exogenous gene flow in some populations may have induced an increase in their genetic variance independently of natural selection.
## 4.3 No association with spatial or temporal heterogeneity

Populations from spatially heterogeneous environments did not show higher genetic variation for any trait (Fig. 2), which was also the case in the validation analysis (Fig. 3). This goes against a previous estimate in lodgepole pine suggesting that up to 20% of the genetic variation in growth within populations is explained by spatial heterogeneity (Yeaman and Jarvis 2006). A potential explanation of this discrepancy is the smaller experiment size in our study compared to that of Yeaman and Jarvis (103 populations with an average of 28 planting sites per population). However, in our study, we obtained reasonable credible intervals for most traits (allowing the detection of associations with other drivers) and data simulations suggested that our models have adequate power, rendering this explanation unlikely. The discrepancy with Yeaman and Jarvis (2006) may also stem from the use of younger trees in our study. Indeed, the effect of environmental heterogeneity on the within-population quantitative genetic variation might be age-dependent, as might be the case also for the associations we detected with cold winters (see the section below for potential explanations of this phenomena).

Another explanation is that genetic variation within populations is not affected by the environmental heterogeneity at the regional scale imposed by the 1 × 1 km resolution of our climate dataset but at finer spatial scales (also discussed in Yeaman and Jarvis 2006). Indeed, populations can adapt along microgeographic environmental gradients despite the homogenizing effect of gene flow (Richardson et al. 2014), even for forest tree populations with their long-generation times and large effective population sizes (Scotti et al. 2016). However, a correlation between regional and microgeographic environmental heterogeneity across the maritime pine range is very likely: populations showing the highest environmental heterogeneity in our study were located in mountainous areas in which we also expect higher microgeographic variation, e.g. the Cúmpeta population (COM) located in the Tejeda and Almijara mountains (southern Spain), the Arenas de San Pedro population (ARN) located in the Sierra de Gredos (central Spain) or the Pineta population (PIE) located close to the Punta di Forchelli (Corsican mountains), while populations with the lowest environmental heterogeneity were located on flat plateaus, e.g. populations from the Landes plateau and the Atlantic coastal regions in France (HOU, MIM, PET, VER, OLO, STJ, PLE), and populations from the central Spain plateau near to Segovia (CUE, COC, CAR). Thus, even if genetic variation was maintained by migration-selection balance at microgeographic scales, we would have been able to detect the effect of environmental heterogeneity at the regional scale. Last, genetic variation within populations might also be affected by environmental heterogeneity at even greater scales than the one considered in our study, since, for instance, some potential environmental drivers of adaptation (e.g. wind, precipitation) vary at large scales. However, for these factors to affect genetic variation at the local scale substantial gene flow among distant populations would be needed, which does not seem the case in maritime pine, a species characterised by strong genetic differentiation of gene pools.

Finally, a last explanation is related to the different biological features between lodgepole pine and maritime pine. Lodgepole pine has extensive gene flow and low population structure ($F_{ST}$ = 0.016 in Yeaman et al. 2016) while maritime pine shows restricted gene flow with
strong population structure (at least six distinct gene pools and F<sub>ST</sub> = 0.112; Jaramillo-Correa et al. 2015; our study) and fragmented distribution (Alberto et al. 2013). Pollen dispersal kernels in maritime pine are highly leptokurtic, as for other wind-pollinated pines (Schuster and Mitton 2000, Robledo-Arnuncio and Gil 2005), with estimated mean dispersal distances from 78.4 to 174.4m (de-Lucas et al. 2008). Interestingly, McDonald and Yeaman (2018) showed that high levels of quantitative genetic variance can be maintained when a trait is under stabilizing selection only at intermediate levels of migration. Migration rates in maritime pine may therefore not be strong enough to compensate for the purifying effect of natural selection in spatially heterogeneous environments, especially in mountainous areas which may represent barriers to gene flow and where populations are more isolated (see González-Martínez et al. 2007 for maritime pine). Meanwhile, in the homogeneous plateaus of the Landes forest and central Spain, natural selection may be low because conditions are more favorable, and these populations are less isolated, which may maintain genetic variation at levels similar to those of populations in heterogeneous landscapes. Investigating local adaptation and gene flow at microgeographic scales in natural populations of maritime pine located in both homogeneous and heterogeneous environments would be highly valuable to understand why spatial environmental heterogeneity does not seem to play a major role in maintaining genetic variation in this species. Moreover, conducting similar analyses in sister species such as Scots pine, with low population genetic structure and continuous populations (Alberto et al. 2013), could help to determine whether genetic variation in forest tree populations experiencing higher migration rates are more prone to be impacted by spatial environmental heterogeneity.

The lack of robust support for an association between within-population quantitative genetic variation and temporal climatic heterogeneity was more expected than for spatial heterogeneity. Indeed, theory predicts lower levels of genetic variation for fitness under temporal heterogeneity than under spatial heterogeneity, though still higher than under constant environments (Levene 1953, Felsenstein 1976). These predictions were supported by empirical studies (e.g. Huang et al. 2015; but see Mackay 1981) and individual-based simulations showing that, in contrast to spatially heterogeneous environments, temporally heterogeneous environments can only marginally increase the additive genetic variance within populations (McDonald and Yeaman 2018). In forest trees, Ramírez-Valiente et al. (2019) found no association between quantitative genetic variation measured in 1-year-old seedlings of *Q. oleoides* and temporal environmental heterogeneity. One possible explanation is that long-lived, sessile species, such as forest trees, respond to temporal heterogeneity primarily through phenotypic plasticity, resulting in the selection of the most plastic genotypes in temporally fluctuating environments rather than the selection of a wide range of genotypes adapted to different conditions.

## 4.4 Potential causes underlying contrasting patterns across traits

The consistent association between within-population genetic variation and cold winter temperatures for early height growth was in line with the initial hypothesis that the effect of
natural selection on quantitative genetic variation is likely to be stronger (and thus easier to
detect) for traits more closely related to fitness. Indeed, height growth is an integrative trait
that can be seen as the end-product of multiple ecophysiological processes both genetically
and environmentally determined (Grattapaglia et al. 2009). Taller trees compete more effi-
ciently for light, water and nutrients (Morgenstern 1996), and are also more likely to have
high fecundity (Rehfeldt et al. 1999, Wu and Ying 2004, Aitken and Bemmels 2015). As
such, Younginger et al. (2017) reviewed 170 studies in herbs, shrubs and trees and showed
that, for conspecifics of the same age class, larger plants have generally a higher fitness,
i.e. they have more reproductive output and are therefore more likely to produce viable
offspring. In addition, as maritime pine is a pioneer species, selection pressures might be
particularly high in the early growth stages and might be exacerbated by increased com-
petition under harsh conditions (even at the expense of future tree longevity; Bigler and
Veblen 2009). Interestingly, no association between within-population genetic variation and
climate harshness was found for later height growth. It may result from the weaker selection
pressures for later growth as competition for light and nutrients is probably less intense after
the removal of most seedlings in the early selection phase.

Furthermore, the different patterns obtained for early and later height growth highlight that ontogenic stage is a key parameter to consider when studying the evolutionary forces driving the quantitative genetic variation in forest tree populations. In this vein, genetic parameters in forest trees have been shown to vary with age; e.g. heritability generally increases with age until reaching a plateau, especially for height-related traits (e.g. Johnson et al. 1997, Kroon et al. 2011). In maritime pine, an increase in heritability with age was found in Costa and Durel (1996) but not in Kusnandar et al. (1998). Interestingly, the expression of age-dependent heritabilities for another growth trait, diameter growth, was associated with spacing among trees in common gardens of *Pinus radiata* (Lin et al. 2013).

Leaf phenology-related traits exhibit steep adaptation gradients in forest trees and have a relatively high heritability, e.g. 0.15-0.51 for bud burst in pedunculate oak (Scotti-Saintagne et al. 2004), ~0.36 for bud burst and ~0.24 for the duration of bud burst in our study (Table S20). Gauzere et al. (2020) showed that both the mean and the variance of leaf phenology-related traits varied along an altitudinal gradient in natural oak populations, with populations at high altitude having a narrower fitness peak. We might therefore have expected lower genetic variation for leaf phenology-related traits in populations experiencing cold winters (and at higher altitude; see also Alberto et al. 2011). A biological explanation for this could be the weaker link with fitness in young trees for phenology-related traits when compared to early height growth, implying that they are less prone to a decrease in their genetic variance caused by strong selection pressures. Another explanation is that variation in phenology-related traits is likely to be controlled by genotype-by-environment interactions in forest trees (e.g. opposite genetic clines in common gardens and natural populations; Vitasse et al. 2009) and therefore such association may be hidden in the unique environment of one common garden. Last, since phenology-related traits are derived from discrete measurements, their distribution deviates more from a normal distribution than that of height, requiring more statistical power to detect associations.
We might also have expected lower genetic variation in populations experiencing drier summer conditions for traits considered to be related to drought resistance such as SLA (Mitchell et al. 2008; but see Reynolds et al. 2018) and $\delta^{13}C$ (Warren et al. 2001). In line with our study, Anderegg et al. (2021) found no differences in phenotypic variances for traits related to leaf and stem robustness and allocation along an aridity gradient in eucalypts. Investigating how the genetic variation of drought-resistance traits is impacted by intense drought conditions is worth addressing in face of climate change, as a lack of association may reflect a limited adaptive potential for these traits. However, genotype-by-environment interactions but also insufficient statistical power is likely to explain the lack of detected associations for these traits in our study since only a small proportion of their variation is genetically determined (Corcuera et al. 2010, Alía et al. 2014, $H_{SLA}^2 \sim 0.09$ and $H_{\delta^{13}C}^2 \sim 0.09$ in our study; Table S20), resulting in smaller effect sizes that are more difficult to detect.

Importantly, theoretical work suggests that much of the genetic variation associated with a trait is likely maintained by pleiotropic effects, which are independent of the selection on that trait, implying that stabilizing selection can only act on a reduced number of independent dimensions in the trait space (Barton 1990, Walsh and Lynch 2018). As we used univariate models, we cannot exclude that the associations detected for early height growth originate from genetic correlation with other traits under selection, or that the lack of association with other traits (notably functional traits such as $\delta^{13}C$) does not originate from genetic constraints (Walsh and Blows 2009). For example, in maritime pine, trait canalisation and genetic constraints may explain low quantitative genetic differentiation for hydraulic traits (e.g. P50, the xylem pressure inducing 50% loss of hydraulic conductance; Lamy et al. 2014). In this line, Benavides et al. (2021) found higher trait covariance in Scots pine populations under harsher environments, thus suggesting that strong selection pressures may not only reduce the genetic variance of individual traits but may also generate more coordinated and efficient phenotypes ('tighter' phenotypes), a process that deserves further investigation.

# 5 Conclusion

Our manuscript contributes to the current debate on the maintenance of quantitative genetic variation within populations by providing empirical support for the role of natural selection in decreasing genetic variation for a key trait of early survival in forest trees. Indeed, our results consistently showed that genetic variation for early height growth is lower in maritime pine populations experiencing colder winters (i.e. experiencing stronger climate-related selection). Surprisingly, we did not find any association between spatial or temporal environmental heterogeneity and within-population genetic variation for any trait; whether for technical reasons (e.g. sample size, spatial scale considered) or for genuine biological reasons (e.g. too low migration, ontogenic stage), it would be worth further exploration. Importantly, this work contributes to understanding how evolutionary forces shape genetic variation in key adaptive traits within populations and ultimately how populations adapt to their local environment. However, whether the reduction in quantitative genetic varia-
tion potentially induced by cold temperatures has significant evolutionary implications for the adaptive capacity of populations remains unknown and would require more prediction-oriented studies.

## References

Aase K, Jensen H and Muff S (2022). Genomic estimation of quantitative genetic parameters in wild admixed populations. Methods in Ecology and Evolution 13:1014–1026.

Aitken SN and Bemmels JB (2015). Time to get moving: assisted gene flow of forest trees. Evolutionary Applications 9:271–290.

Alberto F, Aitken SN, Alía R, González-Martínez SC, Hänninen H, Kremer A et al. (2013). Potential for evolutionary responses to climate change – evidence from tree populations. Global Change Biology 19:1645–1661.

Alberto F, Bouffier L, Louvet JM, Lamy JB, Delzon S and Kremer A (2011). Adaptive responses for seed and leaf phenology in natural populations of sessile oak along an altitudinal gradient. Journal of Evolutionary Biology 24:1442–1454.

Alía R, Chambel R, Notivol E, Climent J and González-Martínez SC (2014). Environment-dependent microevolution in a Mediterranean pine (Pinus pinaster Aiton). BMC Evolutionary Biology 14:1–12.

Anderegg LDL, Loy X, Markham IP, Elmer CM, Hovenden MJ, HilleRisLambers J et al. (2021). Aridity drives coordinated trait shifts but not decreased trait variance across the geographic range of eight Australian trees. New Phytologist 229:1375–1387.

Archambeau J, Garzón MB, Barraquand F, Miguel M de, Plomion C and González-Martínez SC (2022). Combining climatic and genomic data improves range-wide tree height growth prediction in a forest tree. The American Naturalist 200:E141–E159.

Barton NH (1990). Pleiotropic models of quantitative variation. Genetics 124:773–782.

Benavides R, Carvalho B, Matesanz S, Bastias CC, Cavers S, Escudero A et al. (2021). Phenotypes of Pinus sylvestris are more coordinated under local harsher conditions across Europe. Journal of Ecology 109:2580–2596.

Benito Garzón M, Robson TM and Hampe A (2019). ΔTraitSDMs: species distribution models that account for local adaptation and phenotypic plasticity. New Phytologist 222:1757–1765.

Bigler C and Veblen TT (2009). Increased early growth rates decrease longevities of conifers in subalpine forests. Oikos 118:1130–1138.

Bucci G, González-Martínez SC, Provost GL, Plomion C, Ribeiro MM, Sebastiani F et al. (2007). Range-wide phylogeography and gene zones in Pinus pinaster Ait. revealed by chloroplast microsatellite markers. Molecular Ecology 16:2137–2153.

Buffalo V (2021). Quantifying the relationship between genetic diversity and population size suggests natural selection cannot explain Lewontin's Paradox. eLife 10:e67509.
Byers DL (2005). Evolution in heterogeneous environments and the potential of maintenance of genetic variation in traits of adaptive significance. *Genetics of Adaptation:107–124.

Carley LN, Morris WF, Walsh R, Riebe D and Mitchell-Olds T (2022). Are genetic variation and demographic performance linked? *Evolutionary Applications* 15:1888–1906.

Carpenter B, Gelman A, Hoffman MD, Lee D, Goodrich B, Betancourt M et al. (2017). *Stan: A Probabilistic Programming Language*. Journal of Statistical Software 76.

Charmantier A, Kruuk LEB, Blondel J and Lambrechts MM (2004). Testing for microevolution in body size in three blue tit populations. *Journal of Evolutionary Biology* 17:732–743.

Colautti RI, Eckert CG and Barrett SCH (2010). Evolutionary constraints on adaptive evolution during range expansion in an invasive plant. *Proceedings of the Royal Society B: Biological Sciences* 277:1799–1806.

Conrad O, Bechtel B, Bock M, Dietrich H, Fischer E, Gerlitz L et al. (2015). System for Automated Geoscientific Analyses (SAGA) v. 2.1.4. *Geoscientific Model Development* 8:1991–2007.

Corcuera L, Gil-Pelegrin E and Notivol E (2010). Phenotypic plasticity in Pinus pinaster *δ*13C: environment modulates genetic variation. *Annals of Forest Science* 67:812–812.

Costa P and Durel CE (1996). Time trends in genetic control over height and diameter in maritime pine. *Canadian Journal of Forest Research* 26:1209–1217.

De-Lucas AI, Robledo-Arnuncio JJ, Hidalgo E and González-Martínez SC (2008). Mating system and pollen gene flow in Mediterranean maritime pine. *Heredity* 100:390–399.

De Miguel M, Rodríguez-Quilón I, Heuertz M, Hurel A, Grivet D, Jaramillo-Correa JP et al. (2022). Polygenic adaptation and negative selection across traits, years and environments in a long-lived plant species (Pinus pinaster Ait., Pinaceae). *Molecular Ecology*.

Falush D, Stephens M and Pritchard JK (2003). Inference of population structure using multilocus genotype data: linked loci and correlated allele frequencies. *Genetics* 164:1567–1587.

Felsenstein J (1976). The theoretical population genetics of variable selection and migration. *Annual Review of Genetics* 10:253–280.

Fisher RA (1930). *The genetical theory of natural selection*. Oxford: Clarendon Press. 289 pp.

Frankham R (1996). Relationship of genetic variation to population size in wildlife. *Conservation Biology* 10:1500–1508.

Gaspar MJ, Velasco T, Feito I, Alía R and Majada J (2013). Genetic variation of drought tolerance in Pinus pinaster at three hierarchical levels: a comparison of induced osmotic stress and field testing. *PLOS ONE* 8:e79094.

Gauzere J, Teuf B, Davi H, Chevin LM, Caignard T, Leys B et al. (2020). Where is the optimum? Predicting the variation of selection along climatic gradients and the adaptive value of plasticity. A case study on tree phenology. *Evolution Letters* 4:109–123.

Gibson G and Dworkin I (2004). Uncovering cryptic genetic variation. *Nature Reviews Genetics* 5:681–690.
González-Martínez SC, Alía R and Gil L (2002). Population genetic structure in a Mediterranean pine (Pinus pinaster Ait.): a comparison of allozyme markers and quantitative traits. Heredity 89:199-206.

González-Martínez SC, Gómez A, Carrión JS, Agúndez D, Alía R and Gil L (2007). Spatial genetic structure of an explicit glacial refugium of maritime pine (Pinus pinaster Aiton) insoutheastern Spain. *Phylogeography of Southern European Refugia: Evolutionary Perspectives on the Origins and Conservation of European Biodiversity*. Ed. by Weiss S and Ferrand N. Dordrecht: Springer Netherlands. 257-269.

Grattapaglia D, Plomion C, Kirst M and Sederoff RR (2009). Genomics of growth traits in forest trees. Current Opinion in Plant Biology 12:148-156.

Grivet D, Sebastiani F, Alía R, Bataillon T, Torre S, Zabal-Aguirre M et al. (2011). Molecular footprints of local adaptation in two mediterranean conifers. Molecular Biology and Evolution 28:101-116.

Hamrick JL (2004). Response of forest trees to global environmental changes. Forest Ecology and Management 197:323-335.

Hedrick PW (1986). Genetic polymorphism in heterogeneous environments: a decade later. Annual Review of Ecology and Systematics 17:535-566.

(— 2006). Genetic polymorphism in heterogeneous environments: the age of genomics. Annual Review of Ecology, Evolution, and Systematics 37:67-93.

Hereford J, Hansen TF and Houle D (2004). Comparing strengths of directional selection: how strong is strong? Evolution 58:2133-2143.

Hiederer R, European Commission, Joint Research Centre and Institute for Environment and Sustainability (2013). *Mapping soil properties for Europe spatial representation of soil database attributes*. Luxembourg: Publications Office of the European Union.

Hoffmann AA and Parsons PA (1991). *Evolutionary genetics and environmental stress*. Oxford University Press.

Huang Y, Stinchcombe JR and Agrawal AF (2015). Quantitative genetic variance in experimental fly populations evolving with or without environmental heterogeneity. Evolution 69:2735-2746.

Hurel A, de Miguel M, Dutech C, Desprez-Loustau ML, Plomion C, Rodríguez-Quilón I et al. (2021). Genetic basis of growth, spring phenology, and susceptibility to biotic stressors in maritime pine. Evolutionary Applications 14:2750-2772.

Jaramillo-Correa JP, Rodríguez-Quilón I, Grivet D, Lepoittevin C, Sebastiani F, Heuertz M et al. (2015). Molecular proxies for climate maladaptation in a long-lived tree (Pinus pinaster Aiton, Pinaceae). Genetics 199:793-807.

Johnson GR, Sniezko RA and Mandel NL (1997). Age trends in Douglas-fir genetic parameters and implications for optimum selection age. Silvae Genetica. 46: 349-358.

Johnson T and Barton NH (2005). Theoretical models of selection and mutation on quantitative traits. Philosophical Transactions of the Royal Society B: Biological Sciences 360:1411-1425.

Jump AS and Peñuelas J (2005). Running to stand still: adaptation and the response of plants to rapid climate change. Ecology Letters 8:1010-1020.

Körner C (1995). Alpine plant diversity: a global survey and functional interpretations. *Arctic and Alpine Biodiversity: Patterns, Causes and Ecosystem Consequences*. Ed. by Chapin FS and Körner C. Berlin, Heidelberg: Springer. 45-62.
Kremer A, Ronce O, Robledo-Arnuncio JJ, Guillaume F, Bohrer G, Nathan R et al. (2012). Long-distance gene flow and adaptation of forest trees to rapid climate change. Ecology Letters 15:378-392.

Kroon J, Ericsson T, Jansson G and Andersson B (2011). Patterns of genetic parameters for height in field genetic tests of Picea abies and Pinus sylvestris in Sweden. Tree Genetics & Genomes 7:1099-1111.

Kurt Y, González-Martínez SC, Alía R and Isik K (2012). Genetic differentiation in Pinus brutia Ten. using molecular markers and quantitative traits: the role of altitude. Annals of Forest Science 69:345-351.

Kusnandar D, Galwey NW, Hertzler GL and Butcher TB (1998). Age trends in variances and heritabilities for diameter and height in maritime pine (Pinus pinaster Ait.) in western Australia. Silvae Genetica 47:136-141.

Lamy JB, Delzon S, Bouche PS, Alia R, Vendramin GG, Cochard H et al. (2014). Limited genetic variability and phenotypic plasticity detected for cavitation resistance in a Mediterranean pine. New Phytologist 201:874-886.

Lande R and Shannon S (1996). The role of genetic variation in adaptation and population persistence in a changing environment. Evolution 50:434-437.

Lawton-Rauh A (2008). Demographic processes shaping genetic variation. Current Opinion in Plant Biology 11:103-109.

Ledón-Rettig CC, Pfennig DW, Chunco AJ and Dworkin I (2014). Cryptic genetic variation in natural populations: a predictive framework. Integrative and Comparative Biology 54:783-793.

Leimu R, Mutikainen P, Koricheva J and Fischer M (2006). How general are positive relationships between plant population size, fitness and genetic variation? Journal of Ecology 94:942-952.

Leites LP, Robinson AP, Rehfeldt GE, Marshall JD and Crookston NL (2012). Height-growth response to climatic changes differs among populations of Douglas-fir: a novel analysis of historic data. Ecological Applications 22:154-165.

Lemoine NP (2019). Moving beyond noninformative priors: why and how to choose weakly informative priors in Bayesian analyses. Oikos 128:912-928.

Lepoittevin C, Rousseau JP, Guillemin A, Gauvrit C, Besson F, Hubert F et al. (2011). Genetic parameters of growth, straightness and wood chemistry traits in Pinus pinaster. Annals of Forest Science 68:873-884.

Lerner IM (1958). The genetic basis of selection. *The genetic basis of selection*. Chapman & Hall, Ltd., London & John Wiley & Sons, Inc., New York.

Levene H (1953). Genetic equilibrium when more than one ecological niche is available. The American Naturalist 87:331-333.

Lin Y, Yang H, Ivković M, Gapare WJ, Colin Matheson A and Wu HX (2013). Effect of genotype by spacing interaction on radiata pine genetic parameters for height and diameter growth. Forest Ecology and Management 304:204-211.

Linhart YB and Grant MC (1996). Evolutionary Significance of Local Genetic Differentiation in Plants. Annual Review of Ecology and Systematics 27:237-277.

Mackay TFC (1981). Genetic variation in varying environments. Genetics Research 37:79-93.
Marchi M, Castellanos-Acuña D, Hamann A, Wang T, Ray D and Menzel A (2020). ClimateEU, scale-free climate normals, historical time series, and future projections for Europe. Scientific Data 7:428.

McDonald TK and Yeaman S (2018). Effect of migration and environmental heterogeneity on the maintenance of quantitative genetic variation: a simulation study. Journal of Evolutionary Biology 31:1386-1399.

Merilä J, Sheldon B and Kruuk L (2001). Explaining stasis: microevolutionary studies in natural populations. Genetica 112:199-222.

Merilä J, Söderman F, O'Hara R, Räsänen K and Laurila A (2004). Local adaptation and genetics of acid-stress tolerance in the moor frog, Rana arvalis. Conservation Genetics 5:513-527.

Mitchell PJ, Veneklaas EJ, Lambers H and Burgess SSO (2008). Leaf water relations during summer water deficit: differential responses in turgor maintenance and variation in leaf structure among different plant communities in south-western Australia. Plant, Cell & Environment 31:1791-1802.

Mitchell-Olds T, Willis JH and Goldstein DB (2007). Which evolutionary processes influence natural genetic variation for phenotypic traits? Nature Reviews Genetics 8:845-856.

Morgenstern M (1996). *Geographic variation in forest trees: genetic basis and application of knowledge in silviculture*. Vancouver: Univ. B.C. Press. 224 pp.

Neale DB and Savolainen O (2004). Association genetics of complex traits in conifers. Trends in Plant Science 9:325-330.

Petit RJ and Hampe A (2006). Some evolutionary consequences of being a tree. Annual Review of Ecology, Evolution, and Systematics 37:187-214.

Plomion C, Bartholomé J, Lesur I, Boury C, Rodríguez-Quilón I, Lagraulet H et al. (2016). High-density SNP assay development for genetic analysis in maritime pine (Pinus pinaster). Molecular Ecology Resources 16:574-587.

Pritchard JK, Stephens M and Donnelly P (2000). Inference of population structure using multilocus genotype data. Genetics 155:945-959.

R Core Team (2020). *R: A Language and Environment for Statistical Computing*. Vienna, Austria: R Foundation for Statistical Computing.

Ramírez-Valiente JA, Etterson JR, Deacon NJ and Cavender-Bares J (2019). Evolutionary potential varies across populations and traits in the neotropical oak Quercus oleoides. Tree Physiology 39:427-439.

Rehfeldt GE, Ying CC, Spittlehouse DL and Hamilton DA (1999). Genetic responses to climate in Pinus contorta: niche breadth, climate change, and reforestation. Ecological Monographs 69:375-407.

Reid JM and Arcese P (2020). Recent immigrants alter the quantitative genetic architecture of paternity in song sparrows. Evolution Letters 4:124-136.

Reynolds VA, Anderegg LDL, Loy X, HilleRisLambers J and Mayfield MM (2018). Unexpected drought resistance strategies in seedlings of four Brachychiton species. Tree Physiology 38:664-677.

Richardson JL, Urban MC, Bolnick DI and Skelly DK (2014). Microgeographic adaptation and the spatial scale of evolution. Trends in Ecology & Evolution 29:165-176.
Riley SJ, DeGloria SD and Elliot R (1999). Index that quantifies topographic heterogeneity. Intermountain Journal of Sciences 5:23-27.

Robledo-Arnuncio JJ and Gil L (2005). Patterns of pollen dispersal in a small population of Pinus sylvestris L. revealed by total-exclusion paternity analysis. Heredity 94:13-22.

Rodríguez-Quilón I, Santos-del-Blanco L, Grivet D, Jaramillo-Correa JP, Majada J, Vendramin GG et al. (2015). Local effects drive heterozygosity-fitness correlations in an outcrossing long-lived tree. Proceedings of the Royal Society B: Biological Sciences 282:20152230.

Rodríguez-Quilón I, Santos-del-Blanco L, Serra-Varela MJ, Koskela J, González-Martínez SC and Alía R (2016). Capturing neutral and adaptive genetic diversity for conservation in a highly structured tree species. Ecological Applications 26:2254-2266.

Santos-del-Blanco L, Olsson S, Budde KB, Grivet D, González-Martínez SC, Alía R et al. (2022). On the feasibility of estimating contemporary effective population size (Ne) for genetic conservation and monitoring of forest trees. Biological Conservation 273:109704.

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

Schlichting CD (2008). Hidden reaction norms, cryptic genetic variation, and evolvability. Annals of the New York Academy of Sciences 1133:187-203.

Schuster WSF and Mitton JB (2000). Paternity and gene dispersal in limber pine (Pinus flexilis James). Heredity 84:348-361.

Scotti I, González-Martínez SC, Budde KB and Lalagüe H (2016). Fifty years of genetic studies: what to make of the large amounts of variation found within populations? Annals of Forest Science 73:69-75.

Scotti-Saintagne C, Bodénès C, Barreneche T, Bertocchi E, Plomion C and Kremer A (2004). Detection of quantitative trait loci controlling bud burst and height growth in Quercus robur L. Theoretical and Applied Genetics 109:1648-1659.

Spichtig M and Kawecki TJ (2004). The maintenance (or not) of polygenic variation by soft selection in heterogeneous environments. The American Naturalist 164:70-84.

Stein A, Gerstner K and Kreft H (2014). Environmental heterogeneity as a universal driver of species richness across taxa, biomes and spatial scales. Ecology Letters 17:866-880.

Stock AJ, Campitelli BE and Stinchcombe JR (2014). Quantitative genetic variance and multivariate clines in the Ivyleaf morning glory, Ipomoea hederacea. Philosophical Transactions of the Royal Society B: Biological Sciences 369:20130259.

Storfer A (1996). Quantitative genetics: a promising approach for the assessment of genetic variation in endangered species. Trends in Ecology & Evolution 11:343-348.

Thibault E, Soolanayakanahally R and Keller SR (2020). Latitudinal clines in bud flush phenology reflect genetic variation in chilling requirements in balsam poplar, Populus balsamifera. American Journal of Botany 107:1597-1605.

Triantis KA, Mylonas M, Lika K and Vardinoyannis K (2003). A model for the species-area-habitat relationship. Journal of Biogeography 30:19-27.

Van Heerwaarden B, Kellermann V, Schiffer M, Blacket M, Sgrò CM and Hoffmann AA (2009). Testing evolutionary hypotheses about species borders: patterns of genetic variation towards the southern borders of two rainforest Drosophila and a related
habitat generalist. Proceedings of the Royal Society B: Biological Sciences 276:1517-1526.

Viñas RA, Caudullo G, Oliveira S and de Rigo D (2016). Pinus pinaster in Europe: distribution, habitat, usage and threats.

Vitasse YV, Delzon SD, Bresson CCBC, Michalet RM and Kremer AK (2009). Altitudinal differentiation in growth and phenology among populations of temperate-zone tree species growing in a common garden. Canadian Journal of Forest Research.

Waddington CH (1953). Genetic assimilation of an acquired character. Evolution 7:118-126.

Walsh B and Blows MW (2009). Abundant genetic variation + strong selection = multivariate genetic constraints: a geometric view of adaptation. Annual Review of Ecology, Evolution, and Systematics 40:41-59.

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

Warren CR, McGrath JF and Adams MA (2001). Water availability and carbon isotope discrimination in conifers. Oecologia 127:476-486.

Willi Y, van Buskirk J and Hoffmann AA (2006). Limits to the adaptive potential of small populations. Annual Review of Ecology, Evolution, and Systematics 37:433-458.

Wood JLA, Tezel D, Joyal D and Fraser DJ (2015). Population size is weakly related to quantitative genetic variation and trait differentiation in a stream fish. Evolution 69:2303-2318.

Wu HX and Ying CC (2004). Geographic pattern of local optimality in natural populations of lodgepole pine. Forest Ecology and Management 194:177-198.

Yeaman S, Chen Y and Whitlock MC (2010). No effect of environmental heterogeneity on the maintenance of genetic variation in wing shape in Drosophila melanogaster. Evolution 64:3398-3408.

Yeaman S, Hodgins KA, Lotterhos KE, Suren H, Nadeau S, Degner JC et al. (2016). Convergent local adaptation to climate in distantly related conifers. Science 353:1431-1433.

Yeaman S and Jarvis A (2006). Regional heterogeneity and gene flow maintain variance in a quantitative trait within populations of lodgepole pine. Proceedings of the Royal Society B: Biological Sciences 273:1587-1593.

Younginger BS, Sirová D, Cruzan MB and Ballhorn DJ (2017). Is biomass a reliable estimate of plant fitness? Applications in Plant Sciences 5:1600094.