Nanying Liang

Laurent Bougrain

Averaging techniques for single-trial analysis of oddball event-related potentials

More and more effort is done in BCI research to improve its usability for patients, with respect to its communication speed and transmission accuracy. In this contribution, we experiment with BCI speller based on P300 evoked potential. More precisely, the typical form of event-related potential (ERP) inspires us to devise classification methods based on the similarity/dissimilarity in the time domain between single trials and one or several estimated ERP templates derived from subject recordings. The reliable estimation of template is difficult in a single trial due to the low signal-to-noise ratio (SNR) of electroencephalographic (EEG) signals. We first explicitly estimate the template using several averaging techniques: pointto-point averaging, cross-correlation alignment and dynamic time warping. Then we inexplicitly estimate several ERP templates using learning vector quantization algorithm combined with an extreme learning machine. Finally classification is realized based on the similarity/dissimilarity between the single trials and the template. Simulation is carried out using a BCI competition III data set acquired with the P300 speller paradigm. The experiments show that template-based classifiers can also obtain high accuracy.

Introduction

Brain-computer interface (BCI) system is a potentially powerful new communication and control option for those with severe motor disabilities. BCI system translates brain activity into commands for a computer or other devices. Electroencephalography (EEG) is the most studied potential noninvasive interface mainly due to its fine temporal resolution, ease of use, portability and low setup cost. Unfortunately, non-invasive implants produce a noisy signal because the skull dampens signals. Another substantial barrier to use EEG as a brain-computer interface is the extensive training required before users can work the technology.

Oddball paradigms are used by brain-computer interfaces to generate event-related potential (ERP) on targets selected by the user. The well-known P300 speller is based on this principle [START_REF] Farwell | Talking off the top of your head: toward a mental prosthesis utilizing event-related brain potentials[END_REF]. The main problem is to be able to detect ERP in a noisy electroencephalographic signal recorded on human scalp. Literature proposes a set of methods to detect ERPs based on averaging. Linear or nonlinear alignment is used to deal with the variable latencies in ERPs. A good alignment allows to obtain a better averaged ERP response, named "ERP template" here, and then to be able to detect ERP faster.

Section 2 of this paper presents a set of methods of template-based classifiers. Section 3 presents the results of these techniques on a well-known third BCI2000 competition data set acquired from a P300 speller. In section 4, we conduct the conclusion and give recommendation for future works.

Methods

Methods design for ERP detection use features extracted from the responses in the time domain. To extract features from the frequency domain after any transformation is difficult in this situation because ERPs are short-time events with a local peak. Here, we discuss two types of methods using one or several templates estimated explicitly and inexplicitly in the time domain respectively.

Classifiers based one ERP template

In oddball paradigm, one deviant rare stimulus is presented to subjects amongst a train of standard ones. Only the deviant stimulus induces an ERP. The most common evoked potential used in BCI is the P300. Its main characteristic is a positif latency that occurs around 300ms after the stimulus presentation. This leads us to derive a set of classifiers. All of them first explicitly estimate one ERP template by averaging the ERP responses in the time domain, and then use the distance between the response and the ERP template as the discriminant criterion.

point-to-point averaging (P2P) classifier: P2P classifier discriminates the ERP responses from non-ERP responses in a straightforward way: first it estimates one ERP template by simply averaging the measurements of ERP responses point-to-point, and then calculates the Euclidean distance between responses and the ERP template. For the P300 speller, the ERP response corresponds with the one producing the minimum distance within a set of column or row responses. Then the other responses are considered as non-ERP responses.

P2P classifier is simple for implementation. However the latency of the components of ERP varies with the external factors (such as target-to-target interval). Simply averaging the measurements point-to-point could blur the component. We then seek the solution of alignment before averaging. Such alignment includes linear alignment (such as, cross-correlation alignment) and nonlinear alignment (such as, dynamic time warping). cross-correlation averaging (CC) classifier: Linear alignments can be done by shifting the measurements before the point-to-point averaging. One of the well-known linear alignments is cross-correlation. The best shift of one sequence in order to align it to the other ones is given by the one producing the highest correlation [START_REF] Woody | Characterization of an adaptive filter for the analysis of variable latency neuroelectric signals[END_REF][START_REF] Wastell | Statistical detection of individual evoked responses: an evaluation of woody's adaptative filter[END_REF]. CC classifier first searches the best shift for pairs of ERP responses before averaging them point-to-point to obtain one template. Then it calculates the cross-correlation between the testing response and the ERP template. The predicted ERP response corresponds with the one giving the maximum correlation value with the template. Such scheme accommodates variable linear latency. However, if one looks the ERP measurements more closely, it is found that the latency of the ERP component is nonlinear. In this case, the nonlinear alignment methods are better than the linear ones for template estimation. dynamic time warping averaging (DTW) classifier: Dynamic Time warping distance is widely used in speech domain. It is able to deal with compression/expansion by comparing the distance of each point of the first sequence to every point of the second one using the information of the amplitude and/or the first derivative of amplitude of the signals [START_REF] Keogh | Derivative dynamic time warping[END_REF]. Hence, DTW classifier first finds out the best nonlinear alignment path giving the closest distance, and then averages the ERP responses according to the alignment path to obtain one ERP template [START_REF] Picton | Evaluation of brain-stem auditory evoked potentials using dynamic time warping[END_REF][START_REF] Wang | Warp-averaging event-related potentials[END_REF][START_REF] Casarotto | Dynamic time warping in the analysis of event-related potentials[END_REF]. Finally, it recognizes the ERP responses as the one producing the minimum DTW distance.

Methods with several templates of ERP & non-ERP

On one hand, considering the complexity of ERP responses, one template could be insufficient to capture the variability of ERP components. On the other hand, one would like to introduce non-ERP templates into decision making process done by the classifiers. Clustering methods can produce several templates. The aim of clustering methods is to split a set of patterns into clusters as homogeneous as possible, i.e., patterns are gathered by similarity. A classical rule to assign a new pattern to a cluster is based on distances between the input pattern and templates. LVQ is the most famous algorithm in this category.

learning vector quantization (LVQ):

The first version of the learning vector quantization algorithm proposed by Kohonen [START_REF] Kohonen | The self-organizing map[END_REF] is a supervised clustering method designed for classification. Input weights (IW) of connections linking inputs to each neuron is a vector m with the same dimension than input x. Each neuron is assigned to a class according to the pre-defined layer weights (LW). A classic competitive step is used to determine which neuron is the closest one to the current input with the formula: c = arg min i {||xm i ||}, and only the winner neuron is updated according to the following rules:

m c (t + 1) = m c (t) + α(t)[x(t) -m c (t)] if x and m c belong to the same class m c (t + 1) = m c (t) -α(t)[x(t) -m c (t)] if
x and m c belong to different classes where α is learning rate. To apply this algorithm, one needs to decide the value of learning rate and the number of neurons assigned to each class. Further considering the multichannel recording of EEG measurements, if one simply merges features from each channel, the dimension of IW increases with the number of channels, and it thus causes memory problem in simulation. We then consider multichannel LVQ. multichannel learning vector quantization (mLVQ): For multichannel LVQs, one LVQ model is created for each channel, and trained according to the above rules to update IWs. Thus, there are several ERP and non-ERP templates generated for each channels. In order to combine the templates in an optimal way, we further introduce an update scheme for LWs: LW = (-||X -IW ||) † T , where X is the input matrix and T is the target matrix. The symbol † represents the pseudo-inverse operator. Such a solution of LW is known as minimum norm least square solution. It is originally found being used in extreme learning machine (ELM) in [START_REF] Huang | Extreme learning machine: A new learning scheme of feedforward neural networks[END_REF].

Results

Wadsworth BCI data set

The results presented here are based on the P300 speller data set from the BCI competition III [START_REF] Krusienski | Wadsworth BCI Dataset (P300 Evoked Potentials)[END_REF] 1 . The data set contains training and testing EEG data from two subjects. There are 85 letters for training and 100 letters for testing. For each letter, the recording consists of 15 epochs, and within each epoch, there are 12 flashings. For each epoch, a random permutation is chosen to highlight rows and columns. There is a 6x6 grid in this application containing 26 letters, 9 digits and one dash character. We are interested in the measurement window starting from the onset of the flashing and till one second, which consists of 240 samples at a sampling rate of 240Hz. In case the row/column is highlighted exactly where the target letter is located, we got the ERP response; otherwise the non-ERP response.

Experimental Results

ERP Templates

Figure 1 shows the comparison between the ERP templates estimated by P2P averaging, CC and DTW alignment averaging, respectively. Templates are built using the first 128 trials of the training set of subject A. Indeed, the algorithms using CC and DTW prefer to pair trials in a binary tree structure. The latency of the P300 component in the estimated templates by CC and DTW alignment averaging methods is different from P2P averaging. This is the effect of aligning the component in the time domain before averaging. And we found that the template estimated by DTW was smoother than by cross-correlation.

Classifier Comparison

Simulation is carried out to compare the performance of the template-based classifiers using P2P, CC, DTW and LVQ. The result of linear discriminant analysis (LDA) is also provided for its simplicity and good performance in the application of BCIs. Table 1 shows the testing accuracy for subjects A and B. The results are based on the raw measurements from channel Cz (i.e., no pre-processing is involved). We used the full training set to estimate the ERP templates for P2P, CC and DTW; while for LVQ and LDA, they are used to estimate the parameters of the models. The full testing set is used to evaluate the performance of the classifiers. From the table, we can see that LVQ and LDA achieve similar performance and are better than the classifiers using P2P, CC, and DTW.

