Genomic data and common garden experiments reveal climate‐driven selection on ecophysiological traits in two Mediterranean oaks

Improving our knowledge of how past climate‐driven selection has acted on present‐day trait population divergence is essential to understand local adaptation processes and improve our predictions of evolutionary trajectories in the face of altered selection pressures resulting from climate change. In this study, we investigated signals of selection on traits related to drought tolerance and growth rates in two Mediterranean oak species (Quercus faginea and Q. lusitanica) with contrasting distribution ranges and climatic niches. We genotyped 182 individuals from 24 natural populations of the two species using restriction‐site‐associated DNA sequencing and conducted a thorough functional characterization in 1602 seedlings from 21 populations cultivated in common garden experiments under contrasting watering treatments. Our genomic data revealed that both Q. faginea and Q. lusitanica have very weak population genetic structure, probably as a result of high rates of pollen‐mediated gene flow among populations and large effective population sizes. In contrast, common garden experiments showed evidence of climate‐driven divergent selection among populations on traits related to leaf morphology, physiology and growth in both species. Overall, our study suggests that climate is an important selective factor for Mediterranean oaks and that ecophysiological traits have evolved in drought‐prone environments even in a context of very high rates of gene flow among populations.

patterns of genetic divergence among populations are also largely impacted by neutral demographic processes, such as range shifts, gene flow and genetic drift, thus confounding inference of selection Villemereuil et al., 2022). In forest tree species, palynological and ancient DNA studies have provided strong evidence of range shifts triggered by Holocene postglacial warming, particularly in northern temperate regions (Davis & Shaw, 2001;Saltré et al., 2013;Giesecke & Brewer, 2018;Wagner et al., 2018).
The rate and scale of species migrations have differed across continents, mainly due to factors such as climate and the geographical configuration of corridors and barriers to dispersal. Specifically, the dynamics of temperate species were governed largely by longdistance expansions in central Europe, Asia and most areas from North America Petit, Brewer, et al., 2002;de Lafontaine et al., 2010;Cao et al., 2015). In contrast, many Mediterranean species from Europe and North America expanded at lower scales from multiple refugia (Grivet et al., 2006;López de Heredia et al., 2007). Migration routes are generally reflected in the geographical distribution of maternally inherited organelle haplotypes, as a result of demographic expansions and founding events via seed dispersal Pham et al., 2017;Sork, Fitz-Gibbon, et al., 2016;. Subsequent among-population gene flow, especially via long-distance pollen dispersal, tends to erode nuclear genetic differentiation produced by ancient migration (Savolainen et al., 2007;Kremer & Le Corre, 2012), interacting with divergent selection and drift to produce presentday patterns of population genetic divergence.
Despite extensive long-distance gene flow among forest tree populations, there is strong evidence that suggests the emergence and maintenance of adaptive population divergence after postglacial colonization , and references therein). For example, widely distributed European oaks exhibit strong population differences in traits related to growth rates, vegetative phenology and drought tolerance (Alberto et al., 2013;Torres-Ruiz et al., 2019;Sáenz-Romero et al., 2017) that have presumably occurred as a result of divergent selection during the last 8000 years despite high gene flow (Alberto et al., 2013;. Climate has been considered the main selective agent driving population genetic differentiation in forest tree species (Alberto et al., 2013;Aitken & Bemmels, 2016;Baughman et al., 2019;Ramírez-Valiente et al., 2022).
In particular, temperature and precipitation regimes have been related to the evolution of traits important for fitness (i.e., functional traits) across species' ranges (Alberto et al., 2013. Conversely, there are also well-documented examples of functional trait population divergence in animal and plant species, including forest trees, of a magnitude that is consistent with neutral demographic processes alone. For example, López-Goldar et al. (2019) showed that demographic processes could explain by themselves population differentiation in constitutive secondary metabolites in the Mediterranean pine Pinus pinaster, a species showing a marked genetic structure. Similarly, Vázquez-González et al. (2019) in a study on the same species found that genetic differentiation in stem anatomical traits could be largely explained by genetic drift.
Other studies have reported coincident geographical patterns in neutral and adaptive population genetic variation for species with weaker spatial genetic structure as well. This is the case, among others, of the Mediterranean oak Quercus suber, for which a longitudinal pattern of population variation in fitness-related traits, phenology and growth architecture (e.g., Gandour et al., 2007;Ramírez-Valiente et al., 2014, Sampaio et al., 2016 mirrors the neutral genetic structure in both nuclear and organelle genomes (Magri et al., 2007;Pina-Martins et al., 2019).
Mediterranean regions are characterized by summers with low precipitation, which creates a period of water stress for plants (Larcher, 2000). Differences in the length and severity of the dry season have been suggested to exert differential selective pressures driving genetic divergence in traits related to drought tolerance and resource-use strategies (e.g., Gratani et al., 2003;Ramírez-Valiente, Lorenzo, et al., 2009;Ramírez-Valiente 2010;Andivia et al., 2012;Ghouil et al., 2020;Solé-Medina et al., 2022). However, for most Mediterranean species the extent to which nuclear genomic divergence and functional trait divergence may have been influenced by population demographic dynamics is still not well understood. In this study, we integrate genomic anal-  et al., 2009), and exhibits a range of leaf habits-from evergreen to deciduous-and growth forms-from shrub-like to large treesobservable across natural populations. Q. lusitanica is a shrub or small tree (seldom taller than 1 m), with evergreen or brevideciduous leaf habit and a distribution range restricted to the Atlantic coast of the Iberian Peninsula. It inhabits areas with oceanic influence and dry summers, except a northern isolated population that grows under a rainy temperate climate (Figure 1). By studying both species in the same common-garden experiment, potential interspecific differences in observed patterns of trait population divergence and in selection signals may be ascribed to species-specific genetic and selective factors, rather than to contemporary growing-environment effects. Both Q. faginea and Q. lusitanica are closely related phylogenetically to widespread temperate European white oaks such as Q. robur and Q. petraea . Analyses of chloroplast DNA (cpDNA) of white oaks have revealed four geographically structured haplotypic groups within the Iberian Peninsula shared across all species, including Q. faginea, which have been related to the existence of different glacial refugia followed by subsequent postglacial expansions and intra-and interspecific gene flow (Olalde et al., 2002;. Here, we test the hypothesis that climate has driven the divergent evolution of functional traits related to drought tolerance and growth among populations of both Q. faginea and Q. lusitanica, even in a probable scenario of high gene flow among conspecific populations inhabiting contrasting environmental conditions (Savolainen et al., 2007;. As observed in other oak species, we expect patterns of divergent selection to be particularly evident on morphological traits and growth rates Ramírez-Valiente et al., 2018). Specifically, we expect populations from drier areas to exhibit conservative strategies, with more sclerophyllous and smaller leaves, providing increased desiccation tolerance under low water potentials (Ramírez-Valiente et al., 2010. We also test the effect of experimental environment on inferred signals of selection. Despite genotype-by-environment interactions being ubiquitous for many traits and species (Matesanz & Ramírez-Valiente, 2019), only a few studies, mainly on animals, have tested the influence of experimental growing environments on selection inference (Palo et al., 2003;Hangartner et al., 2012;Goodrich et al., 2016). In our study, we expect that experimentally increasing water availability should help to reveal an acquisitive strategy for seedlings from mesic populations, enhancing the expression of faster growth relative to seedlings from xeric populations, which would suggest an adaptive increase in competitive ability under favourable conditions. Finally, we expect patterns of population divergence and trait selection to be stronger in Q. faginea than in Q. lusitanica, as a result of the more restricted climatic niche and lower climatic variability across the distributional range of the latter. Such potential differences between two closely related species with contrasting climatic niches and ecological strategies might improve our understanding of the role of environmental heterogeneity in driving local adaptation. We address the following specific questions: Do Q. faginea and Q. lusitanica exhibit population F I G U R E 1 Top panels show the location of the studied populations of Quercus faginea (a) and Q. lusitanica (b) in the Iberian Peninsula. Light and dark grey shaded areas in (a) indicate the Iberian distribution range of Q. faginea according to Caudullo et al. (2017) andDe Rigo et al. (2016), respectively. The shaded area in (b) shows the distribution range of Q. lusitanica estimated using species distribution modelling, where a higher probability of occurrence is represented by darker grey colours. Circles indicate populations with genomic and phenotypic data, squares indicate populations with only phenotypic data, and diamonds indicate populations with only genomic data. Middle panels represent the genetic structure of populations of (c) Q. faginea and (d) Q. lusitanica based on the Bayesian clustering method implemented in structure for the K values with the highest ln Pr(X|K). Admixture proportions are represented using pie charts, with each colour indicating a different genetic cluster. Pie chart size is proportional to the number of genotyped individuals at each location. Population codes are described in Table 1 We used stacks version 1.35 (Hohenlohe et al., 2010;Catchen et al., 2011Catchen et al., , 2013 to assemble our sequences into de novo loci and call genotypes, a pipeline that produces data sets consistent with those obtained using reference genome-based mapping approaches (Shafer et al., 2017;see also Fitz-Gibbon et al., 2017). Briefly, we demultiplexed and filtered reads for overall quality using the program process_radtags, retaining reads with a Phred score > 10 (using a sliding window of 15%), no adaptor contamination, and having an unambiguous barcode and restriction cut site. We screened raw reads for quality with fastqc version 0.11.5 (http://www.bioin forma tics.babra ham.ac.uk/proje cts/fastq c/) and trimmed all sequences to 129 bp using seqtk (Heng Li, https://github.com/lh3/seqtk) in order to remove low-quality reads near the 3′ ends. We assembled filtered reads de novo into putative loci with the program ustacks. We set the minimum stack depth (m) to three and allowed a maximum distance (M) of two nucleotide mismatches to group reads into a "stack." We used the "removal" (r) and "deleveraging" (d) algorithms to eliminate highly repetitive stacks and resolve over-merged loci, respectively.
We identified single nucleotide polymorphisms (SNPs) at each locus and called genotypes using a multinomial-based likelihood model that accounts for sequencing errors, with the upper bound of the error rate (ε) set to 0.2 (Hohenlohe et al., 2010;Catchen et al., 2011Catchen et al., , 2013. Then, we built a catalogue of loci using the cstacks program, with loci recognized as homologous across individuals if the number of nucleotide mismatches between consensus sequences (n) was ≤2.
Finally, we matched each individual data against this catalogue using the program sstacks and exported output files in different formats for subsequent analyses using the program populations. For all subsequent analyses, we exported only the first SNP per RAD locus and retained loci with a minimum stack depth ≥5 (m = 5), a minimum minor allele frequency (MAF) ≥0.01 (min_maf = 0.01) and represented in at least ̴ 80% of the populations (p = 12 for Q. faginea and p = 7 for Q. lusitanica) and 50% of the individuals within each population (r = .5).

| Common garden experiments
For the ecophysiological characterization, we used a common gar-  (Table 1), covering the distribution range of both species. Seeds were stored at 4°C. In spring 2018, we sowed the seeds in 2.3-L pots (CC-1028R; Sansan) filled with fine sand and placed them randomly within the greenhouse. Seedlings were maintained under ambient light and 10-25°C of temperature for 2 months after emergence. Then, seedlings were transplanted into double pots (i.e., two 2.3-L CC-1028R Sansan pots stacked on top of each other, after removing the bottom of the upper pot) with a total volume of 4.5 L and a length of 55 cm. We used a soil mixture of 1:1 (v/v) peat and fine sand. The experiment followed a complete block design with six blocks. Between one and five seedlings per family of each species were placed within each block, depending on seedling availability.
In total, the experiment included 1350 seedlings, 966 of Q. faginea and 384 of Q. lusitanica. After watering all seedlings every other day for 2 weeks, watering was withdrawn in three blocks to induce gradual water stress (dry treatment), whereas seedlings in the other three blocks were watered to field capacity every 2-3 days (wellwatered treatment). Volumetric water content (VWC) was measured once a week in seven randomly selected pots per block using time domain reflectometry (ProCheck, Decagon Devices, Inc.). In the well-watered treatment, VWC was maintained between 28.3% and 34.3% throughout the experiment, while it was gradually reduced down to an average low of 6.7% in the dry treatment at the end of the experiment. Additionally, for Q. faginea, we established a common garden experiment under field conditions (outdoor common garden experiment) to contrast inferred patterns of divergent selection based on phenotypes expressed in the greenhouse vs. in more natural environmental conditions (Methods S1).

| Ecophysiological measurements
We measured stem height and basal diameter of all seedlings at the time of the start of the watering treatments (1 month after transplantation) and the end of the experiment (5 months after transplantation). Absolute growth rate (AGR) and relative growth rate were highly correlated (r = .75, p < .001 and r = .65, p < .001, respectively), so only height growth rates were used for further analyses.
One month after the start of the watering treatments, we measured a set of physiological parameters related to leaf photochemical efficiency, carbon assimilation rates, photoprotection and transpiration on sun leaves of two to five seedlings per maternal family per treatment. For this, we measured chlorophyll a fluorescence using a portable pulse-modulated fluorimeter (FMS2, Hansatech Instruments). Minimal (F 0 ) and maximal (F m ) fluorescence were measured at predawn (n = 1088). Around solar noon (12-14 p.m.), we measured maximal (F m ′) and steady-state (F s ) fluorescence and applied a far-red pulse (740 nm) for 5 s to obtain minimum light-adapted fluorescence (F 0 ′) (n = 1149). For further analyses, we calculated which represents the excess energy dissipation via heat; the ef- represents the proportion of absorbed energy used in photochemistry; and maximum quantum yield of Photosystem II in light Note: Garden indicates the common garden experiment in which populations were assayed (GH = greenhouse experiment, OE = outdoor experiment). Elevation is given in metres a.s.l. n Fam is the number of open-pollinated maternal families in the common garden, and n Gen indicates the number of genotyped individuals. solar noon (12-14 p.m.) using a leaf porometer (SC-1; Decagon Devices) (n = 1300). Stomatal conductance was calculated on both area (g s,area ) and mass (g s,mass ) basis multiplying g s,area by the specific leaf area (see paragraph below).
Leaves used for measuring physiological parameters were subsequently collected for analyses of leaf morphology. Leaf thickness was measured using a micrometre when leaves were still fresh. Then, we scanned them to determine leaf lamina area and perimeter using WinFOLIA software (Regent Instruments Inc.).
We calculated the perimeter-to-area ratio, which combines information on leaf area and the degree of lobulation of the leaf. Finally, we oven-dried the leaves at 60°C for 48 h and weighed them to calculate specific leaf area (SLA), a trait related to carbon, nitrogen and phosphorous economics, by dividing fresh leaf area by leaf dry mass.

| Population genetic structure
We analysed patterns of SNP genetic structure for each species separately using the Bayesian Markov chain Monte Carlo (MCMC) clustering method implemented in the program structure version 2.3.3 (Pritchard et al., 2000;Falush et al., 2003;Hubisz et al., 2009) and discriminant analysis of principal components (DAPC; Jombart et al., 2010). Neither of the two clustering methods requires a priori population delimitation, but they differ in their analytical approaches and assumptions. structure builds clusters by minimizing Hardy-Weinberg and gametic disequilibrium within clusters and typically fails to detect isolation-by-distance (IBD). By contrast, the multivariate DAPC does not rely on the assumptions of structure and could be more efficient to detect complex patterns of genetic differentiation (Jombart et al., 2010).
For the Bayesian clustering analyses in structure, we considered correlated allele frequencies and an admixture model without prior information on population origin (Hubisz et al., 2009).
We performed 15 independent runs for each value of K (K = 1-8 for each species) with a burn-in period of 200,000 steps and a run length of 1,000,000 MCMC cycles. We retained the 10 runs having the highest likelihood for each value of K and estimated the bestsupported number of genetic clusters with the log probability of the data [Ln Pr (X|K)] (Pritchard et al., 2000) and the ΔK method (Evanno et al., 2005) as implemented in structure harvester (Earl & VonHoldt, 2012). For visualization, we used the "full search" algorithm in the program clumpp version 1.1.2 to align replicated runs and average individual assignment probabilities for the same K value (Jakobsson & Rosenberg, 2007). DAPCs were implemented using the package "adegenet" (Jombart & Ahmed 2011) in r version 3.6.1 (R Core Team, 2022).
We additionally estimated the magnitude of neutral genetic structure of each species by calculating overall and pairwise F ST values, using the estimator of Weir & Cockerham (1984) as implemented in the r package "hierfstat" (Goudet 2005

| Population divergence in ecophysiological traits and signals of selection
We used the approach developed by Ovaskainen et al. (2011) and Karhunen et al. (2013) with modifications by  to investigate if traits were potentially under amongpopulation divergent selection. This method appears to be prefer- For the analyses, we first estimated a co-ancestry matrix (i.e., genetic distances based on genomic markers) between all pairs of populations using SNP allele frequencies assuming an admixture F-model (AFM) and using a Metropolis-Hastings algorithm implemented in the r package "RAFM" (Karhunen & Ovaskainen 2012).
This model assumes that populations diverged from a common ancestral genetic pool and that genetic differences among populations have emerged from neutral drift. We conducted these analyses for each species and common garden, independently. We ran 30 independent Markov chains with a burn-in of 20,000 iterations followed by 10,000 iterations and a thinning interval of 10. To test for chain convergence, we estimated the mean of the off-diagonal values of F ST of each matrix, ran Heidelberger's test (Heidelberger & Welch, 1981) and calculated Geweke's z-scores (Geweke, 1991) using the "coda" r package. Chains that did not converge were excluded from further analyses.
Second, we evaluated the overall evidence of selection across all populations by comparing mean additive trait values for the study populations with those estimated for the assumed ancestral population. For this purpose, we used the MH and S-test functions of the "driftsel" r package (Ovaskainen et al., 2011;Karhunen et al., 2013).
MH conducts a Bayesian mixed-effects animal model that takes into account the co-ancestry matrix obtained from AFM (i.e., the genetic structure from SNP markers), the family structure among and  (Lewontin & Krakauer, 1973), which accounts for deviations in F ST due to demographic factors (Whitlock, 2008;Whitlock & Guillaume, 2009).

| Environmental drivers of population genetic divergence in ecophysiological traits
For traits showing at least weak evidence of being under divergent selection (estimated S > 0.5), we also implemented an H-test (sensu Karhunen et al., 2014) to identify the potential environmental drivers of population divergence, using the modification of this test developed by . Briefly, this test estimates a statistic (H* sensu    & Samani (1985). I m,summer was estimated as the difference between precipitation and potential evapotranspiration for July-September.
The environmental variables were centred and scaled before performing the PCA using the prcomp function in r.
We used the individual PC scores from the PCA-phenotype of Q. faginea to test for associations between multivariate phenotypes and climate using H*-tests and Pearson correlations.  Accordingly, several pairwise F ST estimates were significantly different from zero but consistently low in both Q. faginea (range = 0.002-0.073; Table S3) and Q. lusitanica (range = 0.006-0.089; Table S4).

| Genomic data
Analyses based on data sets excluding outlier loci (i.e., putatively under selection; see details in Methods S2) produced virtually identical results in structure (Figures S1 and S3) and DAPC ( Figure S4) and yielded similar F ST estimates (Tables S3 and S4).

| Signals of divergent selection on ecophysiological traits
Our results provided evidence of divergent selection in different leaf and growth traits. All traits for Q. faginea and seven out of 11 traits for Q. lusitanica exhibited differences among populations in at least one environment (Figure 2; Figure S5). Overall, more inland Q. faginea populations (MOL, QUI, SOT, SAL) exhibited lower specific leaf area and higher stomatal conductance (g s,mass , g s,area ) in the greenhouse experiment, particularly under dry conditions, and lower absolute growth rates under well-watered conditions ( Figure S6). In the outdoor common garden experiment, CUE, SOT and SAL had the lowest specific leaf area ( Figure S7) and dry (red dots) conditions. The 95% credible intervals were small and are hidden behind the dots. Mass-based stomatal conductance (g s,mass ) in the well-watered and dry conditions had the same signal of selection (S = 1). Asterisks (*) indicate at least one population pair with >95% posterior probability of having different additive population means, as estimated in well-watered conditions (blue) and dry conditions (red) in the greenhouse common garden. Traits include specific leaf area (SLA), leaf lamina area (area), leaf thickness (thickness), leaf perimeter-to-area ratio (PA), mass-based stomatal conductance (g s,mass ), area-based stomatal conductance (g s,area ), maximum quantum yield of photosystem II in light (F v '/F m ′), effective quantum yield of photosystem II (Φ PSII ), nonphotochemical quenching (NPQ), and absolute and relative height growth rate (AGR and RGR).

(a) (b)
in growth traits (AGR and RGR) in the well-watered treatment, but only for Q. faginea (Figure 2 and relative growth rates in the dry treatment and for most physiological traits in the well-watered treatment.

| Potential role of climate on divergent selection
The

| Weak genetic structure in Q. faginea and Q. lusitanica
Previous studies have shown strong genetic structure of European white oaks in maternally inherited cpDNA .

F I G U R E 3
Principal component analyses (PCAs) of 19 bioclimatic variables (data from www.world clim.org), soil pH (data from Trabucco & Zomer, 2010), and annual and summer index of moisture (I m,annual and I m,summer , respectively) for Quercus faginea and Q. lusitanica populations (left panels). Small grey circles and large coloured symbols indicate variable loadings and population scores for the first two PC axes (PC-climate1, PC-climate2), respectively. Population codes are described in Table 1. The colour scale used for populations reflects their latitudinal ranking following Figure 1. Circles indicate populations with genomic and phenotypic data, squares indicate populations only with phenotypic data and diamonds indicate populations only with genomic data. Right panels show the results of the H*-test (Csilléry, Ovaskainen, et al., 2020;Karhunen et al., 2014) and r-coefficients from Pearson correlations between traits and the first two climatic PC-axes for each species. H* values close to 1 (dark blue) indicate evidence for climatic adaptation (see text for details). Colour intensity is proportional to H* values and correlation coefficients, in a blue scale for positive values and in a red scale for negative values. Significant (p < .05) r-coefficients are underlined. Only traits showing a signal of selection greater than 0.5 were considered for these analyses. Trait abbreviations as in Figure 3. In the experiment column (Exp): Outdoor, outdoor common garden experiment; GH-WW, greenhouse experiment under well-watered treatment; and GH-D, greenhouse experiment under dry treatment.
In the Iberian Peninsula, Olalde et al. (2002) found two main maternal lineages differentiating northwestern from southeastern white oak populations, as well as two other minor lineages specific to northeastern populations. The existence of multiple Iberian glacial refugia during the last glacial maximum and posterior recolonization and hybridization between species have been proposed as the main processes driving the current genetic structure in cpDNA in white oaks (Olalde et al., 2002;. Using biparentally inherited nuclear markers, we did not observe the same pattern. In fact, structure and DAPC analyses revealed a very weak genetic structure in both Q. faginea and Q. lusitanica (Figure 1). Similar contrasting geographical patterns in nuclear versus plastidial markers have been commonly reported in other tree species, including oaks (e.g., Magri et al., 2007;Ortego et al., 2015;Pina-Martins et al., 2019).
Chloroplast and nuclear genomes have different forms of inheritance, mutation rates and vectors of dispersion (Petit et al., 2005). The chloroplast is maternally inherited in oaks and dispersed via seeds (Dumolin et al., 1995). Thus, the stronger genetic structure in chloroplast vs. nuclear loci probably reflects shorter dispersal distances of seeds by animals and gravity, relative to pollen dispersal distances by wind (e.g., Dow & Ashley, 1998;Cavender-Bares & Pahlich, 2009;Ortego et al., 2014;Sork et al., 2015). Greater genetic structure in chloroplast loci may also result from their greater sensitivity to genetic drift than nuclear markers, because of their haploid status and, thus, lower effective population sizes (Ortego et al., 2015).
The low population divergence in both Q. faginea and Q. lusitanica is consistent with that reported in other oak species from temperate regions (e.g., Cavender-Bares et al., 2011Pina-Martins et al., 2019). In wind-pollinated trees, only strong geographical isolation might effectively restrict gene flow and create marked neutral genetic divergence provided such isolation is maintained through extended periods of time (e.g., Cavender-Bares et al., 2011. Although the persistence and stability of oak populations from the Iberian Peninsula through the Pleistocene might have potentially contributed to the long-term isolation of local/regional gene pools, large effective population sizes are expected to have buffered genetic drift and resulted in a limited genetic structure across space and time (e.g., Lesser et al., 2013).
Thus, the interplay between the demographic history and life-history characteristics of oaks (i.e., long life, overlapping generations, wind pollination, outcrossing mating and large historical effective population size; Petit & Hampe, 2006) probably underlies the weak genetic structure observed in our two focal species.

| Consistent evidence of divergent selection on leaf morphology
Our results revealed that SLA was the studied trait that exhibited the strongest evidence of divergent selection, for any combination F I G U R E 4 Principal component analysis (PCA) of phenotypic traits measured in the greenhouse experiment under dry (red) and wellwatered (blue) treatments for the additive mean populations and the estimated ancestral population of Quercus faginea (left panel). Small grey circles and large blue/red circles indicate phenotypic trait loadings and population scores for the first two PC axes (PC-phenotype1, PC-phenotype2), respectively. Population codes are described in Table 1. "A" indicates the ancestral population mean estimated following Ovaskainen et al. (2011) and. Phenotypic traits: Specific leaf area (SLA), leaf lamina area (area), leaf thickness (thickness), leaf perimeter-to-area ratio (PA), mass-based stomatal conductance (g s,mass ), area-based stomatal conductance (g s,area ), maximum quantum yield of photosystem II in light (F v '/F m ′), effective quantum yield of photosystem II (Φ PSII ), nonphotochemical quenching (NPQ), and relative and absolute height growth rate (RGR and AGR). The right panel shows linear regressions between the second axis of the PCA of environmental variables (PC-climate2) and the second axis of the PCA of phenotypic traits (PC-phenotype2) for Q. faginea populations. Dots indicate additive population means for PC-phenotype2 in the dry (red) and well-watered (blue) treatments in the greenhouse experiment. H* values were obtained following , and R 2 coefficients from Pearson correlations are shown for each treatment. Asterisks for R 2 coefficients indicate p < .001.
of species, experimental treatment and garden. In particular, we detected genetically-based trait differentiation among populations, a signal of selection (S) close to 1, Q ST significantly higher than F ST and a strong association with climate (Figures 2 and 3; Figure S5, Table S5).
The importance of SLA in ecological processes and adaptation to resource gradients and drought tolerance has been demonstrated in numerous studies. SLA is part of the so-called Leaf Economics Spectrum, a suite of correlated leaf traits that define the strategies of investment and return of carbon, nitrogen and phosphorous at global scale (Wright et al., 2004;Shipley et al., 2006). In oaks, SLA has been found to be related to growth rates and other fitness components at different organizational scales, supporting its role as a key functional trait (Ramírez-Valiente et al., 2010. From an evolutionary perspective, leaf morphological traits, including SLA, have been shown to be evolutionarily labile, evolving in response to climatic gradients (Ramírez-Valiente et al., 2020;Hipp et al., 2020;Sancho-Knapik et al., 2021).
There was evidence of divergent selection in other morphological traits across species and common garden experiments, such as leaf lamina area, leaf thickness and perimeter-to-area ratio, although the signal of selection was lower than that detected for SLA, particularly for Q. lusitanica (Figure 2). Multivariate analyses for Q. faginea also revealed strong correlated selection on leaf morphological traits associated with the climate of the provenances (Figure 4).
Similar results have been reported in other oaks, where SLA but also leaf lamina area or leaf anatomical traits such as size and density of stomata appear to have evolved under divergent selection . Support for the adaptive value of SLA and other leaf morphological traits has been found in different plant species. For example, a study on Pinus canariensis found that population genetic differentiation in leaf mass per area (LMA, a trait inversely related to SLA) was higher than expectations based on neutral evolution (López et al., 2013). Keller et al. (2011)  For both Q. faginea and Q. lusitanica, mass-based stomatal conductance (g s,mass ) also presented a strong signal of divergent selection, although with contrasting species patterns ( Figure 2; Table S5).
In Q. faginea, the highest values of g s,mass under dry conditions were found for populations from dry and continental areas whereas, in Q. lusitanica, the highest values both in dry and well-watered conditions were found in the wettest population, located in the northernmost extreme of the distribution range of the species (Figure 3).
Stomatal control is a key response to water availability conditions.
Closing stomata under water stress allows plants to maintain plant water potentials but at the cost of reducing carbon assimilation (Medrano et al., 2002;Flexas et al., 2004). Under well-watered conditions, growth traits (AGR and RGR) also exhibited strong evidence of divergent selection for Q. faginea.
Populations from rainier and milder climates exhibited higher RGR ( Figure 4). It has been hypothesized that larger growth rates are beneficial under mesic environments because they increase the ability to compete for other resources such as light. Previous studies have demonstrated evidence of divergent selection on growth traits for many species (e.g., Keller et al., 2011;Csilléry, Ovaskainen, et al., 2020;Marin et al., 2020;Gauzere et al., 2020) including oaks (e.g., . Growth is considered a trait that integrates responses from different plant functions and is one of the components of plant performance together with reproduction and survival (Violle et al., 2007).

| Differences between species, common gardens and watering treatments in patterns of selection
In the greenhouse experiment, where Q. faginea and Q. lusitanica plants were grown under the same environmental conditions, the two species exhibited a contrasting number of traits with evidence of selection (S) and associations with climate (Figures 2 and 3).
Specifically, for Q. faginea, seven out of 11 traits exhibited a signal of selection higher than 0.9 and eight of them presented associations with climate (PC-climate2). In contrast, for Q. lusitanica, only SLA and g s,mass had strong evidence of selection (S > 0.9). Similar differences between species were found based on Q ST -F ST analyses (Table S5).
The weaker evidence of divergent selection and weaker associations between trait population divergence and climate in Q. lusitanica could be related to the fact that this species occupies a niche with lower cli-

| CON CLUS IONS
Our study demonstrates the importance of integrating genomic data and phenotypic information from common garden experiments to understand patterns of adaptation to climate. It should be borne in mind that although the method developed by Ovaskainen et al. (2011) and Karhunen et al. (2013) addresses some limitations of the classical Q ST -F ST analysis, it still assumes a normal distribution of traits, low mutation rates and additive genetic effects. Although our study used normally distributed traits and SNP markers, which might help in approximating the two first assumptions, future developments should account for additional sources of variation, such as dominance, G × E interactions and non-normally distributed traits (Karhunen et al., 2013). Collectively, our results support an important role of leaf morphology, stomatal conductance and growth rates to adapt to climate heterogeneity in two Mediterranean species. Lower levels of functional divergence in the narrowly distributed shrub Quercus lusitanica compared to Q. faginea suggest that climatic niche and life history traits might be key aspects in shaping the genetic structure in genes involved in plant functions. These results highlight the need for more multidisciplinary studies that integrate a thorough functional characterization of traits important to drought tolerance and genomic data to understand population responses to climate change. for physiological measurements. We also thank two anonymous reviewers for their constructive and valuable comments on an earlier version of the manuscript.

CO N FLI C T O F I NTE R E S T
The authors have no conflicts of interest to declare.