Harnessing tree-ring phenotypes to disentangle gene by environment interactions and their climate dependencies in a circum-Mediterranean pine

Abstract

Background and Aims

Understanding the genetic basis of adaptation and plasticity in trees constitutes a knowledge gap. We linked dendrochronology and genomics [single nucleotide polymorphisms (SNPs)] for a widespread conifer (Pinus halepensis Mill.) to characterize intraspecific growth differences elicited by climate.

Methods

The analysis comprised 20-year tree-ring series of 130 trees structured in 23 populations evaluated in a common garden. We tested for genotype by environment interactions (G × E) of indexed ring width (RWI) and early- to latewood ratios (ELI) using factorial regression, which describes G × E as differential gene sensitivity to climate.

Key Results

The species’ annual growth was positively influenced by winter temperature and spring moisture and negatively influenced by previous autumn precipitation and warm springs. Four and five climate factors explained 10 % (RWI) and 16 % (ELI) of population-specific interannual variability, respectively, with populations from drought-prone areas and with uneven precipitation experiencing larger growth reductions during dry vegetative periods. Furthermore, four and two SNPs explained 14 % (RWI) and 10 % (ELI) of interannual variability among trees, respectively. Two SNPs played a putative role in adaptation to climate: one identified from transcriptome sequencing of P. halepensis and another involved in response regulation to environmental stressors.

Conclusions

We highlight how tree-ring phenotypes, obtained from a common garden experiment, combined with a candidate-gene approach allow the quantification of genetic and environmental effects determining adaptation for a conifer with a large and complex genome.

Adaptive variation, common garden, dendroecology, factorial regression, Pinus halepensis, phenotypic plasticity, single nucleotide polymorphisms

INTRODUCTION

Tree growth performance is limited in dry areas such as the Mediterranean region, where seasonal drought along with intensification in the length, frequency and severity of extreme climate events are limiting ecosystem functioning (Martín-Benito et al., 2010; Forner et al., 2018; Helluy et al., 2020). Changes in the population dynamics of tree species, which are partly linked to migration processes and genetic adaptedness, are therefore expected under future climate (Housset et al., 2018; Royer-Tardif et al., 2021). The pace of genetic adaptations to new conditions, however, is likely to be too slow to conveniently track global warming (Jezkova et al., 2016). At the same time, responses to environmental factors may differ strongly among populations. This is because trees are usually locally adapted because of inter- and intraspecific genetic differentiation; in turn, disparate selective pressures also trigger differential plasticity among individuals (Benito Garzón et al., 2011; Grivet et al., 2011). Thus, some populations may have higher adaptive potential than others when facing climate instability. Yet, the prediction of tree responses to future climate is constrained by knowledge gaps on the genetic basis of adaptation as well as on the prevailing climate drivers of growth and survival (Housset et al., 2018).

Pinus halepensis Mill. (Aleppo pine) is a very plastic, drought-resistant conifer able to adjust its growth rate during extended drought periods (Camarero et al., 2010; Pasho et al., 2012). It presents large intraspecific variability in life-history traits (Climent et al., 2008; Voltas et al., 2018; Santini et al., 2019), which is reflected in local adaptations to the different niches inhabited by this species (Patsiou et al., 2020; Hevia et al., 2020). Previous studies have demonstrated that the intraspecific variation in traits related to growth (Voltas et al., 2018; Patsiou et al., 2020), wood anatomy (Hevia et al., 2020) or reproduction (Climent et al., 2008) is differentially affected by changing environmental conditions (i.e. differential phenotypic plasticity). To identify the genetic basis of this differentiation, recent studies have targeted single nucleotide polymorphism (SNP) markers (Pinosio et al., 2014) potentially related to adaptive variation (Daniels et al., 2019; Santini et al., 2021). SNPs offer a straight link to gene functions (i.e. candidate gene approach) without the need to investigate the whole genome (Jaramillo-Correa et al., 2015).

Intraspecific variation in Aleppo pine is customarily examined through functional traits measured at a particular age (Klein et al., 2013; Santini et al., 2019). Trait expression can be therefore considered as the outcome of an individual’s response to the environment at a certain time (e.g. physiological attributes), or as the integrative response over an individual’s life measured in that particular moment (e.g. morphometric characters). As an alternative, investigating wood growth traits (i.e. secondary growth) in common gardens, in combination with a candidate-gene approach, constitutes an exceptional opportunity to characterize the importance of genetic and environmental effects determining tree performance dynamics on a multi-annual scale (Depardieu et al., 2021). Basic wood features such as tree-ring width are suitable to assess adaptation and plasticity during climate fluctuations (Montwé et al., 2016; Martínez-Sancho et al., 2021) because tree rings reflect climatic sensitivity across the lifespan of individuals (Housset et al., 2016). Additional information on tree growth sensitivity to climate can also be provided by alternative traits such as the relative proportion of early- to latewood within a ring as the formation of these wood components is related to the prevailing climate occurring during specific seasonal windows (Torbenson et al., 2016). The early- to latewood ratio is particularly informative with regard to the hydraulic and mechanical properties of the main trunk (Domec and Gartner, 2002; De Luis et al., 2011; Novak et al., 2013). A high earlywood to latewood proportion is related to a more conductive wood with wider xylem vessels, whereas a high latewood proportion implies a denser wood. The relative importance of these wood types is mainly determined by rainfall seasonality (i.e. spring to autumn precipitation) and it influences the trade-off among hydraulic conductivity, resistance to cavitation and mechanical stability of trees (De Luis et al., 2011; Camarero et al., 2021).

In this study, we analysed the interannual growth patterns of 23 Aleppo pine populations covering most of the species’ natural distribution range (Fig. 1) that were growing in a common garden experiment representative of the average climate conditions for the species. To explore the genetic basis of secondary growth determination, we attempted to link genes and phenotypes structured in populations with potential climate drivers of growth variability. For this, we tested for genotype by (multi-year) environment interactions (G × E) using ring widths and early- to latewood ratios as phenotypic traits. We used a factorial regression approach as a meaningful method for the analysis and interpretation of G × E (van Eeuwijk et al., 2002). Factorial regression allows the inclusion of explicit environmental and genetic information (i.e. covariates) in G × E models along with a direct evaluation (i.e. quantification) of the importance of these covariates for G × E explanation (Malosetti et al., 2013; Sixto et al., 2016). We tested models in which environmental covariables were represented by monthly climate factors and genetic predictors were defined by particular SNPs localized in known candidate genes, aiming at describing G × E through explicit (biophysical and genetic) underlying factors. Importantly, genetic markers allow for investigating the relevance of additive and non-additive genetic variation in trait expression. Thus, we evaluated different genetic effects (additive, dominant, first-order epistasis) on secondary growth by defining the type of action of candidate genes relevant for the explanation of differential expression across years (Malosetti et al., 2004; Vargas et al., 2006; Calleja-Rodríguez et al., 2021).