Table 2 shows the testing accuracy for subjects A and B based on measurements of all channels. Due to limit of memory and the large size of the input data, one needs to re-sample the signals. Thus, for Table 2, we first smooth the signal using a moving averaging filter with window size of 13, and then re-sample it with a down-sampling factor of 13. Further, multichannel version of LVQ is used here. We find that the performance of mLVQ and LDA are similar again and both increase as more channels are used for classification; while the performance of the explicit template-based classifiers using P2P, CC and DTW is bad.

From the above simulation results, we found that the performance of classifiers (such as LVQ), which take account of both ERP and non-ERP templates is better than those (such as P2P, CC, and DTW) which use only one ERP template. The optimization of the combination of the templates from different channels through assigning different weights is even more critical for improving the accuracy. take into account the specificities of non-ERP responses to catch small differences between noisy ERP and non-ERP responses. The first difficulty exists for any classification method. But, our method mLVQ obtains similar results as LDA which is a reference technique for the P3Speller. It is because mLVQ takes into account the second remark. Thus, template-based classifiers can also obtain good results. Moreover, this approach can extract knowledge about the ERP waveform of patients.

In this article, we also evaluated the promising dynamic time warping distance. This kind of similarity measure is widely used in speech domain and align efficiently two sequences even if they differ each other from compression and/or dilatation on some segments. However, from our opinion, DTW suffers from the low SNR. The real EEG signals used in our simulation are too noisy. Many artificial peaks appear in the recording and thus it is difficult to obtain a strong alignment. In order to mitigate this difficulty, one may reduce the window of investigation to [250ms, 650ms] and constrain the window size of permitted distortion to short time. It should help to focus the attention on the critical zone.

In practice, it is also very interesting to reduce the number of channels in order to equip the patient with a compact system with only several necessary electrodes. Thus, channel selection will be investigated in a further work.

Figure 1 :

 1 Figure 1: Comparison of ERP templates (red solid curves) and averaged background EEG activities (black dot curves) for responses recorded on channel Cz using training data from subject A. From top to bottom, templates are estimated by point-to-point averaging, cross-correlation and dynamic time warping alignment averaging, respectively.

4

  Conclusion and further workBCI applications based on P300 ERP detection, such as the P3speller, need statistical or machine learning techniques with temporal features. ¿From a recording point of view, ERPs are graphicelements. They have a specific waveform clearly different from EEG background activity. Thus, efficient modeling of the waveform can help to detect ERPs. Averaging techniques and clustering methods allow extracting one or several ERP templates. Averaging techniques using point-topoint averaging, cross-correlation and dynamic time warping are not efficient. They produce only one ERP template and suffer from two difficulties: first, responses are too noisy to easily distinguish ERP from non-ERP responses; second, according to the previous remark, they do not

					Method
	Subject Epoch P2P CC DTW LVQ LDA
		05	7%	8%	15%	17% 18%
	A	10	5%	26%	28%	28% 30%
		15	3%	21%	39%	43% 41%
		05	2%	6%	6%	8%	15%
	B	10	5%	9%	10%	13% 19%
		15	3%	16%	15%	21% 26%

Table 1 :

 1 Experimental results based on channel Cz for subjects A and B. The percentage of correct classified letters contained in the testing set is shown regarding using the first 5 epochs, 10 epochs and all epochs.

					Method
	Subject Epoch P2P CC DTW mLVQ LDA
		05	9%	7%	7%	47%	45%
	A	10	13% 6%	12%	77%	78%
		15	14% 7%	15%	87%	88%
		05	3%	2%	4%	72%	76%
	B	10	3%	3%	3%	91%	92%
		15	4%	4%	3%	96%	96%

Table 2 :

 2 Experimental results based on all channels for subjects A and B. The percentage of correct classified letters contained in the testing set is shown regarding using the first 5 epochs, 10 epochs and all epochs.

The data set is available at: http://ida.first.fraunhofer.de/projects/bci/competition_iii/


Dragos Postolache

Sylvie Oddou-Muratorio
email: sylvie.muratorio@inrae.fr

Elia Vajana

Francesca Bagnoli

Erwan Guichoux

Arndt Hampe

Gregoire Le Provost

Isabelle Lesur

Flaviu Popescu

Ivan Scotti

Postolache D A1

Le Provost

Genetic signatures of divergent selection in European beech ( Fagus sylvatica L.) are associated with the variation in temperature and precipitation across its distribution range

Keywords: candidate gene, phenology, drought stress, divergence outlier, genotypeenvironment associations analyses, forest tree, local adaptation

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Local adaptation is pervasive in forest tree populations [START_REF] Alberto | Potential for evolutionary responses to climate change -evidence from tree populations[END_REF][START_REF] Savolainen | Gene Flow and Local Adaptation in Trees[END_REF] and is expected to play a major role in their response to ongoing environmental changes [START_REF] Fady | Evolution-based approach needed for the conservation and silviculture of peripheral forest tree populations[END_REF]. Local adaptation implies that some key adaptive traits are genetically differentiated among populations, and thus that individual populations could evolve differently to the same environmental stress due to their different genetic setup [START_REF] Kawecki | Conceptual issues in local adaptation[END_REF]. Most temperate tree species have developed their present-day geographical patterns of local adaptation following a considerable range expansion from their glacial refugia after the Last Glacial Maximum (19.5-26 kyr BP, [START_REF] Clark | The Last Glacial Maximum[END_REF][START_REF] De Lafontaine | Invoking adaptation to decipher the genetic legacy of past climate change[END_REF]. There is ample concern that population processes going along with this postglacial range expansion, such as founder events, genetic drift, or allele surfing might have left lasting imprints that could compromise the correct identification of adaptive genetic variation in extant natural tree populations (de Villemereuil, [START_REF] De Villemereuil | Genome scan methods against more complex models: When and how much should we trust them?[END_REF][START_REF] Hoban | Finding the Genomic Basis of Local Adaptation: Pitfalls, Practical Solutions, and Future Directions[END_REF][START_REF] Rellstab | A practical guide to environmental association analysis in landscape genomics[END_REF]. Yet very few surveys of adaptive genetic variation in forest trees have to date assembled two prerequisites to properly account for species' postglacial population dynamics: (i) a rangewide perspective with a sampling that thoroughly replicates populations from different regions and (ii) independent paleoecological evidence documenting the species' postglacial range dynamics and expansion [START_REF] De Lafontaine | Invoking adaptation to decipher the genetic legacy of past climate change[END_REF]. Empirical research to elucidate the issue is urgently needed because knowledge on local adaptation is crucial for conceiving conservation and management practices to adapt forest tree species and their ecosystems to ongoing environmental change [START_REF] Aitken | Assisted Gene Flow to Facilitate Local Adaptation to Climate Change[END_REF][START_REF] Fady | Evolution-based approach needed for the conservation and silviculture of peripheral forest tree populations[END_REF][START_REF] Oney | Intraspecific variation buffers projected climate change impacts on Pinus contorta[END_REF].

Reciprocal transplant experiments have been the classic approach to investigate local adaptation in forest trees, highlighting that phenotypes generally match their environment in extant populations [START_REF] Rehfeldt | Intraspecific responses to climate in Pinus sylvestris[END_REF][START_REF] Savolainen | Gene Flow and Local Adaptation in Trees[END_REF]. Provenance trials have also shown that trees generally have high levels of phenotypic plasticity at adaptive traits, and high levels of genetic variability within populations [START_REF] Alberto | Potential for evolutionary responses to climate change -evidence from tree populations[END_REF][START_REF] Gárate-Escamilla | Range-wide variation in local adaptation and phenotypic plasticity of fitness-related traits in Fagus sylvatica and their implications under climate change[END_REF]. This combination of local adaptation (i.e., mean trait values close to the optimum) and high within-population variation at key adaptive traits (i.e., large variance around the means) indicates a high genetic load, which in turn can be an asset when facing a swift environmental change [START_REF] Savolainen | Gene Flow and Local Adaptation in Trees[END_REF]. The development of population genomics has provided complementary approaches to study local adaptation [START_REF] Hoban | Finding the Genomic Basis of Local Adaptation: Pitfalls, Practical Solutions, and Future Directions[END_REF][START_REF] Lind | The genomics of local adaptation in trees: Are we out of the woods yet?[END_REF]. Two approaches, in particular, have become widely used to identify loci involved in local adaptation for non-model organisms: differentiation outlier analyses, which aim at identifying loci with disproportionate allele frequency differentiation among populations, and gene-environment association (GEA) analyses, which aim at identifying loci exhibiting significant correlations with ecological variables. A key strength of such approaches is that they are cost-efficient for targeting large numbers of populations across well-defined environmental gradients, and therefore for investigating ecological hypotheses on the genomic basis of local adaptation [START_REF] Capblancq | Climate associated genetic variation in Fagus sylvatica and potential responses to climate change in the French Alps[END_REF][START_REF] Rellstab | Signatures of local adaptation in candidate genes of oaks ( Quercus spp.) with respect to present and future climatic conditions[END_REF][START_REF] Temunović | Candidate gene SNP variation in floodplain populations of pedunculate oak ( Quercus robur L.) near the species' southern range margin: Weak differentiation yet distinct associations with water availability[END_REF].

Present-day patterns of adaptive phenotypic and genetic differentiation have built up relatively quickly. At most, it is only a few hundreds of generations ago that most temperate forest tree species recolonized large areas becoming available after the Last Glacial Maximum.

Studies combining extensive surveys of fossil records (pollen and macro-remains) and of population genetic variation have been able to provide detailed direct evidence of the number and spatial location of glacial refugia, the trajectory of postglacial expansion routes, the existence/location of hybrid zones, and the timing of expansion events (e.g., de Lafontaine, Amasifuen [START_REF] De Lafontaine | Cryptic no more: Soil macrofossils uncover Pleistocene forest microrefugia within a periglacial desert[END_REF][START_REF] Magri | A new scenario for the quaternary history of European beech populations: Palaeobotanical evidence and genetic consequences[END_REF]. This knowledge represents a highly valuable baseline information for disentangling demographic effects from those of selection, which is a major challenge for differentiation outlier and GEA analyses [START_REF] De Villemereuil | Genome scan methods against more complex models: When and how much should we trust them?[END_REF][START_REF] Frichot | Detecting adaptive evolution based on association with ecological gradients: Orientation matters![END_REF][START_REF] Hoban | Finding the Genomic Basis of Local Adaptation: Pitfalls, Practical Solutions, and Future Directions[END_REF][START_REF] Rellstab | A practical guide to environmental association analysis in landscape genomics[END_REF]. Theoretical studies have shown that population expansions into new areas can go along with repeated founder effects, increasing random fluctuations of allele frequencies, possibly leading to the rise of neutral mutations to high frequencies ("allele surfing"), loss of genetic diversity and strong spatial genetic structure (SGS) along the expansion axis (de Lafontaine, Ducousso, [START_REF] De Lafontaine | Stronger spatial genetic structure in recolonized areas than in refugia in the European beech[END_REF][START_REF] Excoffier | Genetic Consequences of Range Expansions[END_REF][START_REF] Slatkin | Isolation by Distance in Equilibrium and Non-Equilibrium Populations[END_REF]. Allele surfing, in particular, could be mistaken for the increase in allelic frequency of a beneficial mutation propagated by selection [START_REF] Paulose | The impact of long-range dispersal on gene surfing[END_REF][START_REF] Ruiz Daniels | Inferring selection in instances of long-range colonization: The Aleppo pine ( Pinus halepensis ) in the Mediterranean Basin[END_REF]. As both allele surfing and selection typically affect only a subset of loci in the genome, the former must be carefully considered when screening for the latter in expanding populations. Alternatively, colonization by many individuals should result in high genetic diversity at the colonization front and shallow SGS, particularly when founders originate from a variety of source populations. The effective number of founders depends on patterns of long-distance dispersal and on a variety of demographic processes and life-history traits [START_REF] Austerlitz | Effects of Colonization Processes on Genetic Diversity: Differences Between Annual Plants and Tree Species[END_REF][START_REF] Fayard | Long distance dispersal and the fate of a gene from the colonization front[END_REF][START_REF] Roques | Allee effect promotes diversity in traveling waves of colonization[END_REF]. The actual relevance of these processes during the postglacial range expansion of temperate forest trees and their possible traces in the presentday population genetic structures remain, however, under investigated.

This study takes advantage of the outstandingly well-known postglacial population history of the European beech (Fagus sylvatica L.) to investigate the climate-associated genetic variation across its distribution range. Beech putative glacial refugia and colonization routes have been identified based on very detailed pollen records and genetic population surveys with chloroplast and isozyme markers [START_REF] Magri | A new scenario for the quaternary history of European beech populations: Palaeobotanical evidence and genetic consequences[END_REF]. Beech is known to be highly sensitive to summer droughts [START_REF] Aranda | Variation in photosynthetic performance and hydraulic architecture across European beech (Fagus sylvatica L.) populations supports the case for local adaptation to water stress[END_REF][START_REF] Knutzen | Recent climate warming-related growth decline impairs European beech in the center of its distribution range[END_REF], and, to a lesser extent, to late frosts [START_REF] Kreyling | Local adaptations to frost in marginal and central populations of the dominant forest tree Fagus sylvatica L. as affected by temperature and extreme drought in common garden experiments[END_REF][START_REF] Petit-Cailleux | Comparing statistical and mechanistic models to identify the drivers of mortality within a rear-edge beech population. BioRxiv 645747 Ver. 7[END_REF]. Genetic variation has been investigated at various climate-related phenological traits (Gárate-Escamilla, Hampe, [START_REF] Gárate-Escamilla | Range-wide variation in local adaptation and phenotypic plasticity of fitness-related traits in Fagus sylvatica and their implications under climate change[END_REF][START_REF] Gauzere | Microgeographic adaptation and the effect of pollen flow on the adaptive potential of a temperate tree species[END_REF][START_REF] Gömöry | Trade-off between height growth and spring flushing in common beech (Fagus sylvatica L.)[END_REF][START_REF] Kramer | Chilling and forcing requirements for foliage bud burst of European beech ( Fagus sylvatica L.) differ between provenances and are phenotypically plastic[END_REF][START_REF] Vitasse | Altitudinal differentiation in growth and phenology among populations of temperate-zone tree species growing in a common garden[END_REF], physiological or morphological traits [START_REF] Bresson | To what extent is altitudinal variation of functional traits driven by genetic adaptation in European oak and beech?[END_REF][START_REF] Hajek | Intraspecific variation in wood anatomical, hydraulic, and foliar traits in ten European beech provenances differing in growth yield[END_REF][START_REF] Wortemann | Genotypic variability and phenotypic plasticity of cavitation resistance in Fagus sylvatica L. across Europe[END_REF] and performance traits [START_REF] Gárate-Escamilla | Range-wide variation in local adaptation and phenotypic plasticity of fitness-related traits in Fagus sylvatica and their implications under climate change[END_REF]. Phenological traits such as the timing of budburst and leaf senescence show consistent patterns of genetic variation across latitude or elevation at various spatial scales, with populations of higher elevation or latitude flushing earlier than populations from low elevation or latitude in common garden conditions [START_REF] Gauzere | Microgeographic adaptation and the effect of pollen flow on the adaptive potential of a temperate tree species[END_REF][START_REF] Gömöry | Trade-off between height growth and spring flushing in common beech (Fagus sylvatica L.)[END_REF][START_REF] Vitasse | Altitudinal differentiation in growth and phenology among populations of temperate-zone tree species growing in a common garden[END_REF]. Genetic variation for performance traits, such as growth and juvenile survival, also shows a spatial structure, driven by spatial variations of maximal potential evapotranspiration [START_REF] Gárate-Escamilla | Range-wide variation in local adaptation and phenotypic plasticity of fitness-related traits in Fagus sylvatica and their implications under climate change[END_REF]. By contrast, other functional traits involved in photosynthesis and transpiration are usually only weakly differentiated among populations but show instead high withinpopulation variation [START_REF] Hajek | Intraspecific variation in wood anatomical, hydraulic, and foliar traits in ten European beech provenances differing in growth yield[END_REF]. These contrasting patterns of differentiation raise questions about the role of phenological and physiological traits in local adaptation and the spatial scales at which it takes place. Only a few published studies have so far used genetic approaches to investigate local adaptation, and mostly at local to regional scales [START_REF] Capblancq | Climate associated genetic variation in Fagus sylvatica and potential responses to climate change in the French Alps[END_REF][START_REF] Csilléry | Detecting short spatial scale local adaptation and epistatic selection in climate-related candidate genes in European beech ( Fagus sylvatica ) populations[END_REF][START_REF] Cuervo-Alarcon | A candidate gene association analysis identifies SNPs potentially involved in drought tolerance in European beech (Fagus sylvatica L.)[END_REF][START_REF] Krajmerová | Nucleotide polymorphisms associated with climate, phenology and physiological traits in European beech (Fagus sylvatica L.)[END_REF][START_REF] Lalagüe | Nucleotide diversity and linkage disequilibrium at 58 stress response and phenology candidate genes in a European beech (Fagus sylvatica L.) population from southeastern France[END_REF][START_REF] Müller | A candidate gene-based association study reveals SNPs significantly associated with bud burst in European beech (Fagus sylvatica L.)[END_REF][START_REF] Pluess | Genomeenvironment association study suggests local adaptation to climate at the regional scale in Fagus sylvatica[END_REF].

This study investigates range-wide patterns of adaptive genetic variation in beech using 405 SNPs that are located in candidate genes involved either in budburst phenology and dormancy regulation [START_REF] Lalagüe | Nucleotide diversity and linkage disequilibrium at 58 stress response and phenology candidate genes in a European beech (Fagus sylvatica L.) population from southeastern France[END_REF][START_REF] Lesur | A unigene set for European beech ( Fagus sylvatica L.) and its use to decipher the molecular mechanisms involved in dormancy regulation[END_REF] or in response to stresses [START_REF] Lalagüe | Nucleotide diversity and linkage disequilibrium at 58 stress response and phenology candidate genes in a European beech (Fagus sylvatica L.) population from southeastern France[END_REF]. We genotyped 446 individuals from 64 populations covering the entire species range (Fig. 1) to address the three following questions: (Q1) What are the risks that past population demography, including post-glacial recolonization, blur potential selective imprints on genetic structure? We expect limited genetic drift and allele surfing in beech. (Q2) Do certain loci show imprints of local adaptation? We expect genes related to phenological traits to show stronger signals of selection compared to genes related to stress-response traits. (Q3) if the spatial and climatic effects can be separated, what are the respective impacts of temperature-versus precipitation-related variables on adaptive differentiation? In line with Q2, we expected the temperature variables to stand out more than precipitation variables.

Materials and Methods

Sampling design

European beech is a dominant broadleaved tree species of many lowland and mountain forests across Europe, extending from Spain to the Carpathians and from Sicily to southern Sweden. We sampled 446 adult trees in 64 populations across Europe (Fig. 1A), thoroughly covering the geographical range and bioclimatic niche of beech (Fig. 1B) and including all the major glacial refugia identified by [START_REF] Magri | A new scenario for the quaternary history of European beech populations: Palaeobotanical evidence and genetic consequences[END_REF]. Leaves were collected from 4 to 10 (7 on average) haphazardly chosen dominant adult trees (at a minimal distance of 40 m from each other) growing in native beech stands.

SNP development, genotyping and filtering

Nuclear DNA was extracted from 20-30 mg of dry leaf tissue per individual with the DNeasy Plant Mini Kit (QIAGEN) following the manufacturer's instructions. DNA concentration was measured on a ND-8000 NanoDrop spectrophotometer (Thermo Scientific, Wilmington, USA). Samples were genotyped at 405 SNPs distributed in two multiplex assays.

The first assay of 165 SNPs previously developed by [START_REF] Lalagüe | Nucleotide diversity and linkage disequilibrium at 58 stress response and phenology candidate genes in a European beech (Fagus sylvatica L.) population from southeastern France[END_REF] was carried out using Kompetitive Allele Specific PCR (KASP; [START_REF] He | SNP genotyping: The KASP assay[END_REF]. Among those, 37 SNPs were located in 15 genes annotated as phenology-related in Quercus petraea. The other 128 SNPs were located in 37 stress-related genes that were selected from different sources (see [START_REF] Lalagüe | Nucleotide diversity and linkage disequilibrium at 58 stress response and phenology candidate genes in a European beech (Fagus sylvatica L.) population from southeastern France[END_REF] for details): (1) literature on candidate genes involved in plant response to abiotic stress; (2) sequences of three proteins involved in cavitation resistance in beech); (3) annotated amplicons from Q. petraea and Q. robur.

The second assay targeted 240 SNPs in six multiplexes of 40 SNPs each. This assay was developed for this study from available genomic resources [START_REF] Lesur | A unigene set for European beech ( Fagus sylvatica L.) and its use to decipher the molecular mechanisms involved in dormancy regulation[END_REF] and included 104 SNPs located in 51 genes differentially expressed in quiescent buds (QB) as well as 116 SNPs located in 58 genes differentially expressed in swelling buds (SB). This assay also included 20 unrelated control SNPs located in 15 housekeeping genes. Genotyping was performed on a MassARRAY System (Agena Bioscience, USA) using the iPLEX Gold chemistry following [START_REF] Gabriel | SNP Genotyping Using the Sequenom MassARRAY iPLEX Platform[END_REF]. Data analysis was performed with Typer Analyzer 4.0.26.75 (Agena Bioscience). We filtered out all monomorphic SNPs, as well as loci with a weak or ambiguous signal (i.e., displaying more than three clusters of genotypes or unclear cluster delimitation).

Raw variant data were filtered with Plink (v.1.9; [START_REF] Chang | Second-generation PLINK: Rising to the challenge of larger and richer datasets[END_REF]. Individuals and SNPs with >15% missing data were filtered out, together with variants with MAF<1%.

Linkage disequilibrium

Pairwise linkage disequilibrium (LD) estimates were obtained within each of the genetic clusters identified by genetic clustering analyses (see below) using the LD() function in the R package genetics (Warnes, Gorjanc, Leisch, & Man, 2021). The p.adjust() function was used to correct p-values for multiple testing with the Bonferroni method, and association between allele frequencies was deemed as statistically significant at a nominal significance threshold equal to 1✕10 -3 . An ad hoc algorithm was devised to iteratively identify and remove the loci most frequently involved in pairwise significant tests, leading to cluster-specific lists of SNPs showing statistical evidence of linkage. The markers shared across all analyses were removed to obtain an LD-pruned version of the SNP dataset to be used in the analyses assuming linkage equilibrium among loci (i.e., STRUCTURE and pcadapt, see below).

Population genetic structure

The Bayesian clustering analysis implemented in STRUCTURE [START_REF] Pritchard | Inference of Population Structure Using Multilocus Genotype Data[END_REF] and the multivariate method implemented in the Discriminant Analysis of Principal Components (DAPC; [START_REF] Jombart | Discriminant analysis of principal components: A new method for the analysis of genetically structured populations[END_REF] were used to infer patterns of population structuring and admixture among beech populations. The main difference between the two methods is that STRUCTURE builds genetic clusters so to minimise the overall departures from HWE, whereas DAPC is based on maximizing the differentiation between inferred genetic clusters while minimizing variation within them. STRUCTURE 2.3.4 [START_REF] Pritchard | Inference of Population Structure Using Multilocus Genotype Data[END_REF] was run using default settings and parameter values, assuming the admixture model, and the putative number of different genetic clusters (K) ranging from one to 10. Each run consisted of 5×10 4 burn-in iterations and 1×10 5 data collection iterations. Ten independent runs were performed for each value of K. The average likelihood and ΔK statistics described in [START_REF] Evanno | Detecting the number of clusters of individuals using the software structure: A simulation study[END_REF] were calculated for each K and used to identify the most-likely K-value. For informative values of K, distinct runs were averaged using CLUMPAK [START_REF] Kopelman | Clumpak: A program for identifying clustering modes and packaging population structure inferences across K[END_REF] to obtain the final estimates of the membership coefficients (q-values) at individual and population levels. To comply with model assumptions, STRUCTURE was run first with the complete SNPs dataset, and then with the LD-pruned version of the dataset (once LD estimated within each cluster).

DAPC is based on a discriminant analysis (DA) of genetic data preceded by a few analytical steps to meet its requirements, all implemented in the R package adegenet [START_REF] Jombart | adegenet: A R package for the multivariate analysis of genetic markers[END_REF]. Since DA requires a priori definition of clusters, K-means clustering of principal components (PC) on individual allele frequencies was first used to identify both group priors and the most likely number of genetic clusters. K-means was run on 150 PCs with K ranging from one to 40, and the Bayesian Information Criterion (BIC) was used to assess the best supported K-value. Then, as DA requires the variables to be uncorrelated and fewer than the number of observations, we used a principal component analysis (PCA) and the randomization approach implemented in the a.score() function in adegenet to select the number of PCs optimizing the trade-off between power of discrimination and over-fitting. Finally, the DA was run on 23 PCs extracted from the original dataset.

Basic diversity and differentiation statistics

We computed allelic richness (Ar) and mean number of alleles per locus (Na) using the R package diveRsity [START_REF] Keenan | DiveRsity: An R package for the estimation and exploration of population genetics parameters and their associated errors[END_REF]; percentage of polymorphic loci using GenAlEx 6.5 [START_REF] Peakall | GenAlEx 6.5: Genetic analysis in Excel. Population genetic software for teaching and research-an update[END_REF]; observed (HO) and expected (HE) heterozygosity, Wright's inbreeding coefficient (FIS), and βWT (a plot-specific index of genetic differentiation relative to the entire pool; Weir and Goudet 2017) using the R package hierfstat [START_REF] Goudet | Hierfstat, a package for r to compute and test hierarchical F-statistics[END_REF]. Parameters of genetic fixation (GST; [START_REF] Nei | F-statistics and analysis of gene diversity in subdivided populations[END_REF] and differentiation (Jost's D;[START_REF] Jost | GST and its relatives do not measure differentiation[END_REF] among populations/clusters were also calculated with GenAlEx and their statistical significance assessed with 999 permutations.

Isolation by distance and barriers to gene flow

We estimated spatial genetic structure (SGS) among populations and tested whether geographic distance significantly shaped the patterns of genetic differentiation (as estimated by FST/1-FST) using the software SpaGeDi 1.4c [START_REF] Hardy | Spagedi : a Versatile Computer Program To Analyse Spatial[END_REF]. To test for isolation by distance (IBD), the (FST/1-FST) values were regressed on ln(dij), where dij is the spatial distance between populations i and j. Then, we tested the regression slope (null hypothesis: blogFST = 0) using 5,000 permutations of genotypes over populations. These analyses were run both on the 64 populations and within each cluster identified with DAPC. For within-cluster analyses, we retained all the individuals successfully assigned to a given cluster (i.e., those having a q-value above the nominal threshold of 0.6). Spatial variation in genetic diversity and gene flow rates were estimated using Estimated Effective Migration Surfaces (EEMS; [START_REF] Petkova | Visualizing spatial population structure with estimated effective migration surfaces[END_REF]. This method tests for regional departures from the IBD model: areas where the decay of genetic differences across geographical distance is higher than expected under an IBD model are considered as suggestive of barriers to gene flow. A user-selected number of demes determines the geographic grid size and possible migration paths between all populations, and the EEMS are calculated by adjusting the migration rates so that the genetic differences obtained under a stepping-stone model match as closely as possible the observed genetic differences. The estimates are subsequently interpolated over the geographic space to provide a surface of observed genetic dissimilarities. We ran the runeems_snps executable with 500,000 burn-in MCMC steps and 2✕10 6 subsequent iterations. To reduce the potential influence of grid size, we averaged the results over nine independent runs with different numbers of demes (nDemes=800, 1200 and 1600 with three repetitions each) and combined the results across the three independent analyses. A higher value for nDemes produces a denser grid, which increases the flexibility of the EEMS model but also the computational cost; so we chose nDemes values for which estimations remained tractable. We assessed convergence of runs, plotted geographic distance and genetic dissimilarity across demes, and generated effective diversity (q) and effective migration rates (m) surfaces using the R package reemsplots [START_REF] Petkova | Visualizing spatial population structure with estimated effective migration surfaces[END_REF].

Bioclimatic data

Each sampling site was characterized by a set of 19 bioclimatic variables extracted from the WorldClim database v.1.4 [START_REF] Hijmans | Very high resolution interpolated climate surfaces for global land areas[END_REF] with a grid cell resolution of 30-arc second (ca 1×1 km) using DIVA-GIS v.7.5. These bioclimatic variables represent annual trends (e.g., mean annual temperature, annual precipitation), seasonality (e.g., annual range in temperature and precipitation) and extreme or limiting climatic factors (e.g., temperature of the coldest and warmest month, and precipitation of the wettest and driest quarter). To calculate predictors in the following analyses, we used the mean values of these 19 bioclimatic variables (Table S1) over the period from 1950 to 2000.

To reduce the multidimensional bioclimatic data set to a few uncorrelated factors, we performed two PCAs using the R package FactoMineR [START_REF] Lê | FactoMineR : An R Package for Multivariate Analysis[END_REF], one focusing on the temperature-related variables (BIO1 to 11), and the other focusing on the precipitation-related predictors (BIO12 to 19). The selected principal components of each PCA were used as individual climate variables in lfmm and Samβada analyses, and combined in matrices to represent the climatic structure in the variance partitioning analyses (see sections below).

Detection of signatures of selection

To detect loci carrying putative signatures of divergent selection (i.e., outliers), we used two differentiation outlier search methods, pcadapt [START_REF] Luu | pcadapt: An R package to perform genome scans for selection based on principal component analysis[END_REF] and the FSTbased method by [START_REF] Martins | Identifying outlier loci in admixed and in continuous populations using ancestral population differentiation statistics[END_REF] as implemented in lea (Frichot & François, 2015). Moreover, we used two genotype-environment associations analyses, lfmm [START_REF] Frichot | Testing for Associations between Loci and Environmental Gradients Using Latent Factor Mixed Models[END_REF] and Samβada [START_REF] Stucki | High performance computation of landscape genomic models including local indicators of spatial association[END_REF]; these approaches are detailed in Supplementary Online Appendix 2. Note that each method allows the correction of outlier detection for the confounding effects of population structure. For all four methods, p-values were corrected across multiple tests using the same local False Discovery Rate (FDR) algorithm.

Finally, we annotated the identified outliers. For the loci obtained from Lalagüe et al. ( 2014), we queried against The Arabidopsis Information Resource (TAIR 11) database using BlastX with an E-value cut-off of 10 -5 . For the loci obtained from [START_REF] Lesur | A unigene set for European beech ( Fagus sylvatica L.) and its use to decipher the molecular mechanisms involved in dormancy regulation[END_REF], the functional annotation of gene sequences containing outlier SNP was also reported (from their Table S2) .

Isolation by distance and environment

To evaluate the respective importance of IBD versus isolation by environment (IBE), we used a variance partitioning approach [START_REF] Legendre | Should the Mantel test be used in spatial analysis?[END_REF]. We partitioned the explanatory power (as expressed by the adjusted R²) of the climatic and spatial structures on the genetic structure. Genetic structure was obtained through a Principal Coordinate Analysis (PCA) on the matrix of pairwise population genetic distances as returned by GenAlEx [START_REF] Peakall | GenAlEx 6.5: Genetic analysis in Excel. Population genetic software for teaching and research-an update[END_REF]. Principal coordinates explaining up to 80% of the total variance were included in the response data table. The spatial structure was modelled by distance-based Moran's eigenvector maps (dbMEM; [START_REF] Dray | Spatial modelling: A comprehensive framework for principal coordinate analysis of neighbour matrices (PCNM)[END_REF], as suggested by [START_REF] Legendre | Should the Mantel test be used in spatial analysis?[END_REF], and estimated by the mem() function in the R package adespatial [START_REF] Dray | adespatial: Multivariate Multiscale Spatial Analysis (Version R package version 0.1-1[END_REF]. We retained only the statistically significant eigenvectors modelling positive spatial autocorrelation and, therefore, describing global patterns sensu [START_REF] Jombart | adegenet: A R package for the multivariate analysis of genetic markers[END_REF]. The climatic structure was summarized by temperature-or precipitation-related synthetic climatic variables (see "Bioclimatic data" section above). We assessed the relative contribution of climatic and spatial structure in explaining the genetic structure of populations using the function varpart() of the R package vegan (Oksanen et al., 2019). Significance of the variance components was calculated through an ANOVA-like permutation test for redundancy analysis (RDA) and partial ReDundancy Analysis (pRDA) based on 10,000 permutations [START_REF] Legendre | Numerical ecology[END_REF].

All statistical analyses were conducted using R 3.6.2 (R Core Team, 2019) unless otherwise indicated. Raw data and R scripts for analyses are publicly available [START_REF] Oddou-Muratorio | Raw data and scripts for the manuscript : " Genetic signatures of divergent selection in European beech (Fagus sylvatica L.) are associated with the variation in temperature and precipitation across its distribution range[END_REF].

Results

SNP dataset

Of the 405 SNPs from the two multiplex arrays, 135 were discarded from the analyses: 9 failed to amplify in all samples; 50 were monomorphic; 75 were of poor quality (based on visual inspection of genotyping outputs); and one had a call rate <85%. Sixteen individuals were discarded due to a low call rate. The resulting dataset comprised 430 individuals and 270 SNPs (Supplementary Online Appendix A1), with 3-10 individuals per population (6.7 on average). The 270 SNPs included 150 SNPs in 93 phenology genes, 109 SNPs in 38 stressrelated genes, and 11 SNPs in 8 housekeeping genes. After removing linked markers, the LDpruned version of the dataset comprised 212 SNPs.

Population genetic structure

Both the STRUCTURE and the DAPC analysis revealed an optimal grouping at K=3 (Fig. 2 and Fig. S1) with clear geographic boundaries among the three genetic clusters (as represented by the green, red and blue colors in Fig. 2) and admixture zones (areas with cross symbols in Fig. 2). Populations were considered as admixed when none of the population qvalues exceeded the threshold of 0.60 for any inferred genetic cluster (Table S2).

The green (dominant in Northern Europe) and red (dominant in Western Europe) clusters were separated by a main genetic boundary extending from the northern Tyrrhenian coast to southern England, which was also an area of admixture. In southern Europe, populations from the Apennines clustered with those from the Balkan peninsula and Carpathian Mountains, forming the blue cluster. In central-eastern Europe, the main boundary between the blue and green clusters was located between Slovenia and Croatia and between the Western and Eastern Carpathian Mountains. The area between the southern Carpathian Mountains and the Baltic sea represented the second largest admixture area detected. The GST and Jost's D pair-wise values among the three genetic clusters spanned from 0.024 to 0.025, and from 0.022 to 0.023, respectively.

Spatial patterns of genetic diversity and differentiation

Genetic diversity and differentiation estimates are detailed in Table S2. The percentage of polymorphic SNPs within a population (%polloc) ranged from 62% to 87 % (mean %polloc = 77%), corresponding to an allelic richness ranging from 1.48 to 1.65 (mean Ar = 1.58). Observed and expected heterozygosities per population ranged, respectively, from 0.262 to 0.347 and from 0.248 to 0.327 (mean Ho = 0.305 and mean HE = 0.300), indicating a small heterozygote excess (mean FIS = -0.025). All genetic diversity indices showed the same patterns of variation across the range (Fig. 3 a,b), where south-eastern populations had below-average values, whereas western populations showed above-average values. The indices of genetic diversity (He, %polloc) also varied among clusters (Fig. S2). Populations assigned to red or green clusters had a higher HE compared to populations assigned to the blue cluster (p = 0.003 and 0.074 respectively), while populations assigned to the red cluster showed a higher %polloc value than populations assigned to the green or blue clusters (p < 10 -3 ). Spatial patterns of genetic diversity were consistent when assessed independently with SNPs from the two different arrays (Fig. S3). Although other array-specific representation problems may occur, such a finding rules out a major distortion due to ascertainment bias. Genetic differentiation of each population from the entire gene pool ranged from -0.034 to 0.20 (mean βWT = 0.05). Patterns of βWT variation across Europe were opposed to diversity patterns (Fig. 3c), with south-eastern populations characterized by higher average values than western populations. Accordingly, βWT of the populations assigned to the red or green clusters was lower than βWT of the populations assigned to the blue cluster (p = 0.005 and 0.091 respectively) (Fig. S2c).

Genetic differentiation between all population pairs revealed a significant signal of IBD (Fig. 4a). Pairwise FST/1-FST values increased with increasing geographic distance (blog = 0.027, p < 0.001). This IBD signal was also observed within the blue and green clusters (Table S3, Fig. S4), but not within the red cluster.

The EEMS analyses highlighted several barriers to gene flow corresponding to biogeographical barriers (Fig. 4b). The first barrier separated the UK and Scandinavia from the European mainland, and corresponded to the English Channel, the Baltic sea and plains in Western Germany and North-Western France. A second barrier corresponded to the Alps (Northern Italy, Austria) and extended to Slovakia and the Carpathians in the east. Weaker barriers to gene flow also occurred in Southern Italy and the Balkans.

Spatial patterns of effective diversity estimated with EEMS (Fig. 4c) partially contrasted with the maps of HE and %polloc (Fig. 3). Indeed, areas of higher-than-average diversity were found with EEMS in Western Europe (UK, Pyrenees mountains, Germany) but also in Central and South-Eastern Europe. Diversity was lower than average in Spain, Southern Italy, and in an area from Eastern Scandinavia to Poland.

Climate data analysis

The first three principal components of the temperature-focussed PCA (Temp1, Temp2, Temp3) were retained, and accounted for 90.9% of the total variance of the dataset (Online Appendix A3). Temp1 is an axis of mean temperatures, opposing hot (southern) to cold (northern) climates. Temp2 can be interpreted as an axis of climate continentality, opposing climates with strong versus weak variation of temperatures among years and seasons (e.g. continental vs. oceanic climates). Temp3 can be interpreted as an axis of climate xericity, opposing climates with a high diurnal range and the wettest season corresponding to the coldest months (i.e., Mediterranean climates) to climates with a low diurnal range and the wettest season corresponding to the warmest months (i.e. temperate mesic climates).

The first three principal components of the precipitation-focussed PCA (Precip1, Precip2, Precip3) were also retained, and accounted for 98.6% of the total variance of the dataset. Precip1 is a precipitation abundance axis, opposing wet (Great-Britain, Northern Italy) to dry climates (Greece, Spain). Precip2 is a precipitation variability axis, opposing climates with strong (Greece, Italy) versus weak (France) variation of precipitation. Precip3 captures the coupling between precipitation and seasonal temperatures, opposing climate where high precipitation occurs during the vegetation period (Poland, Romania) to those where high precipitation occurs in winter (Greece, Italy). populations in the blue areas are connected by higher migration rates than expected under isolation by distance (IBD) while the ones in the orange areas have lower migration rates than expected and are interpreted as migration barriers. In white areas, the effective migration surface is close to the one expected under IBD. (c) Contour maps representing the posterior mean of effective diversity (q) surface; populations in the orange (respectively blue) areas have lower-than-expected (respectively higher-than-expected) genetic diversity than the average. On maps (b) and (c), black dots represent the studied populations, aggregated per grid cell (with size proportional to the number of genotyped individuals).

Selection signatures pcadapt:

The first two PCs were retained to represent population structure in pcadapt analysis based on the Cattle's rule (Online Appendix A3). One candidate SNP under selection (0.3%) was identified after controlling for FDR (Table S4). lea: The lowest cross-entropy criterion value was found at K = 3. Five SNPs (1.85%) were identified as potentially under divergent selection after controlling for FDR and after calibrating p-values using the calculated genomic inflation factor (λ = 6.0; Table S4). lfmm: After controlling for FDR, 44 SNPs (16.3%) were found to be associated with temperature or precipitation-related climatic variables (60 significant associations in total, Table S5): seven SNPs showed correlations with Temp1, 13 with Temp2 and 11 SNPs with Temp3; four SNPs showed correlations with Precip1, nine with Precip2, and 16 with Precip3.

Samβada: Thirty-two genotypes at 22 SNPs (8.1%) were associated with temperature and precipitation-related variables after controlling for FDR. In particular, four loci showed correlations with Precip1 and 11 with Precip3, while 5 SNPs were associated withTemp1, seven with Temp2 and five with Temp3 (Table S6). Overlapping signatures of selection and ontology of genes bearing outlier: Overlapping signatures of selection from population divergence and GEA analyses were detected at genes 154_1 and QB_c10512 (Table 1, Fig. 5). The GEA analyses shared signatures of local adaptation in 12 additional genes (13 common SNPs). The population divergence analyses also shared signatures of local adaptation at one additional gene. Finally, 28 and six outlier SNPs were detected by lfmm and Samβada alone, respectively. While our panel of 270 SNPs included 67.9% of putatively phenology-related, 26.4% of stressrelated and 5.7% of control-related genes, respectively, outliers included 61.4% of phenologyrelated, 34.1% of stress-related, and 4.5% of control genes (Table 1). Hence, there was no difference of category (stress, phenology, control) among initial and outlier genes (χ² = 0.99, p = 0.61).

Isolation by distance and environment

The variance partitioning and partial RDA analyses revealed a greater effect of spatial structure than of climatic structure on the spatial distribution of genetic variation (Figure 6, Table S7). Considering the 218 putatively neutral SNPs only, the climatic structure alone explained ~1% of variance in the genetic structure, which is not statistically different from the null expectation of 0% variance explained (F6,52 = 1.17, p = 0.16). On the contrary, the contribution of the spatial structure alone was much larger (R² = 16%) and its effect was statistically significant (F5,52 = 4.12, p < 0.001). The joint effect of climatic and spatial structures contributed significantly in explaining the genetic structure (R² = 25%). In this case, such a joint effect is not equivalent to a standard interaction term, and relates to the intrinsic covariation of climatic and spatial effects.

When considering the set of 52 SNPs putatively under selection, the contribution of the climatic structure to the genetic structure was significant (R² = 3%; F6,52 = 1.65, p = 0.006), but still smaller than that of the spatial structure (14%; F5,52 = 4.008, p < 0.001). The joint effects of climatic and spatial structures explained 32% of the genetic structure. Thus, the variance contributed by climatic and spatio-climatic structures combined was higher for the 52 outliers (~35%) than for the supposedly neutral loci (~26%; Fig. 6). Moreover, the temperature and precipitation components of the sole climatic effects explained a similar and statistically significant amount of variance (1.5% and 1.2%, respectively; p < 0.001; Table S7B) in the supposedly adaptive genetic structure. The variance contributed by shared spatial and temperature structures (21%) was higher than that contributed by shared spatial and precipitation structures (6%). These Venn diagrams summarize the portions of variance (adjusted R²) in the genetic structure explained by spatial and climatic effect for both sets of 218 putatively neutral SNPs (Neutral) and 52 outliers (Outliers). In left column, the effects of Temperature and Precipitation were considered together ("Climate" effect) while they were separated in right column. The intersections illustrate the contribution of joint effects, but note that this intersection is not equivalent to an interaction in an ANOVA for instance, and relates to the portion explained by effects which cannot be disentangled from each other.

Discussion

Our study supported our first expectation: we observed weak founder effects in beech, which indicates that past population demography is not likely to blur the detection of selection signatures. On the contrary, our second expectation was not met as we identified loci with the signature of divergent selection as often in genes involved in phenology as in genes involved in stress response. And counter to our third expectation, temperature and precipitation related variables were equally represented in the significant genotype-climate associations. Overall, our results suggest a balanced contribution of traits related with phenology and with stress responses to local adaptation in beech.

Impact of past recolonization history on genetic diversity

Our results revealed a clear spatial disjunction between three main gene pools. A first pool (blue cluster in Fig. 2) corresponds to the area harbouring beech glacial refugia in Southeastern Europe. Its geographical distribution matches well with results of genetic and paleoecological studies showing that lineages from these refugia expanded only as far as the Northern Apennines and Central Carpathians [START_REF] Leonardi | Genetic variability of Fagus sylvatica L. in Italy: The role of postglacial recolonization[END_REF][START_REF] Magri | Patterns of post-glacial spread and the extent of glacial refugia of European beech (Fagus sylvatica)[END_REF][START_REF] Magri | A new scenario for the quaternary history of European beech populations: Palaeobotanical evidence and genetic consequences[END_REF]. A second gene pool (red cluster) includes the Southwest European glacial refugia located on the Iberian Peninsula and in southern France, where beech persisted throughout several glacial cycles [START_REF] De Lafontaine | Cryptic no more: Soil macrofossils uncover Pleistocene forest microrefugia within a periglacial desert[END_REF] and was even more abundant than in the south-eastern refugia during the middle and upper Pleistocene [START_REF] Magri | A new scenario for the quaternary history of European beech populations: Palaeobotanical evidence and genetic consequences[END_REF]. A third gene pool (green cluster) mostly corresponds to recently recolonized areas in Northern Europe [START_REF] Sjölund | Understanding the legacy of widespread population translocations on the post-glacial genetic structure of the European beech, Fagus sylvatica L[END_REF] and likely originates from glacial refugia located in the Eastern Alps-Slovenia and in Slovakia-Moravia [START_REF] Magri | Patterns of post-glacial spread and the extent of glacial refugia of European beech (Fagus sylvatica)[END_REF]. Boundaries between the three clusters were associated with strong admixture, which would be consistent with the relatively high number of recolonization routes known for beech as compared to other tree species [START_REF] De Lafontaine | Cryptic no more: Soil macrofossils uncover Pleistocene forest microrefugia within a periglacial desert[END_REF][START_REF] Magri | Patterns of post-glacial spread and the extent of glacial refugia of European beech (Fagus sylvatica)[END_REF].

Consistently with previous studies based on allozymes or microsatellites [START_REF] Comps | Diverging Trends Between Heterozygosity and Allelic Richness During Postglacial Colonization in the European Beech[END_REF][START_REF] De Lafontaine | Stronger spatial genetic structure in recolonized areas than in refugia in the European beech[END_REF], spatial patterns of SNP genetic diversity did not reflect signals of founder effects resulting from the post-glacial expansion. In particular, the northern populations (green cluster) showed values of Nei's heterozygosity (He) and genetic differentiation (βWT) similar to those of the south-western populations (red cluster), while the south-eastern populations (blue cluster) showed lower HE and higher βWT. Hence, diversity across the 270 studied SNPs was higher, and differentiation lower, both in recently recolonized areas, and in areas where beech was more abundant in the past. These patterns likely result from the combination of several processes and lifehistory traits specific to trees in general and beech in particular: first, the long juvenile phase of forest trees strongly attenuates founder effects during colonization in a diffusive dispersal model [START_REF] Austerlitz | Effects of Colonization Processes on Genetic Diversity: Differences Between Annual Plants and Tree Species[END_REF]. Moreover, long-distance pollen dispersal is frequent in beech [START_REF] Gauzere | Ecological determinants of mating system within and between three Fagus sylvatica populations along an elevational gradient[END_REF][START_REF] Piotti | Comparison of pollen gene flow among four European beech (Fagus sylvatica L.) populations characterized by different management regimes[END_REF], which is expected to increase the number of founders and the mixing of genes from distant sources, resulting in a rapid increase of genetic diversity after the initial colonization [START_REF] Fayard | Long distance dispersal and the fate of a gene from the colonization front[END_REF][START_REF] Lander | Weak founder effects but significant spatial genetic imprint of recent contraction and expansion of European beech populations[END_REF][START_REF] Paulose | The impact of long-range dispersal on gene surfing[END_REF]. Finally, beech is one of the tree species that recolonised northern Europe the latest, and the factors that limited its ability to migrate probably also contributed to its retaining a high level of diversity along the expansion front [START_REF] Roques | Allee effect promotes diversity in traveling waves of colonization[END_REF][START_REF] Saltré | Climate or migration: What limited European beech post-glacial colonization?: European beech post-glacial migration[END_REF] Range-wide spatial genetic structure (SGS) was statistically significant but weak. The strongest signal was found in the south-eastern genetic cluster. This is consistent with the theoretical work of [START_REF] Slatkin | Isolation by Distance in Equilibrium and Non-Equilibrium Populations[END_REF] on IBD, who showed that a species having restricted dispersal should exhibit SGS if enough time has elapsed after establishment. Since the southeastern European populations (blue cluster) have undergone a relatively early and shortdistance post-glacial expansion [START_REF] Magri | A new scenario for the quaternary history of European beech populations: Palaeobotanical evidence and genetic consequences[END_REF], they would have had the longest time for the establishment of SGS.

Altogether, our results hence agree on the absence of a marked signature of genetic drift and allele surfing in beech due to recolonization. It appears therefore warranted to assume that our analysis of genetic signatures of local adaptation is little burdened with such sources of uncertainty, in line with previous studies on temperate forest trees that have explicitly tested for such effects [START_REF] Eckert | Patterns of Population Structure and Environmental Associations to Aridity Across the Range of Loblolly Pine ( Pinus taeda L., Pinaceae)[END_REF][START_REF] Ruiz Daniels | Inferring selection in instances of long-range colonization: The Aleppo pine ( Pinus halepensis ) in the Mediterranean Basin[END_REF][START_REF] Temunović | Candidate gene SNP variation in floodplain populations of pedunculate oak ( Quercus robur L.) near the species' southern range margin: Weak differentiation yet distinct associations with water availability[END_REF].

Genomic signatures of local adaptation along climatic gradients

Two genes showed convergent signatures of selection using GEA and differentiation outlier analyses. GEA detected more outliers than differentiation outlier analyses, with 12 genes showing convergent signatures of divergent selection using lfmm and Samβada. The 52 outliers identified with at least one of the methods displayed significant IBE patterns (while the putatively neutral markers did not), consistently with the fact that allele frequencies covary with climatic variables at the loci under selection.

According to GEA, 50 associations were attributable to the temperature variables and 42 to the precipitation variables. Among the temperature variables, associations with climate continentality (22) were found more often than associations with mean temperature (12) or climate xericity ( 16). This finding may indicate that the risk of late frosts could represent a major constraint for the evolution of phenology-related traits in beech [START_REF] Gauzere | Microgeographic adaptation and the effect of pollen flow on the adaptive potential of a temperate tree species[END_REF][START_REF] Kreyling | Local adaptations to frost in marginal and central populations of the dominant forest tree Fagus sylvatica L. as affected by temperature and extreme drought in common garden experiments[END_REF]. Among the precipitation variables, associations with the coupling between precipitation and seasonal temperatures (25) were found more often than associations with precipitation abundance (8) or variability (9). This result could be due to genetic differentiation between locations where high precipitation occurs during the vegetation period (coupling) versus those where a precipitation deficit occurs during the vegetation period (decoupling). This would highlight the major role of low precipitation in driving patterns of local adaptation, in agreement with the known sensitivity of beech to drought [START_REF] Aranda | Variation in photosynthetic performance and hydraulic architecture across European beech (Fagus sylvatica L.) populations supports the case for local adaptation to water stress[END_REF]Cuevo-Alarcon et al., 2021), and with the major role of maximal potential evapotranspiration as a driver of genetic differentiation for growth and survival [START_REF] Gárate-Escamilla | Range-wide variation in local adaptation and phenotypic plasticity of fitness-related traits in Fagus sylvatica and their implications under climate change[END_REF]. Partial RDA analyses also indicate significant effects of temperature and precipitation on the genetic structure at all the 52 outlier loci together, even though the portion of genetic variance contributed by "pure" temperature or precipitation effects was low in both cases (1%). Considering that phenology-related candidate genes were slightly over-represented in our set of 270 SNPs, and assuming that optimal values of phenological traits are likely to vary primarily with temperature and photoperiod rather than with precipitation [START_REF] Metcalf | Life history in a model system: Opening the black box with Arabidopsis thaliana[END_REF], the balanced contribution of precipitation and temperature variables to the genetic-climate associations suggests that range-wide local adaptation in beech is driven by traits related to various climate components (see also Garate-Escamilla et al., 2019).

Although we cannot completely rule out some false positives, the large number of outliers detected using GEA analyses compared to differentiation outlier analyses is methodologically consistent. Indeed, our sampling design with a high number of populations covering beech distribution range, at the cost of the local sampling intensity, proves to be appropriate to maximize the environmental and geographical representativeness of sampling sites, control for population structure, and ultimately optimize statistical power in GEA analyses [START_REF] De Mita | Detecting selection along environmental gradients: Analysis of eight methods and their effectiveness for outbreeding and selfing populations[END_REF][START_REF] Selmoni | Sampling strategy optimization to increase statistical power in landscape genomics: A simulation-based approach[END_REF]. As long as the total number of populations and individuals are large enough (i.e. >40 and >400 respectively, [START_REF] Selmoni | Sampling strategy optimization to increase statistical power in landscape genomics: A simulation-based approach[END_REF] and sampled across steep environment gradients, sampling a low number of individual per population is not expected to limit the power of outlier tests and GEA, especially under range expansion scenarios [START_REF] Lotterhos | The relative power of genome scans to detect local adaptation depends on sampling design and statistical method[END_REF]. Moreover, the major post-glacial expansion axes of beech align with steep ecological gradients: the South-to-North axis of expansion opposes hot to cold climates, while the axis from Central Europe to Great Britain opposes continental to oceanic climates [START_REF] Magri | A new scenario for the quaternary history of European beech populations: Palaeobotanical evidence and genetic consequences[END_REF]. Such a configuration is expected to minimize the number of false positive with GEA, as compared to the opposite case where ecological gradients are orthogonal to the expansion axis (Frichot et al., 2015).

Functional role of the genes under selection

Among the 139 candidate genes investigated, eighteen showed signatures of divergent selection with at least two methods. Eight of them had also shown signatures of divergent selection in previous studies (Table 1 andTable S8). Among these "best candidates", the outlier gene 92 encodes for ACO, an enzyme involved in the production of the plant hormone ethylene, which regulates many plant developmental processes and stress responses. This study found a putative signature of selection at the non-synonymous locus at position 352, (coding for histidine or glutamine; Table 1), where the frequency of the homozygous genotype TT decreased with drought stress (Fig. 7a). This locus was also associated with annual and growth season temperatures in the study of [START_REF] Cuervo-Alarcon | Genetic variation and signatures of natural selection in populations of European beech (Fagus sylvatica L.) along precipitation gradients[END_REF]. Two other loci within this gene (although non-coding or synonymous) were also detected as outliers by [START_REF] Pluess | Genomeenvironment association study suggests local adaptation to climate at the regional scale in Fagus sylvatica[END_REF], where their frequencies correlated with drought indices. These two previous studies were conducted in Switzerland, along drought and precipitation gradients using GEA approaches. Hence, our combined results suggest that the ACO gene could be under divergent selection at various hierarchical scales across Europe.

Another interesting example is the outlier gene QB_c10512, which encodes for the NAC domain-containing protein 72, a transcription factor responsive to desiccation. This gene was found to be differentially expressed in beech quiescent buds by [START_REF] Lesur | A unigene set for European beech ( Fagus sylvatica L.) and its use to decipher the molecular mechanisms involved in dormancy regulation[END_REF]. At synonymous position 206 (in a Leucine-coding codon), the frequency of the AA genotype increased with drought stress (Fig. 7b). Two other variants within this gene (including a nonsynonymous one) were also detected as outliers by [START_REF] Cuervo-Alarcon | Genetic variation and signatures of natural selection in populations of European beech (Fagus sylvatica L.) along precipitation gradients[END_REF], where their frequencies also correlated with precipitation during the growing season. Moreover, the gene QB13549, detected as an outlier by pcadapt and lea, also encodes for the NAC domaincontaining protein 72. At position 857, allele frequency showed a strong variation from Eastern to Western Europe (Fig. 7c).

The signatures of divergent selection were not particularly enriched in genes related to phenology. This is a counter-intuitive result, as previous quantitative genetic approaches failed to detect divergent selection at various physiological traits related to drought stress, but did detect significant differentiation in phenological [START_REF] Gauzere | Microgeographic adaptation and the effect of pollen flow on the adaptive potential of a temperate tree species[END_REF][START_REF] Hajek | Intraspecific variation in wood anatomical, hydraulic, and foliar traits in ten European beech provenances differing in growth yield[END_REF], growth and survival traits [START_REF] Gárate-Escamilla | Range-wide variation in local adaptation and phenotypic plasticity of fitness-related traits in Fagus sylvatica and their implications under climate change[END_REF][START_REF] Gauzere | Microgeographic adaptation and the effect of pollen flow on the adaptive potential of a temperate tree species[END_REF]. This is likely because the stress-related genes genotyped in this study are involved in the response to multiple stresses varying across climate gradients. Moreover, our panel of candidate genes probably determines a larger number of stress-related traits than usually phenotyped in quantitative genetic approaches [START_REF] Gauzere | Microgeographic adaptation and the effect of pollen flow on the adaptive potential of a temperate tree species[END_REF][START_REF] Hajek | Intraspecific variation in wood anatomical, hydraulic, and foliar traits in ten European beech provenances differing in growth yield[END_REF]. Another limitation of quantitative genetic approaches is that they are usually conducted on a few populations and do hence not adequately cover the range-wide diversity of stress gradients (see also Garate-Escamilla et al., 2019). This comparison illustrates the complementarity of quantitative genetic and molecular approaches to investigate local adaptation [START_REF] Rudman | What genomic data can reveal about eco-evolutionary dynamics[END_REF].

Implications for conservation and management

The rates of expected species range shifts are likely to be insufficient for trees to track ongoing climate change [START_REF] Saltré | Climate or migration: What limited European beech post-glacial colonization?: European beech post-glacial migration[END_REF][START_REF] Savolainen | Gene Flow and Local Adaptation in Trees[END_REF]. In this context, there is an increasing interest in evolutionary-oriented management strategies, relying on the high genetic diversity observed within and among tree populations to adapt forest to ongoing climate change [START_REF] Aitken | Assisted Gene Flow to Facilitate Local Adaptation to Climate Change[END_REF][START_REF] Lefèvre | Considering evolutionary processes in adaptive forestry[END_REF][START_REF] Oney | Intraspecific variation buffers projected climate change impacts on Pinus contorta[END_REF].

We showed that the spatial distribution of genetic diversity across beech distribution range reflects both biogeographical history and adaptive processes. This has consequences for conservation, where the importance of maintaining adaptive genetic diversity -in addition to preserving as many lineages as possible -cannot longer be overlooked [START_REF] De Lafontaine | Invoking adaptation to decipher the genetic legacy of past climate change[END_REF][START_REF] Ouborg | Conservation genetics in transition to conservation genomics[END_REF][START_REF] Shafer | Genomics and the challenging translation into conservation practice[END_REF]. While the conservation of lineages relies on the assessment of genetic boundaries, the conservation of adaptive diversity may require the identification of relevant loci and the targeted conservation of specific alleles, genotypes or combinations thereof. Our results also have consequences for management, where knowledge of the genetic variants under selection, as combined with the estimation of their current spatial distribution and the prediction of future climate, may help inform decisions about assisted migration, perhaps under the form of an "enrichment" of existing stands with potentially favourable genotypes [START_REF] Rellstab | Signatures of local adaptation in candidate genes of oaks ( Quercus spp.) with respect to present and future climatic conditions[END_REF][START_REF] Rochat | Spatial areas of genotype probability: Predicting the spatial distribution of adaptive genetic variants under future climatic conditions[END_REF]. This would however require a reliable validation of the adaptive meaning of our best candidates by independent proof, as well as the assessment of genotype✕genotype interactions (to make sure there is no outbreeding depression; [START_REF] Grummer | The genetics of assisted gene flow: immediate costs and long-term benefits[END_REF] and of genotype✕environment interactions (to avoid undesired, unforeseen under-performances of the introduced genotypes and their progeny in the new environments; [START_REF] Cooper | Genotypic variation in phenological plasticity: Reciprocal common gardens reveal adaptive responses to warmer springs but not to fall frost[END_REF].

At the intersection of management and conservation lies the possibility to favour natural migration and regeneration dynamics, which could result in efficient mixing of genotypes in multiple environments, thus exposing them to natural selection and adaptive processes [START_REF] Lefèvre | Considering evolutionary processes in adaptive forestry[END_REF]. Here, the detailed analysis of barriers to gene flow is of the essence, to understand whether the barriers and "corridors" we detected have been caused by geographical features or isolation by adaptation: while in the former case it may be sensible to manage such barriers and corridors to shape gene flow, in the latter in may be difficult -or even detrimental -to force the modification of gene flow patterns.

Table 1. Outlier SNPs showing signature of divergent selection with at least one of the four methods (pcadapt, lea, lfmm and Samβada). We applied analysis-specific FDR cut-offs granting no expected false positives. For lfmm and Samβada, we report the climatic variable for which the genetic-environment association was found (T1-3: Temp1-3; P1-3: Precip 1-3). Converging selection signatures from at least two methods are underlined in grey. For each SNP, we code the study where it was described first (1: [START_REF] Lalagüe | Nucleotide diversity and linkage disequilibrium at 58 stress response and phenology candidate genes in a European beech (Fagus sylvatica L.) population from southeastern France[END_REF]2: Lesur, et al., 2015) and give the gene sequence (in Genbank for 1; in [START_REF] Lesur | A unigene set for European beech ( Fagus sylvatica L.) and its use to decipher the molecular mechanisms involved in dormancy regulation[END_REF] for 2). We finally provide for each candidate gene its category (stress-related, phenology-related or control genes), its annotation based on the homology with Arabidopsis thaliana sequence (with the TAIR ID and probability of matching E)

Figure 1 :

 1 Figure 1: Distribution of the 64 studied populations (red dots) (a) in the geographical space, overlaid on beech distribution range (in grey) ; (b) in the bioclimatic niche defined by annual precipitations and temperature, with the grey colour intensity indicating increasing density of beech stands.

Figure 2 :

 2 Figure 2: Spatial interpolates of the admixture coefficients estimated with STRUCTURE for K=3. Each color corresponds to one cluster and the colour intensity indicates the probability to belong to the cluster at a given position in space, based on spatial kriging of the individual qmatrix. Only areas belonging to beech distribution range are considered. Crosses indicate admixed populations, not assigned to a single cluster.

Figure 3 :

 3 Figure 3: Estimates of diversity (HE), percentage of polymorphic loci (%polloc) and genetic differentiation relative to the entire pool (βWT) in the 64 studied populations, overlaid on beech distribution range (in grey). Each circle is a population, values above average are in red or orange, and values below the average are in blue. Circle sizes are proportional to the deviation from the mean.

Figure 4 :

 4 Figure 4: Patterns of isolation by distance and barriers to gene flow among the 64 studied populations. (a) Spatial genetic structure as depicted by the variation of genetic differentiation against geographic distance (on a log-scale). The grey envelope represents expected FST/1-FST values under complete spatial randomness and bars represent standard error at 95% level within each distance class. (b) Contour maps representing the posterior mean of effective migration (m) surface;populations in the blue areas are connected by higher migration rates than expected under isolation by distance (IBD) while the ones in the orange areas have lower migration rates than expected and are interpreted as migration barriers. In white areas, the effective migration surface is close to the one expected under IBD. (c) Contour maps representing the posterior mean of effective diversity (q) surface; populations in the orange (respectively blue) areas have lower-than-expected (respectively higher-than-expected) genetic diversity than the average. On maps (b) and (c), black dots represent the studied populations, aggregated per grid cell (with size proportional to the number of genotyped individuals).Selection signatures pcadapt:The first two PCs were retained to represent population structure in pcadapt analysis based on the Cattle's rule (Online Appendix A3). One candidate SNP under selection (0.3%) was identified after controlling for FDR (TableS4). lea: The lowest cross-entropy criterion value was found at K = 3. Five SNPs (1.85%) were identified as potentially under divergent selection after controlling for FDR and after calibrating p-values using the calculated genomic inflation factor (λ = 6.0; TableS4). lfmm: After controlling for FDR, 44 SNPs (16.3%) were found to be associated with temperature or precipitation-related climatic variables (60 significant associations in total, TableS5): seven SNPs showed correlations with Temp1, 13 with Temp2 and 11 SNPs with Temp3; four SNPs showed correlations with Precip1, nine with Precip2, and 16 with Precip3.

Figure 5 :

 5 Venn diagram of the private and common outliers identified by the different methods to detect signatures of divergent selection.

Figure 6 :

 6 Figure 6: Respective contributions of climatic and spatial structure to the genetic structure.These Venn diagrams summarize the portions of variance (adjusted R²) in the genetic structure explained by spatial and climatic effect for both sets of 218 putatively neutral SNPs (Neutral) and 52 outliers (Outliers). In left column, the effects of Temperature and Precipitation were considered together ("Climate" effect) while they were separated in right column. The intersections illustrate the contribution of joint effects, but note that this intersection is not equivalent to an interaction in an ANOVA for instance, and relates to the portion explained by effects which cannot be disentangled from each other.

Figure 7

 7 Figure 7: Variation in allelic/genotypic frequencies at three outlier SNPs (panel a: 92_352; b: QB_c10512; c: QB13549_857) across the studied geographic range. For panels (a) and (b), the graph on the left represents the predicted variation in genotype occurrence probability (GOP) across the environmental gradient as estimated by Samβada. Maps show the observed genotypic/allelic frequency superimposed on the environmental gradient (panels a and b) and spatial genetic structure (c).

Acknowledgements

We thank Dalibor Ballian, Peter Zhelev and Aristotelis C. Papageorgiou for providing beech samples from Bosnia and Herzegovina, Bulgaria and Greece, as well as Ilaria Spanu, Mariaceleste Labriola, Catia Boggi and Adline Delcamp for technical assistance during lab work. SOM and AH were supported by the French National Research Agency Biodiversa project TipTree (ANR-12-EBID-0003); AP, FB and GGV by the Italian Ministry of University and Research (FOE-2019) under the projects "Climate Change" (CNR DTA.AD003.474) and "Green & Circular Economy" (CNR DTA.AD003.139); and AH by the French National Research Agency Cluster of Excellence project COTE (ANR-10-LABX-45, CLIMBEECH). The MassARRAY genotyping was performed at the Genome Transcriptome Facility of Bordeaux (grants from EU FP7 through Trees4Future research infrastructure n° 284181, from the Conseil Régional d'Aquitaine n°20030304002FA and 20040305003FA, from the European Union FEDER n°2003227 and from Investissements d'Avenir ANR-10-EQPX-16-01).

Data Accessibility

Raw SNPs genotypes, detailed information on SNPs and sampling locations, as well as R scripts for LD analyses, for differentiation outlier analyses (with pcadapt and lea) and for GEA analyses (with lfmm and Samβada) are available at Portail Data INRAE: https://doi.org/10.15454/OMNGXR.

Author Contributions

Conceptualization: S.O.M., G.G.V., D.P., E.V., A.P. Population sampling: D.P., F.B., G.G.V, A.H and F.P. Selection of candidate genes, SNPs and associated bioinformatics tools: D.P., I.L., G.L.P, S.O.M, G.G.V. SNP genotyping: D.P. and E.G. Statistical analyses: D.P., S.O.M, E.V., A.P. Writing-Original Draft : D.P., S.O.M, E.V., A.P. All the authors reviewed, edited and approved the final manuscript. Funding acquisition: S.O.M. and G.G.V.

Supporting Information

This manuscript includes four files as supplementary materials:

• MEC-Supplemental-Information-Postolache_etal.pdf contains additional Tables andFigures •