Fig. 1.

Geographical origin of 23 Aleppo pine populations (coloured dots) growing in a common garden (black dot) located in Altura (Castellón province, Spain). Population acronyms (described in Supplementary Data Table S1) are shown on top of each population origin. The dark green areas represent the natural distribution of Aleppo pine according to EUFORGEN (http://www.euforgen.org/species/pinus-halepensis/). An aerial picture of the trial is also included, where the red rectangle exemplifies an experimental unit consisting of four trees of the same population.

We expected a diverse set of interannual growth patterns among populations of Aleppo pine as a result of adaptive divergence, as previously described for many functional traits (e.g. Climent et al., 2008; Martín-Sanz et al., 2019; Lombardi et al., 2021), including secondary growth (Hevia et al., 2020). In particular, we expected populations from warmer and drier origins to be less affected by hot and dry growth periods (i.e. having wider rings and larger early- to latewood ratios) compared with their wetter and colder counterparts as the result of local adaptation, in line with previous results in other conifers (Depardieu et al., 2020). We also hypothesized that interannual divergence in secondary growth can be explained, at least in part, by some SNP loci previously linked to morpho-physiological variability of Aleppo pine in the same trial (Santini et al., 2021), thereby modulating adaptive differentiation elicited by climate. This study is novel in the sense that it models G × E interactions of tree growth for a widespread conifer by integrating information on climate and genomics (SNPs associated with candidate genes) at an interannual scale using tree-ring phenotypes. It benefits from dendrochronological series of secondary growth which constitute archived and readily available phenotypic information in forest trees.

MATERIAL AND METHODS

Genetic trial

The genetic trial (common garden) was established in 1997 and is located in the municipality of Altura (39°49ʹ29ʹʹN, 00°34ʹ22ʹʹW, 640 m a.s.l; Castellón Province, eastern Spain). The site is representative of the average climate condition of the species distribution range, with a mean annual temperature of 13.8 °C and mean annual precipitation of 468 mm (Patsiou et al., 2020; Lombardi et al., 2021). It consists of 896 individual adults of Aleppo pine belonging to 56 populations that cover most of the species’ natural distribution range across the Mediterranean basin. Seeds from these populations were collected in 1995 from ~25 trees spaced at least 100 m apart to minimize the kinship of individuals within populations. Seeds were sown in a forest nursery in Spain and seedlings were planted at the study site at 1-year-old at a spacing of 2.5 m following row and column directions. Each experimental unit comprised four individual trees of the same population planted across the same row. The trial was set up according to a Latinized row–column design with four replicates. Out of the 412 surviving trees in November 2019, 130 trees belonging to 23 populations representative of the natural habitat of the species were used because of a thinning treatment conducted in November 2019, where about half of the trees originally presented in the trial had been systematically cut.

Meteorological data

Values of monthly maximum, minimum and mean temperatures (T_max, T_min, T_mean, respectively) as well as monthly precipitation (P) were retrieved from the nearest grid point to the trial location of the gridded climate dataset (Climate Research Unit, CRU TS 4.04 data set, Harris et al., 2020) for 2000–2019. The CRU provides monthly climate series on a 0.5° × 0.5° grid-box basis, interpolated from meteorological stations across the globe. We used this interpolated dataset due to the lack of continuous local climate records in the area. In addition, the Standardized Precipitation Evapotranspiration Index (SPEI) was used as a measure of drought intensity and was calculated at 1- and 6-month scales (SPEI1 and SPEI6, respectively) using the package SPEI in R (Vicente-Serrano et al., 2010).

For each population origin, mean historical climate records were instead obtained from the higher spatial resolution WorldClim database (Fick and Hijmans, 2017). WorldClim provides averages for 1970–2000 of 19 bioclimatic variables derived from temperature and precipitation records at 30-arcsecond resolution (~1 km²). In particular, and based on previous studies of climatic drivers of ecotypic variation in Mediterranean pines (Tapias et al., 2004; Climent et al., 2008), we focused on the mean annual temperature (MAT), mean temperature of the warmest quarter (or summer temperature; MST), mean temperature of the coldest quarter (or winter temperature; MWT), maximum temperature of the warmest month (T_maxW), minimum temperature of the coldest month (T_minC), mean annual precipitation (MAP), precipitation seasonality (or the coefficient of variation in monthly precipitation over the year, PS), precipitation of the driest quarter (PDQ) and the summer to annual precipitation ratio (PsP).

Tree-ring records

For the subset of 130 individuals, we collected cross-sections of about 5 cm thick from the basal part of the trunk. The cross-sections were cut into planks of about 5 cm width, which comprised both ends of each slice, with the pith centred longitudinally. The planks were dried and sanded with progressively finer sandpaper until the surface was smooth enough and the rings were clearly visible.

Tree-ring width was measured with a precision of 0.01 mm through a semi-automatic process in WinDendro 2014a (Regent instrument Inc., Quebec, Canada), and manual corrections were made when necessary. Each tree-ring series was cross-dated and the cross-dating was quality-checked with COFECHA (Holmes, 1983). In addition, earlywood width (EW) and latewood width (LW) were estimated using WinDendro coupled with a computer-integrated Leica binocular microscope (5× zoom). The transition between EW and LW was visually discerned based on changes in wood colour within each ring, where EW corresponded to the light-coloured portion and LW to the darker coloured part of the ring (Griffin et al., 2011; Cabral-Alemán et al., 2017). The EW to LW ratio (EL) was then calculated for each tree ring. Ring-width and EL indices (RWI and ELI hereafter) were obtained for each tree-ring series using the detrend function of the R package dplR (Bunn, 2008; Bunn et al., 2021). To this end, we applied a cubic spline with a 50 % frequency cut-off of a wavelength of half the total number of years to remove non-climatic (i.e. ontogenic) trends from each tree-ring series. This procedure generated a stationary (mean = 1) series of dimensionless indices that preserved a common variance encompassing interannual time scales. Afterwards, the master chronology of both indices (i.e. across all available trees and populations) was built through the crn function (R package dplR) using Tukey’s biweight robust mean.

Genetic data

We used a subset of 20 SNPs derived from a dataset of 294 SNPs originally disclosed from comprehensive transcriptome analyses of P. halepensis and re-sequenced loci identified in P. taeda (Pinosio et al., 2014). The subset of 20 SNPs (Supplementary Data Table S2) was selected from a previous genome-wide association study (GWAS) which tested associations among morphological and multispectral-derived physiological traits and genotypes at single loci (Santini et al., 2021). The GWAS analysis was carried out taking into account the neutral genetic structure of populations, and thus the subset of SNPs used in our study was previously detected after correction for this structure. As a result, SNPs were found as being directly related to tree growth or indirectly related to (remotely sensed) leaf area, photosynthetic pigments and leaf water content in Aleppo pine (Santini et al., 2021) (Table S2).

For those SNPs significantly explaining RWI and ELI variability (see Factorial regression subsection), we re-examined their associated gene functions (previously described in Santini et al., 2021) through Blast tools as sequence databases are continuously updated. The original SNP sequence (Pinosio et al., 2014) was implemented in the Blastn tool and the best match from Blastn was then used as the input in the Blastx tool. The best result given by Blastx was then used to identify the protein function of every SNP based on the protein knowledgebase UniProtKB (Boutet et al., 2016).

Statistical analyses

General analyses of variance.

To test for differences in tree size among populations, the following mixed model was fitted to trunk diameter records:

$${\displaystyle \begin{array}{l}{\displaystyle y_{iklm}=\mu +P_i+B_k+{(PB)}_{ik}+T{\left(P\right)}_{il}+C_m+e_{iklm}}\end{array}}$$

(1)

where y_iklm is the observation of the lth tree of the ith population in the kth replicate and mth column, μ is the general mean, P_i is the fixed effect of the ith population, B_k is the fixed effect of replicate k, (PB)_ik is the random interaction between the ith population and kth replicate, T(P)_il is the fixed effect of the lth tree nested to the ith population, C_m is the fixed effect of column m and e_iklm is the random residual effect of the interaction between the lth tree nested to the ith population, the mth column and the kth replicate. Equation (1) was fitted to the 130 sampled trees plus those remaining trees belonging to the 23 studied populations after the thinning treatment (i.e. the total of four trees per experimental unit). In all cases, tree diameter at breast height (DBH; 1.30 m height) was measured at the time of thinning.

A second (fixed) model (herein general model) was fitted to each tree growth trait (RWI, ELI) testing for population differentiation in annual radial growth patterns as follows:

$${\displaystyle \begin{array}{l}{\displaystyle y_{ijl}=\mu +P_i+T{\left(P\right)}_{il}+Y_j+{\left(PY\right)}_{ij}+e_{ijl}}\end{array}}$$

(2)

where y_ijl is the observation of the lth tree of the ith population in the jth year, μ is the general mean, P_i is the fixed effect of the ith population, T(P)_il is the fixed effect of the lth tree nested to the ith population, Y_j is the fixed effect of the jth year, (PY)_ij is the fixed effect of interaction between the ith population and the jth year, and e_ijl is the random residual effect of the interaction between the lth tree nested to the ith population and the jth year. Because all tree-ring data (RWI and ELI) were detrended (i.e. they shared the same mean value among trees), the terms P_i and T(P)_il do not capture any significant variation (as neither do the design factors block and column, which can be saved). This implies that the variability in tree-ring records can be attributed exclusively to interannual fluctuations, Y_j, the interaction between the effects of population and year, (PY)_ij, and the residual term, e_ijl.

Factorial regression.

The general Equation (2) was expanded to include explicit genetic and environmental (i.e. climatic) covariables to the levels of the tree and year effects, respectively. To this end, we fitted two alternative G × E models using factorial regression approaches (van Eeuwijk et al., 2005,Voltas et al., 2005), which principally focused on the partition of the population by year effects. The first model investigated the genetic basis [SNP (or gene) g₁, …., g_m] of interannual variability in RWI and ELI (herein genetic model; Equation 3). The second model aimed at identifying differential population sensitivities of radial growth and early- to latewood ratios to climate drivers (z₁, … , z_p) (herein climate model; Equation 4). The order in which covariables are included in the model is relevant for the amount of G × E accounted for by each covariable (provided the explanatory variables are not completely orthogonal). Therefore, we progressively included the different covariables in the models sorted by the amount of G × E explained by each of them until the last significant covariable entered the model, ideally leaving a non-significant population by year residual. Type I (sequential) sum of squares was used for hypothesis testing (Nelder, 1994). In this way, population by year effects were adjusted for genetic effects (Equation 3) or climate factors (Equation 4) using covariables as follows (van Eeuwijk et al., 2005):

$${\displaystyle \begin{array}{l}{\displaystyle y_{ijl}=\mu +P_i+T{\left(P\right)}_{il}+Y_j+{\displaystyle \underset{n=1}{{\displaystyle \sum ^N}}}{g}_{ln}{Y}_{jn}+{\left(PY\right)}_{ij}+e_{ijl}}\end{array}}$$

(3)

$${\displaystyle \begin{array}{l}{\displaystyle y_{ijl}=\mu +P_i+T{\left(P\right)}_{il}+Y_j+{\displaystyle \underset{p=1}{{\displaystyle \sum ^P}}}{\beta }_{ip}{z}_{jp}+{\left(PY\right)}_{ij}+e_{ijl}}\end{array}}$$

(4)

where y_ijl is the observation of the lth tree of the ith population in the jth year, μ is the general mean, P_i is the fixed effect of the ith population, T(P)_il is the fixed effect of the lth tree nested to the ith population, Y_j is the fixed effect of the jth year, g_lnY_jn refers to the nth gene or covariable (SNP loci; qualitative factor with three levels or genotypes: aa, ab and bb) of the lth tree of the ith population interacting with the jth year in Equation (3), β_ip is the ith population sensitivity to the pth climatic covariable z for year j in Equation (4), (PY)_ij is the (residual) fixed effect of interaction between the ith population and the jth year, and e_ijl is the random residual effect of the interaction between the lth tree nested to the ith population and the jth year. To facilitate interpretation of the β_i terms, the climate covariables were centred to zero means. As an extension of Equation (3), the population structure (Pstr, K = 2; Santini et al., 2021) was also incorporated as a genetic covariable as an additional Pstr × Y term prior to SNP testing.

The genetic model (Equation 3) was further expanded to determine the type of gene action of relevant SNPs underlying G × E, that is, additive (A) or dominance (D) effects as well as possible first-order epistasis (van Eeuwijk et al., 2002; Calleja-Rodríguez et al., 2021). Additive and dominance effects were sequentially tested for each SNP individually in Equation (3) by creating dummy variables as follows: each SNP loci was coded as (–1, 0, +1) (accounting for additive effects) and (0, +1, 0) or (0, –1, 0) (accounting for dominance effects) at the tree level (van Eeuwijk et al., 2002). Once the type of gene action was identified as significant (either additive or dominant), relevant SNP loci were coded as (0, 1, 2) (additive) or otherwise (1, 1, 0) or (0, 1, 1) (dominant) and added to the final Equation (5) below.

First-order epistatic interactions were also evaluated independently for each significant SNP by considering potential gene interactions with other SNPs that were not necessarily relevant themselves for G × E explanation. For this, the SNP term was initially considered a qualitative factor. Additive and dominance genetic effects were then coded as above per SNP locus and successively tested by considering any possible effect combination between the two SNPs involved in the epistatic effect (additive × additive, additive × dominant, dominant × additive or dominant × dominant), as follows:

$$\begin{array}{cc}\underset{n=1}{\sum ^N}{g}_{ln}{Y}_{jn}& =\underset{n=1}{\sum ^N}{g}_{\left(A\right)ln}{Y}_{jn}+\underset{n=1}{\sum ^N}{g}_{\left(D\right)ln}{Y}_{jn}\\ & +\underset{n=1}{\sum ^N}{g}_{\left(A\right)ln}{g}_{\left(A\right)l{n}^{\prime }}{Y}_{jn}\\ & +\underset{n=1}{\sum ^N}{g}_{\left(D\right)ln}{g}_{\left(A\right)l{n}^{\prime }}{Y}_{jn}\\ & +\underset{n=1}{\sum ^N}{g}_{\left(A\right)ln}{g}_{\left(D\right)l{n}^{\prime }}{Y}_{jn}\\ & +\underset{n=1}{\sum ^N}{g}_{\left(D\right)ln}{g}_{\left(D\right)l{n}^{\prime }}{Y}_{jn}\end{array}$$

(5)

where g_ln refers to the nth gene (SNP, quantitative factor, coded as additive, g_(A), or dominant, g_(D)) of the lth tree interacting with the nth gene of the same tree (also coded as additive, g_(A), or dominant, g_(D)) and the jth year. Non-additive and epistatic genetic effects can largely influence phenotypes (Holliday et al., 2012; Du et al., 2015; Calleja-Rodríguez et al., 2021) and it is, therefore, important to consider those effects in G × E models. Since all the genetic factors present in the models follow a hierarchical order (sequential testing), the significance of genetic effects (both single SNP and epistatic) is contingent on the significance of the additive type of gene action. For example, if the SNP effect is dominant, both the additive and the dominant effect would result in significance in the model.

Additionally, the different SNPs and climate factors previously identified in the genetic and climate models [Equations (3) and (4) respectively] were tested as pair combinations of one SNP and one climate factor. If significant, the common variability explained by these factors was included – as an extension of Equations (3) and (4) – into models that incorporated a cross-product involving one genetic and one climate variable underlying the population by year and the tree nested to the population by year interaction terms. In this case, population sensitivities β_i to an environmental variable z_j were replaced by a constant, c, times a genetic (SNP) variable, x_i, (β_i = c x_i), where the constant c was estimated from the data. These cross-products allowed us to estimate how particular allelic substitutions affected either the interannual tree growth or the proportion of early- to latewood for every unit change in relevant climate factors.

All models were fitted using the MIXED procedure of SAS/STAT (Littell et al., 1998).

Climate–growth associations at the species level and their dependencies on climate at the origin of populations.

To determine the main climate drivers of interannual growth for Aleppo pine across populations (i.e. at the species level) we calculated bootstrapped Pearson correlations between the RWI or ELI master chronology and monthly climate data during 2000–2019. These relationships were analysed from September of the previous year to October of the current year using the treeclim R package (Zang and Biondi, 2015). Moreover, population sensitivities to relevant climate factors [β_i terms as defined in Equation (4)] were correlated with climatic characteristics at the origin of populations.

RESULTS

Growth differences among populations

We found significant differences (P < 0.05) in DBH among populations at the trial site (Supplementary Data Table S3). The populations showing the highest DBH values originated from Greece (Kassandra), southern Italy (Litorale Tarantino) and southern Spain (Monovar), with mean (± s.e.) trunk diameters of 18.4 ± 1.4, 16.9 ± 1.4 and 16.7 ± 1.3 cm, respectively. On the other hand, the populations with the lowest DBH values were from Tunisia (Tabarka) and Mallorca Island (Santanyí and Palma de Mallorca), with mean values of 11.3 ± 1.4, 11.8 ± 1.5 and 12.1 ± 1.2 cm respectively (Table S4). The relative difference between the most and least grown populations was 63 %.

Climate–growth relationships

At the species level (i.e. across populations), RWI was negatively correlated with previous September and current May and June temperatures, as well as with previous November precipitation and SPEI1. On the other hand, January temperature, May precipitation, May and June precipitation as well as the SPEI1, and SPEI6 of late spring and summer (May to August) had positive effects on RWI (Fig. 2). Similar climate factors affected the proportion of latewood relative to earlywood (ELI). In particular, high temperatures of previous September and current May and June, along with high previous November precipitation and SPEI1, decreased ELI; conversely, high current May and June precipitation and SPEI1 and high current summer SPEI6 (June to August) increased ELI (Fig. 2).

Fig. 2.

Bootstrapped correlation coefficients between the master chronology (i.e. obtained across populations) of indexed mean tree ring-width series (green bars) or early- to latewood ratio indices (orange bars) and 6-monthly climate factors calculated from previous September to current October (months of the previous year are indicated with lowercase letters) during the period 2000–2019. Climate factors were mean temperature (T_mean), maximum temperature (T_max), minimum temperature (T_min), precipitation (P), and SPEI at 1- and 6-month scales (SPEI1 and SPEI6, respectively). Filled bars indicate significant bootstrapped correlations (P < 0.05).

Genetic and climate factors involved in G × E

There was significant year-to-year variability in RWI and ELI that was contingent on population (i.e. a significant population by year interaction, Fig. 3). This interaction explained 10 and 16 % of the total sum of squares (SS) for RWI (Table 1) and ELI (Table 2), respectively. Four SNPs together explained 14 % of the interannual differentiation in RWI among trees (or G × E, subsumed in population by year and tree nested to the population by year interaction SS; Table 1). These were SNP₂₀₁, SNP₁₅₁, SNP₁₃₃ and SNP₉, along with two epistatic effects encompassing SNP₁₅₁ and SNP₁₃₃, as well as SNP₁₅₁ and SNP₉. The inclusion of population structure in the genetic models did not substantially change the results (Supplementary Data Tables S5 and S6). Alternatively, four climate factors explained 38 % of the population by year interaction SS for RWI (Table 1): May maximum temperature (accounting for 12 % of the interaction SS), previous year November precipitation (10 %), October SPEI6 (9 %) and previous September SPEI1 (8 %).

Table 1.

Analysis of variance (General model) and factorial regression modelling of genotype by environment (G × E) interaction effects using either molecular markers (Genetic model) or climate information (Climate model) of indexed ring width (RWI) involving 130 individuals belonging to 23 populations of Aleppo pine grown in a common garden in Altura (Spain).

	RWI
Source	d.f.	SS	MS	F-value	Pr (> F)	R2 (%)
	General model
Population (P)	22	0.14	0.007	0.12	1.000	0.1
Tree (within P)	107	0.63	0.006	0.11	1.000	0.2
Year (Y)	19	131.71	6.932	124.92	< 0.001	48.5
P × Y	418	26.34	0.063	1.14	0.043	9.7
Error (Ɛ)	2033	112.82	0.055			41.5
	G × E partition (genetic model)
SNP201(A) × Y	19	2.60	0.137	2.58	< 0.001	1.9
SNP151(A) × Y	19	2.29	0.121	2.28	0.001	1.6
SNP151(D) × Y	19	1.34	0.071	1.34	0.150	1.0
SNP133(A) × Y	19	1.00	0.052	0.99	0.469	0.7
SNP133(D) × Y	19	1.22	0.064	1.21	0.236	0.9
SNP9(A) × Y	19	0.86	0.045	0.85	0.646	0.6
SNP151(A) × SNP133(A) × Y	19	1.60	0.084	1.59	0.050	1.1
SNP151(D) × SNP133(A) × Y	19	1.57	0.083	1.56	0.058	1.1
SNP151(A) × SNP133(D) × Y	19	1.58	0.083	1.57	0.056	1.1
SNP151(D) × SNP133(D) × Y	19	1.68	0.088	1.67	0.035	1.2
SNP151(A) × SNP9(A) × Y	19	1.63	0.086	1.62	0.044	1.2
SNP151(D) × SNP9(A) × Y	19	1.92	0.101	1.91	0.010	1.4
P × Y	418	24.49	0.059	1.11	0.090
Error (Ɛ)	1786	94.63	0.053
	G × E partition (climate model)
P × Tmax_May	22	3.24	0.147	2.68	< 0.001	12.3
SNP201(A) × Tmax Maybetween P	1	0.14	0.141	2.56	0.110	0.5
SNP201(A) × Tmax May within P	1	0.68	0.676	12.28	< 0.001	0.6
P × P Nov(−1)	22	2.52	0.115	2.08	0.002	9.6
SNP201(A) × P Nov(−1) between P	1	0.22	0.225	4.08	0.044	0.9
SNP201(A) × P Nov(−1) within P	1	0.00	0.000	–	–	0.0
P × SPEI6 Oct	22	2.29	0.104	1.89	0.008	8.7
P × SPEI1 Sep(−1)	22	2.03	0.092	1.67	0.026	7.7
P × Y	330	16.26	0.049	0.90	0.899
Error (ɛ)	2031	111.78	0.055

Significant terms are shown in bold (P < 0.05).

Abbreviations: T_max May: maximum temperature in May; P Nov₍₋₁₎: previous November precipitation; SPEI6 Oct: SPEI at 6-month scale of October; SPEI1 Sep₍₋₁₎: previous September SPEI at 1-month scale.

Table 2.

Analysis of variance (General model) and factorial regression modelling of genotype by environment (G × E) interaction effects using either molecular markers (Genetic model) or climate information (Climate model) of indexed early- to latewood ratio (ELI) involving 130 individuals belonging to 23 populations of Aleppo pine grown in a common garden in Altura (Spain).

	ELI
Source	d.f.	SS	MS	F-value	Pr (> F)	R2 (%)
	General model
Population (P)	22	0.05	0.002	0.02	1.000	0.0
Tree (within P)	107	0.18	0.001	0.01	1.000	0.0
Year (Y)	19	86.56	4.556	38.83	< 0.001	22.2
P × Y	418	63.97	0.153	1.30	< 0.001	16.4
Error (Ɛ)	2033	238.51	0.120			61.3
	G × E partition (genetic model)
SNP159(A) × Y	19	4.46	0.244	2.09	0.004	1.5
SNP159(D) × Y	19	3.01	0.160	1.40	0.117	1.0
SNP133(A) × Y	19	2.95	0.160	1.37	0.131	1.0
SNP133(D) × Y	19	2.45	0.130	1.14	0.304	0.8
SNP159(A) × SNP133(A) × Y	19	5.69	0.300	2.64	0.001	1.9
SNP159(D) × SNP133(A) × Y	19	1.92	0.100	0.89	0.593	0.6
SNP159(A) × SNP133(D) × Y	19	4.59	0.240	2.13	0.003	1.5
SNP159(D) × SNP133(D) × Y	19	1.14	0.060	0.53	0.951	0.4
P × Y	399	60.87	0.150	1.35	< 0.001
Error (Ɛ)	1900	215.22	0.110
	G × E partition (climate model)
P × SPEI1 Nov(−1)	22	6.02	0.274	2.34	< 0.001	9.4
P × P Dec(−1)	22	6.00	0.273	2.33	< 0.001	9.4
P × Tmean Jun	22	4.90	0.223	1.91	0.007	7.7
SNP159(A) × Tmean Jun between P	1	0.00	0.000	–	–
SNP159(A) × Tmean Jun within P	1	0.93	0.931	7.96	0.005	0.3
P × Tmin Mar	22	5.32	0.242	2.07	0.003	8.3
P × SPEI1 Sep	22	4.25	0.193	1.65	0.029	6.6
P × Y	309	37.46	0.121	1.04	0.331
Error (ɛ)	2031	237.58	0.117

Significant terms are shown in bold (P < 0.05).

Abbreviations: SPEI1 Nov(-1): SPEI at 1-month scale of previous November; P Dec(-1): previous November precipitation; T_mean Jun: mean temperature of Jun; T_min Mar: minimum temperature of March; SPEI1 Sep: SPEI at 1-month of September.

Fig. 3.

Interannual population variability of raw tree-ring width (RW; mm) (A), tree ring width index (RWI) (B) and indexed early- to latewood ratio (ELI) (C) of 23 populations of Aleppo pine grown in a common garden in Altura (Spain). Each population is represented by a line of different colour; population acronyms are listed in the key and described in Supplementary Data Table S1.

For ELI, SNP₁₅₉ and SNP₁₃₃ and their epistatic effects partly explained interannual tree variability, accounting for 9 % of G × E (Table 2). In turn, populations reacted differently to five climate factors, which explained 41 % of the population by year interaction term (Table 2): previous November SPEI1 (accounting for 9 % of the population by year interaction), previous December precipitation (9 %), June mean temperature (8 %), March minimum temperature (8 %) and September SPEI1 (7 %) (Table 2). Both the genetic and climate models had non-significant population by year interaction residuals (except for the genetic model for ELI), thereby indicating the suitability of the selected molecular markers and climate factors for the explanation of differential population performance (Tables 1 and 2).

Population sensitivities to climate factors and relationships with climate at the origin

Five populations showed a significant RWI sensitivity to May maximum temperature, either positive (two populations) or negative (three populations) relative to a hypothetical average population. Another five populations exhibited a significant (positive or negative) RWI sensitivity to previous November precipitation, and three populations showed a significant RWI sensitivity to October SPEI6 (Supplementary Data Table S7). Regarding ELI, six populations had a significant sensitivity (positive or negative) to previous November SPEI1, while five populations showed a significant sensitivity to previous December precipitation, four populations to previous September SPEI1, three populations to June mean temperature and five populations to March minimum temperature (Table S8).

We found significant associations between population sensitivities to climate at the trial site and climate at the origin of these populations for RWI (Table 5). Temperature was the climate factor most closely related to these sensitivities, with populations from warmer environments (i.e. with higher MAT and MST) and experiencing more severe winter conditions (i.e. lower MWT) being less sensitive to high May temperatures, high previous November precipitation and dry growing seasons (October SPEI6) (Table 5, Fig. 4). On the other hand, populations subject to higher precipitation seasonality (Fig. 4) and lower summer to annual precipitation ratio at origin, as well as those belonging to the southernmost areas of the species distribution range, were more sensitive to dry growing seasons than their wetter counterparts (Table 5). For ELI we did not find relevant associations between population sensitivities to climate and climate at a population’s origin (Supplementary Data Table S9).

Table 5.

Correlation coefficients between population sensitivities (β values) of indexed ring width (RWI) to four climate factors at the trial site [in row order: maximum temperature of May (T_max May), precipitation of previous November (P Nov₍₋₁₎), SPEI at a 6-months scale of October (SPEI6 Oct) and SPEI at a 1-month scale of previous November (SPEI1 Sep₍₋₁₎)] and selected climate factors at population origin.

	Latitude	MAT	MST	MWT	TmaxW	TminC	MAP	PS	PDQ	PsP
Tmax May	−0.22 (0.327)	0.44 (0.039)	0.40 (0.057)	0.41 (0.048)	0.22 (0.302)	0.39 (0.066)	−0.23 (0.300)	0.08 (0.725)	−0.15 (0.497)	−0.07 (0.767)
P Nov(−1)	−0.34 (0.114)	0.39 (0.066)	0.42 (0.046)	0.31 (0.149)	0.32 (0.132)	0.24 (0.271)	−0.15 (0.484)	0.40 (0.059)	−0.32 (0.140)	−0.38 (0.077)
SPEI6 Oct	0.41 (0.053)	−0.62 (0.002)	−0.61 (0.002)	−0.55 (0.007)	−0.37 (0.080)	−0.54 (0.008)	0.48 (0.829)	−0.63 (0.001)	0.48 (0.020)	0.44 (0.037)
SPEI1 Sep(−1)	0.07 (0.758)	−0.21 (0.324)	−0.10 (0.658)	−0.28 (0.196)	0.05 (0.835)	−0.29 (0.179)	0.13 (0.179)	−0.02 (0.940)	<0.01 (0.982)	0.10 (0.643)

Correlation coefficients (r) and their associated probabilities (P, in parentheses) are shown; bold characters indicate probabilities < 0.05.

Abbreviations: MAP: mean annual precipitation; MAT: mean annual temperature; MST: mean summer temperature; MWT: mean winter temperature; T_maxW: maximum temperature of the hottest month; T_minC: minimum temperature of the coldest month; PS: precipitation seasonality; PDQ: mean precipitation of the driest quarter; PsP: summer to annual precipitation ratio.

Fig. 4.

Correlations between population sensitivities of ring width (RWI) to October SPEI-6 at the trial site and selected climate factors (MAT, mean annual temperature; PS, precipitation seasonality) at the origin of 23 Aleppo pine populations grown in a common garden in Altura (Spain). Population acronyms are defined in Supplementary Data Table S1.

Genes interacting with climate for the explanation of G × E

All individual SNPs included in the genetic models showed additive effects for both RWI and ELI, while two epistatic interactions for RWI had dominant × dominant or dominant × additive gene actions and one epistatic interaction for ELI had a dominant × additive gene action. For RWI, we found significant cross-products between SNP₂₀₁ and two climate factors (May maximum temperature and previous November precipitation). However, only the cross-product between SNP₂₀₁ and previous November precipitation significantly explained the interaction between populations and years (albeit only accounting for 1 % of the interaction SS, Table 1). For ELI, SNP₁₅₉ significantly interacted with June mean temperature; however, it explained only interannual variability of trees nested to populations, but not among populations (Table 2).

For the above-mentioned SNPs, the effect of allelic substitutions as related to climate was estimated. In particular, replacing the C allele with the A allele in SNP₂₀₁ decreased RWI by 0.029 and 0.001 (standardized) units per every 1 °C increase in May maximum temperature and every 1 mm increase in previous November precipitation, respectively (Table 4). In turn, replacing the G allele with the A allele in SNP₁₅₉ increased ELI by 0.060 (standardized) units per every 1 °C increase in June mean temperature (Table 4). Known gene functions of every SNP with particular relevance in G × E models are annotated in Table 3. For most SNPs, it was possible to retrieve the related biological functions of their associated proteins.

Table 3.

Annotation of the homologous protein and its biological function of single nucleotide polymorphisms (SNPs) that were informative of either indexed ring width (RWI) or early- to latewood ratio (ELI) variation in the partition of genotype by year effects using factorial regression.

SNP code	Known annotation	Species	Biological function	Trait (GWAS) Santini et al. (2021)
SNP9	Hypothetical protein 2_1014_01	Pinus pinaster	–	Crown area
SNP133	Hypothetical protein 2_7803_01: glycosidase, hydrolase	Pinus taeda	Carbohydrate metabolic process	Water content; height
SNP151	Peroxisome membrane protein 11C	Morella rubra	Peroxisome fission; photomorphogenesis	Leaf area
SNP159	Putative calcium-dependent protein kinase transferase	Cupressus sempervirens	Mediates responses to abiotic and biotic stressor	Height
SNP201	No similarity found	Pinus halepensis	–	Photosynthetic pigments

Table 4.

Effect of allelic substitutions of selected SNPs on indexed ring width (RWI) and early- to latewood ratio (ELI) for every unit change in (standardized) relevant climate factors following the identification of significant cross-products in factorial regression G × E models of Aleppo pine grown in a common garden in Altura (Spain).

Trait	Estimate	SE	t-value	Pr > |t|
Ringwidth (RWI)
SNP201(A) (C→A) × P Nov(−1)	−0.001	0.000	−2.020	0.044
SNP201(A) (C→A) × Tmax May	−0.029	0.009	−3.370	0.001
Early- to latewood ratio (ELI)
SNP159(A) (G→A) × Tmean Jun	0.060	0.021	2.820	0.005

Significant terms are shown in bold (P < 0.05).

Abbreviations: P Nov₍₋₁₎: previous November precipitation; T_max May: maximum temperature of May; T_mean Jun: mean temperature of June.

DISCUSSION

Aleppo pine harbours a large variability in life-history traits related to growth, defence and reproduction (Santos-del-Blanco et al., 2013), which results in a complex adaptive syndrome that tailors individual performances to the array of conditions encountered by the species (Sbay and Zas, 2018; Santini et al., 2019). It is therefore expected that some populations may respond better to climate instability and warming than others because of the combined effects of local adaptation and differential phenotypic plasticity on functional traits (Voltas et al., 2018; Hevia et al., 2020; Patsiou et al., 2020). Factorial regression provided first-hand information on the genetic basis of differentiation in annual growth responses and their climate drivers for this emblematic Mediterranean conifer. The major novelty of this approach is that it quantified with precision growth variability among individuals explained by particular SNPs while characterizing the relationship between the effect of allele substitutions in dendrophenotypes and climate. In this regard, this approach can be regarded as complementary to methods that analyse high-dimensional genomic information in relation to environmental variables but where phenotypic information is absent (e.g. redundancy analysis; Capblancq and Forester, 2021; Varas-Myrik et al., 2022), or that investigate long-term or specific-year tree-ring traits searching for genotype–phenotype associations but disregarding G × E (Housset et al. 2018).

Climate responses at the species level

Populations with larger trunk diameters were from coastal eastern Mediterranean areas and cool climates, while populations that showed smaller sizes originated from warm and dry areas (Voltas et al., 2018). This observation anticipates that low water availability and high temperatures are the main climatic constraints for Aleppo pine growth (de Luis et al., 2013). At the same time, both factors probably had a variable effect on growth depending on the geographical origin of populations (Patsiou et al., 2020), as suggested by significant population by year interactions for ring width and early- to latewood ratio.

At the whole species level (i.e. across populations), annual growth was mainly controlled by late spring (May–June) climate (Pasho et al., 2012; Novak et al., 2013). Water stress in spring, modulated by temperature, is known to negatively affect xylogenesis, thereby inducing a decline in secondary growth as a probable result of meristematic constraints (de Luis et al., 2013; Puri et al., 2015; Gazol et al., 2017; Hevia et al., 2020). In addition, RWI and ELI were influenced negatively by previous September and positively by January temperatures. The former effect is often interpreted as an indication of reserves depletion, hence decreasing secondary growth (Kagawa et al., 2006; Choury et al., 2017) and early- to latewood proportion. The latter suggests that this thermophilic pine is constrained by cold winters, which influences its ability to resume cambial activity earlier in the growing season (Camarero et al., 2010; Housset et al., 2018). The negative response of RWI and ELI to previous autumn precipitation (November), earlier reported for the region (Shestakova et al., 2017; Hevia et al., 2020), might be related to the impairment of photosynthesis and, hence, of carbon storage due to prolonged cloudiness. This would lead to a decrease in earlywood and total ring width (Camarero et al., 2010; Royo-Navascués et al., 2021).

How do populations differentially respond to annual climate variability?

Our findings revealed the existence of different radial growth patterns among Aleppo pine populations. This array of growth responses could be related to particular climate factors at the trial site. Therefore, in addition to genetic variation in tree diameter among populations, adaptation in Aleppo pine appears to be partly driven by differential plasticity in secondary growth (Voltas et al., 2018; Hevia et al., 2020; Patsiou et al., 2020). This differential plasticity, as explained by G × E interaction effects, is better related to water availability than to temperature (compare their relevance in the explanation of G × E; cf. Tables 1 and 2). This observation broadens previous results on natural stands (as opposed to common garden observations), which anticipate such plastic effects in wood traits for this species (de Luis et al., 2013; Novak et al., 2013).

Some populations were particularly sensitive to one or several climate factors at the trial site. A comparison with their climate responses at origin (i.e. in natural stands), as reported in previous studies, provides hints on the adaptive significance of such sensitivities. For example, Tabarka (Tunisia) is a population from a sub-humid environment of the Maghreb that was extremely dependent on water availability at the trial site. This suggests that carbohydrates synthetized under favourable conditions (autumn) may be particularly important for spring growth resumption, and is consistent with in situ observations indicating that drought during the previous autumn through spring is most limiting for Tunisian Aleppo pine forests (Bachtobji Bouachir et al., 2017). Another population from a sub-humid coastal area (Kassandra, Greece) showed high sensitivity to elevated temperatures during late spring. This is in agreement with a previous study indicating that very warm springs strongly limit Aleppo pine growth in Greece due to excessive evapotranspiration rates (Papadopoulos et al., 2001). This performance can be explained by the extreme plasticity of Greek populations adjusting the timing of earlywood formation (Hevia et al., 2020). Conversely, the sensitivity to hot summers in terms of decreased early- to latewood ratio detected for dry Iberian populations (e.g. Benicàssim) points to a conservative water-use strategy: a strongly reduced earlywood in unfavourable years is probably indicative of an earlier growth cessation before peak summer and hence higher drought resistance (Hevia et al., 2020; Royo-Navascués et al., 2021).

Variation in population sensitivities to the above-mentioned climate factors showed significant associations with climate at origin, suggesting distinct adaptations (George et al., 2019; Depardieu et al., 2020). The observation that populations from warm areas present less growth reductions under high May temperatures probably reflects specific adaptation to hot springs (Housset et al., 2018). On the other hand, populations from similarly warm areas but experiencing high precipitation seasonality and low summer rainfall were the most growth-sensitive to dry vegetative periods. This counterintuitive finding may be attributed to the combined effect of traits involved in local adaptation to drought, i.e. a particular realization of the species’ adaptive syndrome as materialized in a particular life-history strategy (Santini et al., 2019). Hence, these populations may present a strong summer cambial dormancy, thereby decreasing tracheid formation and thus secondary growth during dry periods, improving resistance to embolism (Camarero et al., 2010; De Luis et al., 2011). This agrees with intra-annual density fluctuation (IADF) records from a common garden experiment, in which populations from drier areas experienced more IADFs (Hevia et al., 2020).

Interpreting the genetic basis of secondary growth and its climate dependencies

The differential growth responses among individuals were related to allele-specific expressions of a subset of candidate genes, with four and two SNPs partially explaining annual variability in ring width and the early- to latewood ratio, respectively. Despite the limited explanation of G × E effects (<15 %), our candidate gene approach proved to be partly effective for revealing the genetic basis of secondary growth. However, most growth variability related to differential gene expression did not bear a clear population structure, which prevents interpreting such variability in terms of ecotypic differentiation in phenotypic plasticity. In this regard, the traits analysed are very likely to be under strong polygenic control, as shown for tree height in the taxonomically closely related Pinus pinaster (de Miguel et al. 2022). Thus, it is likely that many more SNPs may have been involved in the variability of secondary growth, some with very low effect and thus statistically difficult to detect (de Miguel et al. 2022). Incorporating a substantially larger amount of candidate genes into G × E modelling approaches that could effectively integrate this high-dimensional genetic complexity [e.g. partial least squares estimation (Vargas et al., 2006) or random reaction norm models (Jarquín et al., 2014)] could capitalize upon the increasing wealth of genomic and environmental information available, particularly for pine species with large and complex genomes. On the other hand, factorial regression may be particularly adequate for precise testing and quantification of the role of specific genes of annotated function in the determination of dendrophenotypes (or other longitudinal data) and, hence, for clarifying the existing relationships between genomes and phenomes in forest trees.

The interplay between gene markers and climate indicated that two SNPs had putative roles in tree adaptation to climate. SNP₂₀₁ (related to ring width) interacted with previous November precipitation and May maximum temperature. SNP₂₀₁ has been referred to as the P. halepensis transcriptome (Pinosio et al., 2014), but the molecular and biological function of its associated gene has not been described so far. Notably, a previous GWAS conducted in the same common garden (Santini et al., 2021) linked SNP₂₀₁ to variation in photosynthetic rate among trees. Alternatively, SNP₁₅₉ (related to early- to latewood ratio) interacted with June mean temperature. Calcium-dependent protein kinases (CDPKs) are associated with a candidate gene corresponding to SNP₁₅₉. CDPKs constitute a large multigene family involved in metabolic, ion flux and gene expression alteration and related to phytohormone activity, such as the modulation of GA3 homeostasis (Schulz et al., 2013). CPDKs have been shown to play a role in many physiological processes, including development and growth (Asano et al., 2012; Boudsocq et al., 2013). CDPKs are also recognized as positive regulators to environmental stresses (Schulz et al., 2013), with the G allele of SNP₁₅₉ seemingly providing adaptation to high temperatures during the peak growing season for Aleppo pine.

The remaining SNP loci relevant for G × E explanation could not be related to any explicit climate factor, which may indicate the need to screen for alternative temporal windows or environmental records to fully decipher the genetic basis of differential growth plasticity in the species as related to climate. Regardless, the biological function of these SNPs is known (except SNP₉), and it is related to growth traits in Aleppo pine (Santini et al., 2021; Table 3). In particular, the gene associated with SNP₁₅₁ codes for proteins of the PEX family that are involved in photomorphogenesis, which can influence leaf development (Santini et al., 2020) and photosynthesis (Kaur et al., 2013). These mechanisms may affect the allocation of photosynthetic products to secondary growth throughout the year. On the other hand, the gene corresponding to SNP₁₃₃ (involved in G × E of both RWI and ELI) is associated with polygalacturonases influencing tissue development and defence signals (Rui et al., 2017; Gallego-Giraldo et al., 2020). Interestingly, our study showed that the influence of SNP₁₃₃ on ring width and early- to latewood ratio was in both cases dependent on the joint expression of genes associated with SNP₁₅₁ and SNP₁₅₉, with effects of double dominance and simple dominance, respectively.

To conclude, this study contributes to improving our knowledge on the genetic basis and climate controls of growth variation of a widespread Mediterranean conifer. The different interannual responses among Aleppo pine populations could be described by climate factors mainly related to water availability, and could be successfully related to differential SNP expression of a subset of candidate genes. We therefore show how tree-ring phenotypes obtained from common garden experiments in combination with a candidate-gene approach constitute a great opportunity to disentangle the importance of genetic and environmental effects determining tree adaptation. While genomic data are now widely used as a source of insight into adaptation patterns for non-model species (Sork et al., 2013; Fitzpatrick and Keller, 2015; Wadgymar et al., 2017; Mahony et al., 2020), this is, to the best of our knowledge, the first attempt to integrate genetic and environmental information into statistical G × E models for a forest tree species by considering longitudinal records of ring width and early- to latewood ratio as phenotypic traits. In particular, a major novelty of this study is the possibility to quantify with precision the phenotypic changes attributable to SNPs associated with candidate genes and the associated effect of allelic substitutions in relation to climate variables. Future studies could consider alternative environmental variables that may also drive secondary growth of Aleppo pine (e.g. soil moisture, wind) while integrating a larger set of genetic markers associated with candidate genes. The precise molecular mechanisms underlying the adaptation patterns showed by some individuals carrying certain allele substitutions could also be a target for future research.

SUPPLEMENTARY DATA

Supplementary data are available online at https://academic.oup.com/aob and consist of the following. Table S1. Characteristics of the 23 populations of Aleppo pine used in this study. Table S2. Annotation of the homologous protein and its known biological function of each SNP used in G × E models. Table S3. Analysis of variance of diameter at breast height of the 23 Aleppo pine populations used in this study. Table S4. Least squares means of the 23 Aleppo pine populations obtained from analysis of variance. Table S5. Factorial regression models of genotype by environment interaction effects of indexed ring width, with population structure as the genetic covariable, of 130 individuals belonging to 23 populations of Aleppo pine. Table S6. Factorial regression models of genotype by environment interaction effects of indexed early- to latewood ratio, with population structure as the genetic covariable, of 130 individuals belonging to 23 populations of Aleppo pine. Table S7. Estimates of population sensitivities to relevant climate factors according to factorial regression models. Table S8. Estimates of population sensitivities to relevant climate factors according to factorial regression models. Table S9. Correlation coefficients between population sensitivities to climate at the trial site and different climate factors at origin.

ACKNOWLEDGMENTS

We thank Ricardo Alía and Delphine Grivet (CIFOR-INIA) for providing the genomic data. The data supporting the results of this study are available at the CORA open repository (https://dataverse.csuc.cat/).

CONFLICT OF INTEREST

The authors declare that they have no conflicts of interest.

FUNDING INFORMATION

This work was partly supported by the Spanish Government- Ministerio de Ciencia, Innovación y Universidades; grant numbers AGL2015-68274-C3-3-R (MINECO/FEDER) and RTI2018-094691-B-C31 (MCIU/AEI/FEDER, EU). E. Lombardi was supported by a AGAUR FI-2020 pre-doctoral fellowship (with support from the Secretariat for Universities and Research of the Ministry of Business and Knowledge of the Government of Catalonia and the European Social Fund).