Tiange Lang

Pierre Abadie

Valérie Léger

Thibaut Decourcelle

Jean- Marc Frigerio

Christian Burban

Catherine Bodénès

Erwan Guichoux

Grégoire Le Provost

Cécile Robin

Naoki Tani

Patrick Léger

Camille Lepoittevin

Veronica A El Mujtar

François Hubert

Josquin Tibbits

Jorge Paiva

Alain Franc

Frédéric Raspail

Stéphanie Mariette

Marie-Pierre Reviron

Christophe Plomion

Antoine Kremer

Marie-Laure Desprez-Loustau

Pauline Garnier-Géré
email: pauline.garnier-gere@inra.fr

High-quality SNPs from genic regions highlight introgression patterns among European white oaks (Quercus petraea and Q. robur)

Keywords: SNPs, functional candidate genes, Quercus robur, Q. petraea, Sanger amplicon resequencing, introgression, species differentiation

In the post-genomics era, non-model species like most Fagaceae still lack operational diversity resources for population genomics studies. We completed the analyses of Sanger sequences produced from over 800 gene fragments covering ~530 kb across the genic partition of European oaks in a range-wide sampling of 25 individuals (11 Quercus petraea, 13 Q. robur, one Q. ilex as an outgroup). Regions targeted represented broad functional categories potentially involved in species ecological preferences, and a random set of genes.

Using a high-quality dedicated pipeline, we provide a detailed characterization of over 14500 polymorphisms, including ~12500 SNPs -218 being triallelic-, over 1500 insertion-deletions, and ~200 novel di-and tri-nucleotide SSR loci. This catalog also includes various summary statistics within and among species, gene ontology information, and standard formats to assist loci choice for genotyping projects. The distribution of nucleotide diversity and differentiation across genic regions are also described for the first time in those species (mean  close to ~0.0049 in Q. petraea and to ~0.0045 in Q. robur across random regions, and mean F ST ~0.13 across SNPs). Robust patterns were observed which emphasize a slightly but significantly higher diversity in Q. petraea, across a random gene set and in the abiotic stress functional category, and a heterogeneous landscape of both diversity and differentiation. These patterns are discussed in the context of both species documented introgression history despite strong reproductive barriers. The quality, representativity in terms of species genomic diversity, and usefulness of the resources provided are discussed for possible applications in medium scale landscape ecology projects, and as a reference resource for validation purposes in larger-scale re-sequencing projects. These are preferentially recommended in oaks in contrast to SNP array development, given the large nucleotide variation and low levels of linkage disequilibrium revealed.

Introduction

High-throughput (HT) techniques of the next-generation sequencing (NGS) era and increased genome sequencing efforts in the last decade have greatly improved access to genomic resources in non-model forest tree species, but these have only been applied recently to largescale ecological and population genomics research (Neale andKremer 2011, Neale et al. 2013;[START_REF] Plomion | Decoding the oak genome: public release of sequence data, assembly, annotation and publication strategies[END_REF][START_REF] Holliday | Advances in ecological genomics in forest trees and applications to genetic resources conservation and breeding[END_REF]. One notable exception are studies undertaken in the model genus Populus (e.g. [START_REF] Zhou | Exome resequencing reveals signatures of demographic and adaptive processes across the genome and range of black cottonwood (Populus trichocarpa)[END_REF][START_REF] Geraldes | Landscape genomics of Populus trichocarpa: the role of hybridization, limited gene flow, and natural selection in shaping patterns of population structure[END_REF], Christe et al. 2016b) that benefited from the first genome sequence completed in 2006 in P. trichocarpa [START_REF] Tuskan | The genome of black cottonwood, Populus trichocarpa (Torr. & Gray)[END_REF]. In Fagaceae, previous comparative mapping and "omics" technologies (reviewed in [START_REF] Kremer | Genomics of Fagaceae[END_REF] with recent development of genomic resources (e.g. [START_REF] Faivre-Rampant | Analysis of BAC end sequences in oak, a keystone forest tree species, providing insight into the composition of its genome[END_REF][START_REF] Tarkka | OakContigDF159.1, a reference library for studying differential gene expression in Quercus robur during controlled biotic interactions: use for quantitative transcriptomic profiling of oak roots in ectomycorrhizal symbiosis[END_REF][START_REF] Lesur | The oak gene expression atlas: insights into Fagaceae genome evolution and the discovery of genes regulated during bud dormancy release[END_REF][START_REF] Lepoittevin | Single-nucleotide polymorphism discovery and validation in high density SNP array for genetic analysis in European white oaks[END_REF][START_REF] Bodénès | High-density linkage mapping and distribution of segregation distorsion regions in the oak genome[END_REF] set the path to very recent release of genome sequences to the research community (Quercus lobata, [START_REF] Sork | First draft assembly and annotation of the genome of a California endemic oak Quercus lobata Née (Fagaceae). G3[END_REF]Q. robur, Plomion et al. 2016Q. robur, Plomion et al. , 2018;;Q. suber, Ramos et al. 2018;Fagus sylvatica, Mishra et al. 2018), and these provide great prospects for future evolutionary genomics studies [START_REF] Petit | Fagaceae trees as models to integrate ecology, evolution and genomics[END_REF]Parent et al. 2015;[START_REF] Cannon | Gaining a global perspective on Fagaceae genomic diversification and adaptation[END_REF][START_REF] Lesur | Development of target sequence capture and estimation of genomic relatedness in a mixed oak stand[END_REF].

Recently, building from the European oaks genomic resources (Quercus Portal, htpps://arachne.pierroton.inra.fr, and references therein), natural populations of 4 Quercus species (Q. robur, Q. petraea, Q. pyrenaica, Q. pubescens) were genotyped for ~4000 singlenucleotide polymorphisms (SNPs, from an initial 8K infinium array, [START_REF] Lepoittevin | Single-nucleotide polymorphism discovery and validation in high density SNP array for genetic analysis in European white oaks[END_REF].

The data were further analysed [START_REF] Leroy | Extensive recent secondary contacts between four European white oak species[END_REF], with results extending previous knowledge on their likely diversification during glacial periods, as well as their recolonization history across Europe and recent secondary contacts (SC) after the last glacial maximum [START_REF] Hewitt | The genetic legacy of the Quaternary ice ages[END_REF][START_REF] Petit | Chloroplast DNA variation in European white oaks: phylogeography and patterns of diversity based on data from over 2600 populations[END_REF][START_REF] Brewer | The spread of deciduous Quercus throughout Europe since the last glacial period[END_REF]. Using recent model-based inference allowing for heterogeneity of migration rates [START_REF] Roux | Can we continue to neglect genomic variation in introgression rates when inferring the history of speciation? A case study in a Mytilus hybrid zone[END_REF][START_REF] Tine | European sea bass genome and its variation provide insights into adaptation to euryhalinity and speciation[END_REF], [START_REF] Leroy | Extensive recent secondary contacts between four European white oak species[END_REF] showed that the most strongly supported scenarios for all species pairs included very recent SC, due to a much better fit of the patterns of large heterogeneity of differentiation observed across SNP loci. These recent SC events have been documented in many patchily distributed hybrid zones where current in situ hybridization can occur among European oak species (e.g. [START_REF] Curtu | Evidence for hybridization and introgression within a species-rich oak (Quercus spp.) community[END_REF][START_REF] Jensen | Hybridization between Quercus robur and Q. Petraea in a mixed oak stand in Denmark[END_REF][START_REF] Lepais | Reproductive patterns shape introgression dynamics and species succession within the European white oak species complex[END_REF][START_REF] Guichoux | Outlier loci highlight the direction of introgression in oaks[END_REF]. The resulting low levels of differentiation among Q. robur and Q. petraea in particular is traditionally linked to a model of contrasted colonization dynamics, where the second-insuccession species (Q. petraea) is colonizing populations already occupied by the earlier pioneering Q. robur [START_REF] Petit | Hybridization as a mechanism of invasion in oaks[END_REF]. This model predicts asymmetric introgression towards Q. petraea (see [START_REF] Currat | The hidden side of invasions: massive introgression by local genes[END_REF], as often observed in interspecific gene exchanges [START_REF] Abbott | Plant introductions, hybridization and gene flow[END_REF], and a greater diversity in Q. petraea was documented at SNP loci showing higher differentiation [START_REF] Guichoux | Outlier loci highlight the direction of introgression in oaks[END_REF]. The directionality of introgression in oaks was also shown to depend on species relative abundance during mating periods [START_REF] Lepais | Species relative abundance and direction of introgression in oaks[END_REF][START_REF] Lepais | Reproductive patterns shape introgression dynamics and species succession within the European white oak species complex[END_REF].

Nethertheless, oaks like other hybridizing taxa are known for the integration of their species parental gene pools and strong reproductive isolation barriers [START_REF] Muir | Species status of hybridizing oaks[END_REF][START_REF] Muir | Evidence for shared ancestral polymorphism rather than recurrent gene flow at microsatellite loci differentiating two hybridizing oaks (Quercus spp.)[END_REF][START_REF] Abadie | Strength, diversity and plasticity of postmating reproductive barriers between two hybridizing oak species (Quercus robur L. and Quercus petraea (Matt) Liebl.)[END_REF][START_REF] Lepais | Strength and variability of postmating reproductive isolating barriers between four European white oak species[END_REF][START_REF] Ortiz-Barrientos | Species integrity in trees[END_REF]Christe et al. 2016a), raising essential questions about the interacting roles of divergent (or other types of) selection, gene flow, and recombination rates variation in natural populations, and their imprints on genomic molecular patterns of variation (e.g. [START_REF] Zhang | Recombination rate variation, hitchhiking, and demographic history shape deleterious load in poplar[END_REF]Christe et al. 2016b;[START_REF] Payseur | A genomic perspective on hybridization and speciation[END_REF]. These issues will be better addressed with genome-wide sequence data in many samples [START_REF] Buerkle | The n=1 constraint in population genomics[END_REF]), which will be facilitated in oaks by integrating the newly available genome sequence of Quercus robur to chosen HT resequencing methods (Jones and Good 2016;e.g. Zhou and Holliday 2012;[START_REF] Lesur | Development of target sequence capture and estimation of genomic relatedness in a mixed oak stand[END_REF] for the first target sequence capture in oaks).

However, obtaining high quality haplotype-based data required for nucleotide diversity estimation and more powerful population genetics inferences will likely require the development of complex bioinformatics pipelines dedicated to high heterozygosity genomes and solid validation methods for polymorphism detection (e.g. [START_REF] Geraldes | SNP discovery in black cottonwood (Populus trichocarpa) by population transcriptome resequencing[END_REF]Christe et al. 2016b). Indeed, very stringent filters are often applied in practice to limit error rates and avoid false-positives, hence limiting the impact of variable read depth and possible ascertainment bias risks, which altogether significantly decrease the number of informative loci compared to either initial fixed amounts (in genotyping arrays, e.g. [START_REF] Lepoittevin | Single-nucleotide polymorphism discovery and validation in high density SNP array for genetic analysis in European white oaks[END_REF] or potential amounts (in reference genomes, e.g. [START_REF] Pina-Martins | New insights on adaptation and population structure of cork oak using genotyping by sequencing[END_REF] in Quercus species; see also [START_REF] Van Dijk | Ten years of next-generation sequencing technology[END_REF].

Therefore, the objectives of this work were first to provide a detailed characterization of sequence variation in Quercus petraea and Quercus robur. To that end, we validated previous unpublished Sanger sequence data for fragments of targeted gene regions in a panel of individuals sampled across a large part of both species geographic range. Both functional and expressional candidate genes potentially involved in species ecological preferences, phenology and host-pathogen interactions were targeted, as well as a reference set of fragments randomly chosen across the last oak unigene [START_REF] Lesur | The oak gene expression atlas: insights into Fagaceae genome evolution and the discovery of genes regulated during bud dormancy release[END_REF]. These data were obtained within the framework of the EVOLTREE network activities (http://www.evoltree.eu/). Second, we aimed at estimating the distributions of differentiation and nucleotide diversity across these targeted gene regions for the first time in those species, and further test the robustness of comparative diversity patterns observed in the context of both species contrasted dynamics and introgression asymmetry. We discuss the quality, representativity and usefulness of the resources provided for medium scale genotyping landscape ecology projects or as a reference resource for validation purposes in larger-scale resequencing projects.

Material and methods

Sample collection

The discovery panel (DiP) included 25 individuals from 11 widespread forest stands with 2 to 4 individuals per location (13 from Q. robur, 11 from Q. petraea, 1 from Q. ilex to serve as outgroup, in Table 1). These stands occur across a large part of both Quercus species natural distributions, spanning ~20° in longitude (~2200 km) and ~11° in latitude (~1250 km) in western and central Europe (Fig. S1, Supporting Information). Individuals were chosen either on the basis of their differing leaf morphology among Q. robur and Q. petraea species [START_REF] Kremer | Leaf morphological differentiation between Quercus robur and Quercus petraea is stable across western European mixed oak stands[END_REF], or as parents of mapping pedigrees (e.g.

Choice of genic regions

Genic regions were chosen from over 103 000 Sanger sequences available in expressed sequence tags (EST) databases at the start of the project. These sequences correspond to 14 cDNA libraries obtained from various tissues and developmental stages (bud, leaf, root and wood-forming tissues), and thus likely to target a large range of expressed genes. Overall, 146 individuals were sampled in 3 different French regions (South-West, North-East and North-West). We performed the first working assembly for those sequences, with the main aim of avoiding paralog assembly while limiting split contigs with overlapping homolog sequences, the final assembly including 13477 contigs and 74 singletons (Appendices S1 and S2, Fig. S2-A, Supporting information). The libraries used in this assembly have since been named A, B, F to O, and S, and were included in larger transcriptome resources for Quercus species [START_REF] Ueno | Bioinformatic analysis of ESTs collected by Sanger and pyrosequencing methods for a keystone forest tree species: oak[END_REF].

In parallel, expressional and functional candidate genes information was compiled for targeting those potentially involved in white oaks" local and/or divergent adaptive traits (Fig. S1, Supporting information). Briefly, model species databases were searched for gene accessions by gene ontology (GO) and metabolic pathways keywords. Those sequences were first Blasted against our oak assembly [START_REF] Altschul | Basic local alignment search tool[END_REF][START_REF] Altschul | Gapped BLAST and PSI-BLAST: a new generation of protein database search programs[END_REF]. Second, the sequences from their best hits were extracted (see filtering criteria in Fig. S2-B, Supporting information) and re-Blasted against the non-redundant protein (NR) database at NCBI. Third, their annotation was compared to those of the initial gene accessions, allowing 95% of hits from the oak assembly to be validated (step 2 in Fig. S2-B, Supporting information). Expressional candidate genes sequences from bud tissues or stress treatment libraries and a random set of ESTs were also directly sampled across the oak assembly generated above (see Table S1, column F, Supporting information). Primers were designed with the OSP software [START_REF] Hillier | OSP: a computer program for choosing PCR and DNA sequencing primers[END_REF] by setting up homogenous melting temperatures constraints and excluding low-complexity propositions. Predicted amplicons were Blasted against each other and onto our assembly to exclude those with potential amplification problems and multiband patterns. They were also checked for their depth and presence of polymorphisms in contigs alignment, yielding finally 2000 amplicons for resequencing (Fig. S2-B, Supporting information).

S2-B and Table

Data production and polymorphism discovery

All the sequencing work was performed by Beckman Coulter (Agencourt Bioscience Corporation, Beverly, MA, USA) on ABI3730 capillary sequencers (Applied Biosciences)

after preparing DNA samples according to the company"s guidelines. Data quality steps were designed throughout the process in order to maximize the amount and quality of the sequences finally obtained (Fig. 1-A). Forward and reverse sequences were produced for 981 amplicons across 25 individuals (100+881 in steps 2 and 3, Fig. 1-A), and more than 85% of them yielded at least 12 high-quality sequences (Fig. 1 Supporting information for more details). Simple sequence repeat (SSR) patterns were further detected or confirmed from consensus sequences using the mreps software (Kolpakov et al. 2003, Fig. 1-D). Various additional steps involving the treatment of insertion-deletion polymorphisms (indels) and heterozygote indels (HI) in particular, allowed missing data from polymorphic diploid sequence to be minimized (see Appendix S1, Supporting information).

Functional annotation

BlastN best hits for our working assembly (orict) original contigs and for the amplified expected fragments (orict-cut) were first retrieved using [START_REF] Lesur | The oak gene expression atlas: insights into Fagaceae genome evolution and the discovery of genes regulated during bud dormancy release[END_REF]" most recent oak assembly (ocv4, see Table S2

Characterization of diversity and genetic clustering

Using the SNP-stats script for diploid data (see above), simple statistics were computed across different types of polymorphisms (SNPs, indels, SSRs…) including minimum allele frequencies (maf) and heterozygote counts, Chi-square tests probability for Hardy-Weinberg proportions, G ST [START_REF] Nei | Molecular Evolutionary Genetics[END_REF] and G ST ' standardized measure [START_REF] Hedrick | A standardized genetic differentiation measure[END_REF]. Complex polymorphisms (involving HI and/or SSRs) were also further characterized, and data formatted or analyzed using either Arlequin 3.5 [START_REF] Excoffier | Arlequin suite ver 3.5: a new series of programs to perform population genetics analyses under Linux and Windows[END_REF], SeqQual (e.g.

for Arlequin input file with phase unknown, Fig. 1-C), or R scripts. Nucleotide diversity  [START_REF] Nei | Molecular Evolutionary Genetics[END_REF], based on the average number of pairwise differences between sequences, and its evolutionary variance according to [START_REF] Tajima | Measurement of DNA polymorphism[END_REF], were also estimated and compared among species and across candidate genes grouped by broad functional categories (see column F in Table S1, Supporting information), and [START_REF] Weir | Estimating F-statistics for the analysis of population structure[END_REF] F ST estimates of differentiation were computed among species for SNP data along genic regions using analyses of molecular variance [START_REF] Excoffier | Analysis of population subdivision[END_REF].

The initial morphological species samples were compared to the genetic clusters obtained with the STRUCTURE v2.3.3 inference method [START_REF] Falush | Inference of population structure using multilocus genotype data: linked loci and correlated allele frequencies[END_REF] in order to test possible levels of introgression across individuals. We used the admixture model allowing for mixed ancestry and the correlated allele frequencies assumption for closely related populations as recommended defaults, and since they best represent previous knowledge on both species divergence across their range (e.g. [START_REF] Guichoux | Outlier loci highlight the direction of introgression in oaks[END_REF]. Preliminary replicate runs using the same sample of loci produced very low standard deviation across replicates of the data log likelihood given K (ln Pr(X/K), see Fig. S3-A, Supporting information). We thus resampled loci at random for each of 10 replicate datasets in 3 different manners to add genetic stochasticity: 1) one per region, 2) one per 100 bp block, and 3) one per 200 bp block along [START_REF] Gilbert | Recommendations for utilizing and reporting population genetic analyses: the reproducibility of genetic clustering using the program STRUCTURE[END_REF], along with R scripts for outputs). Statistical independence among loci within each species was verified with

Fisher"s exact tests implemented in Genepop 4.4 [START_REF] Rousset | Genepop'007: a complete re-implementation of the genepop software for Windows and Linux[END_REF].

Results

Polymorphisms typology and counts

Among the amplicons tested, 986 were successful, 13 did not produce any data and 23 were excluded because of paralog amplifications (Fig. 1-C and Table S1, Supporting information).

Around 25% of the successful amplicons overlapped and were merged, consistently with their original design across contigs. Despite the presence of HI patterns due to SSR or indels, most amplicons were entirely recovered with forward and reverse sequencing. Several (5% of the total) were however kept separate, either because of functional annotation inconsistency, or because amplicon overlap was prevented by the presence of SSRs or putative large introns (see "Final gene region ID" column with -F/-R suffix in Table S1, Supporting information).

We finally obtained 852 genic regions covering in total ~529 kilobases (kb), with an average size of 621 bp per region, ranging from 81 to 2009 bp (Table 2, and Appendix S4, Supporting information, for genomic consensus sequences). Compared to the EST-based expected total fragment size of ~ 357 kb, around 187 kb of intron sequence was recovered across 460 resequenced regions (assuming intron presence if an amplicon size was above its expected size by 40 bp). Introns represented ~35% of genic regions in length and ~51% of those including introns.

We observed 14102 polymorphisms in both species across 852 gene regions, 15 of those regions (<2%) being monomorphic (Table 2). This corresponds to 1 polymorphism per ~38 bp, or 1 per ~30 bp when considering the total number of variant positions in both species (17594 bp, Table 2). Remarkably, variant positions involving larger indels, SSRs and mixed complex polymorphism patterns represented ~30% of the total variant positions (Table 2, and see their exhaustive lists with various statistics in Table S3 and S4, Supporting information).

We observed 12478 SNPs (88.5% of all polymorphisms), 1 SNP per 42 bp, and 218 triallelic SNPs (~1.75% of SNPs) were confirmed by visual examination of chromatograms.

Considering only one species, 1 variant position per ~48 bp, 1 polymorphism per ~60 bp, and 1 SNP per ~68 bp were observed on average. Among indels, 1213 (8.6% of all polymorphisms) were single base, 309 ranged from 2 to 10 bp, and 102 had sizes above 10 bp which were mostly shared among species (Table 2). In this range-wide sample, there were 4334 singletons among all single base polymorphisms and 506 of them were indels. Overall, indels were present in 69% of gene regions and non-single base ones across ~30% of them.

Excluding homopolymers (see Appendix S1, Supporting information), we detected 201 SSRs occurring on 163 gene regions by considering a minimum repeat numbers of 4 and a mismatch rate among repeats below 10% (Table 2, Table S1 and Table S5, Supporting information), and 55% (111) were polymorphic in our sample of individuals (Table 2).

Among them, 89 (44%) had dinucleotide repeats and 65 (32%) trinucleotide repeats. The SSRs with the lowest number of repeats (<5) had a majority (59%) of repeat sizes between 4 and 7, the rest being trinucleotides (Table S5, Supporting information).

Given the PCR conditions for sequence acquisition and their high similarity a priori, homologous sequence data were obtained for one individual of the outgroup Quercus ilex across 37% of the gene regions (~197 kb, 397 sequences, 676 heterozygous sites in Table 2, Table S1 column Q, and see Appendix S5, Supporting information, for Q. ilex genomic sequences).

Annotations and GO term distributions

BlastX matches with E-values below 10 -30 were found for ~97% (738/764) of the contig consensus, only 11 sequences (1.4%) having hits with E-values above 10 -10 that were all among the reference random sample (see BlastX criteria in Table S2, Supporting information). The most represented species among the best hits with informative annotations were Prunus persica (111), Theobroma cacao (91), Morus notabilis (57) and Populus trichocarpa (45) (Appendix S6-A, Supporting information), which probably illustrates both the close phylogenetic relationships among Quercus and Prunus genera, consistently with results obtained on the larger ocv4 assembly [START_REF] Lesur | The oak gene expression atlas: insights into Fagaceae genome evolution and the discovery of genes regulated during bud dormancy release[END_REF], and the quality and availability of P. persica genome annotation [START_REF] Verde | The high-quality draft genome of peach (Prunus persica) identifies unique patterns of genetic diversity, domestication and genome evolution[END_REF][START_REF] Verde | The Peach v2.0 release: high-resolution linkage mapping and deep resequencing improve chromosome-scale assembly and contiguity[END_REF].

Between 1 to 30 GO terms could be assigned to 761 sequences, with EC codes and InterProScan identifiers for 343 and 733 of them respectively (Fig. 1, and Table S2, Supporting information). The most relevant GO terms were then retained using the Blast2GO "annotation rule" [START_REF] Conesa | Blast2GO: a universal tool for annotation, visualization and analysis in functional genomics research[END_REF]) that applies filters from the Direct Acyclic Graph (DAG) at different levels (Fig. S4-A-to-F, Supporting information). At biological level (BP)

3, apart from general terms involving "metabolic processes", a large number of sequences (between ~100 and ~150) were mapped to "response to..." either "...stress", "...abiotic stimulus" or "...chemical", and also to categories linked to developmental processes (Fig. S4-D, Supporting information). Enrichment tests also revealed a significant increase at both BP levels 2 and 3 for the following GO categories: "response to stress" or "external stimulus" for bud and biotic gene lists, "response to abiotic stimulus" for the bud list, and "immune" and "biotic stimulus" responses for the biotic list (see Fig. 2-B to 2-D compared to Fig. 2-A, and Fig. S5, Supporting information). Most of these exact tests (>80%) were still significant when selecting genes attributed exclusively to one particular list, which adds to the relevance of our gene lists in targeting particular functional categories.

Species assignment and introgressed individuals

In both species, the proportion of significant association tests among the loci used for clustering (> 2 million within each species) was generally one order of magnitude below the type-I error rates at 5% or 1%. This indicates a very low background LD within species at their range levels, consistently with the underlying model assumptions used in STRUCTURE.

Based on both ln Pr(X/K) and K statistics and as expected, the optimal number of genetic clusters inferred was 2, whatever the number of polymorphisms and type of sampling (Fig. 3, Fig. S3 and S6, Supporting information).

Most individuals (20) clearly belonged to either cluster with a mean probability of cluster assignment above 0.9 , which was not significantly different from 1, based on mean values of 90% Bayesian credible intervals (BCI) bounds across replicates, and for different types of sampling or SNP numbers (Fig. 3 and Fig. S6, Supporting information). Two individuals from Roudsea Wood in UK, the most northerly forest stand of this study, were significantly introgressed, each from a different cluster, with mean probabilities within 0.125 and 0.875.

These values can be considered as typical for back-crosses and later-generation hybrids [START_REF] Guichoux | Outlier loci highlight the direction of introgression in oaks[END_REF], suggesting a mixed ancestry with the other species a small number of generations in the past. In the initial morphological Q. petraea group, two individuals were clearly of recent mixed ancestry: one from the easternmost forest stand of Sopron (S444), and another one (Qs28) from central France, considered previously to be a Q. petraea parental genotype in two oak mapping pedigrees (Bodénès et al. 2012[START_REF] Bodénès | High-density linkage mapping and distribution of segregation distorsion regions in the oak genome[END_REF][START_REF] Lepoittevin | Single-nucleotide polymorphism discovery and validation in high density SNP array for genetic analysis in European white oaks[END_REF].

However, Qs28 shows here a clear F1 hybrid pattern, given its probability values close to 0.5 and its BCI maximum upper and minimum lower bound values of 0.30 and 0.61 respectively across runs (Fig. 3 and Fig. S6-A to S6-J, Supporting information). Testing 3 or 4 possible clusters showed the same ancestry patterns for the introgressed individuals with 2 main clusters and similar Q-values (data not shown), which does not support alternative hypotheses of introgression from different species in those individuals.

Large heterogeneity of diversity and differentiation across genes

Nucleotide diversity was estimated in each parental species after excluding the 4 most introgressed individuals from each initial morphological group (see above). We then checked how the remaining samples represented species" diversity. Starting with one individual, we observe a dramatic drop in the mean proportion of new variant positions brought by each new individual in any species (Mpn) as a function of the initial sample size, followed by a subsequent stabilization (Fig. 4-A, and see Fig. S7-B, Supporting information). Indeed, Mpn was only around 11% when going from 4 to 5 individuals in both species, and stabilized below 5% after 8 individuals in Q. robur (Fig. 4-A). We thus decided to retain 726 gene regions with at least 8 gametes per species (listed in column L in Table S1, Supporting information). The larger Q. robur sample after excluding the most introgressed individuals (24 versus 16 gametes in Q. petraea ) only exhibited slightly higher polymorphism counts than in Q. petraea overall (Table 3). Also, 48% and 52% of the polymorphisms observed were exclusive to Q. petraea and Q. robur respectively in our panel, the rest being shared among species (Table 3). Among exclusive polymorphisms, 46% and 44% were singletons in Q. petraea and Q. robur respectively, suggesting that they might be either rare in both species, or more polymorphic in local populations from which few individuals were sampled across the species wider ranges. Overall and within both species, we observed a large variation in number of segregating sites per gene size (Fig. S7-A, Supporting information).

The mean nucleotide diversity estimates () across genic regions when considering all polymorphisms were 0.00447 and 0.00425 in Q. petraea and Q. robur respectively, with up to a 10-fold variation among polymorphic genes overall and in different functional categories (Fig. 4-B and Table 3). When including SNPs only, mean  decreased overall by more than 10% (Table 3, and see column D in Table S4, Supporting information). The large variation among genes is also illustrated by the absence of significant differences between mean diversity among functional categories within species, in most comparisons using nonparametric Wilcoxon rank sum tests (Wrs) with similar number of genes. Two notable exceptions were observed when considering all polymorphisms: the biotic stress category (358 genes) had on average a lower  in Q. petraea than in the random gene list (211 genes, Wrs Pr<0.042), and the mean  of the reproductive phenology category was significantly lower in both species than that of the Bud phenology category (Wrs Pr<0.040 and Pr<0.013 in Q. petraea and Q. robur respectively, considering exclusive categories from Table S2, Supporting information). Genes with  estimates above 0.02 were found across most categories, whether considering all polymorphisms (Fig. 4-B) or SNPs only. The 8 genic regions showing the highest values in both species were annotated for example as disease resistance, transcription factor or membrane transport proteins, half of them being from the original random list.

Mean  comparison tests between species across all gene regions were not significant (Table 3, Wrs Pr>0.15 for all polymorphisms or SNPs only), nor were they across different categories and between gene pairs, using a 95% confidence interval based on Tajima"s evolutionary variance for  [START_REF] Tajima | Evolutionary relationship of DNA sequences in finite populations[END_REF]) while assuming underlying Gaussian distributions. Indeed for the same genic regions, many examples can be found of higher  estimates in one species or the other. However, comparing diversity estimates across the exact same positions and performing Wilcoxon paired ranked tests (Wpr) across all genes, there was a significant pattern of a slightly higher diversity in Q. petraea (see Table 3 and Fig. 4-B), whether considering all polymorphisms (Wpr Pr<0.028) or SNPs only (Wpr Pr<0.036). This pattern remained significant across the 211 genes chosen randomly (Wpr Pr<0.037, all polymorphisms), even when excluding the 5% or 10% of genes having the highest values, but it was not robust to considering the other 509 genic regions chosen in functional categories, either together or separately in the different categories (Fig. 4-B), except for the Abiotic stress category.

We also observed a very large variation for F ST estimates across gene regions and functional categories, which covered the full range of possible values [0,1], with mean values of ~0.13 whether considering all polymorphisms or SNPs only (Fig. 4-C, and Fig. 4-D for the random genic regions and a representative example of one category). The very few segregating sites with F ST values at 1 had either missing individuals" or strands, possibly caused by polymorphisms within primer regions. Among the sites sequenced for the full sample of gametes, the 20 highest F ST values ranged from 0.6 to 0.9 and belonged to 10 genic regions, many of which also showed null or very low F ST values within 100 bp. This large variation in differentiation was observed between very close variant sites in many genes, suggesting very high recombination rates at genome-wide and range-wide scales, and consistently with the observed very low LD estimates (see above). Additionally, a large variance is expected around F ST estimates due to the relatively low sample size in both species, in particular for biallelic loci [START_REF] Weir | Estimating F-statistics[END_REF][START_REF] Buerkle | The n=1 constraint in population genomics[END_REF]e.g. Eveno et al. 2008).

Discussion

In the NGS era, non-model tree species such as many Fagaceae still lag behind model species for easy access to sequence polymorphism data (but see [START_REF] Gugger | Association of transcriptome-wide sequence variation with climate gradients in valley oak (Quercus lobata)[END_REF] for Quercus lobata). These data are needed for larger scale studies addressing the many diversity issues raised by their combined economic, ecological and conservation interests (Cavender-Bares 2016; [START_REF] Fetter | Landscape Genomics of Angiosperm Trees: From historic Roots to Discovering New Branches of Adaptive Evolution[END_REF][START_REF] Holliday | Advances in ecological genomics in forest trees and applications to genetic resources conservation and breeding[END_REF]). However, recent achievements and data availability from the Q. robur genome sequence project [START_REF] Plomion | Oak genome reveals facets of long lifespan[END_REF]) opens a large range of applications in many related temperate and tropical Fagaceae species due to their conserved synteny [START_REF] Cannon | Gaining a global perspective on Fagaceae genomic diversification and adaptation[END_REF]. In this context, we discuss below the representativity of our data in terms of species genomic diversity as well as the robust patterns observed across genes, and further illustrate their past and future usefulness for Quercus species.

Genic resources content, quality, and representativity

We provide a high-quality polymorphism catalog based on Sanger resequencing data for more than 800 gene regions covering ~530 kb, using a discovery panel (DiP) from mixed Q. robur and Q. petraea populations located across a large part of their geographic range. This catalog details functional annotations, previous published information, allele types, frequencies and various summary statistics within and across species, which can assist in choosing novel polymorphic sites (SNPs, SSRs, indels…) for genotyping studies (Tables S1 to S5 [START_REF] Plomion | Decoding the oak genome: public release of sequence data, assembly, annotation and publication strategies[END_REF]. Standard formats for high-density genotyping arrays and primer information are also provided, making these resources readily operational for medium scale molecular ecology studies while avoiding the burden of bioinformatics work needed for SNP development.

Thanks to a high quality dedicated pipeline, we could perform a quasi-exhaustive characterization of polymorphism types in our DiP and across part of the genic partition of these Quercus species (see Fig. 1). Although base call error rates below 1/1000 were used (as originally developed for Sanger sequencing), most variant sites were located in regions with lower error rates (below 1/10000) so that true singletons could be identified. At the genotypic level, a Sanger genotyping error rate below 1% was previously estimated using a preliminary subset of around 1200 SNPs from this catalog (~5800 data points, [START_REF] Lepoittevin | Single-nucleotide polymorphism discovery and validation in high density SNP array for genetic analysis in European white oaks[END_REF].

This rate can be considered as an upper bound for the present study, given all additional validation and error correction steps performed. Finding an optimal balance between the number of samples and that of loci is critical when aiming to provide accurate estimates of diversity or differentiation in population genetics studies. Given the increasing availability of markers in non-model species (usually SNPs), it has been shown by simulation [START_REF] Willing | Estimates of genetic differentiation measured by F ST do not necessarily require large sample sizes when using many SNP markers[END_REF][START_REF] Hivert | Measuring genetic differentiation from Pool-seq data[END_REF]) and empirical data [START_REF] Nazareno | Minimum sample sizes for population genomics: An empirical study from an Amazonian plant species[END_REF]) that sample sizes as small as 4 to 6 individuals can be sufficient to infer differentiation when a large number of bi-allelic loci (> 1000) are being used. A broadscale geographic sampling is however required if the aim is to better infer genetic structure and complex demographic scenarios involving recolonization and range shifts due to past glacial cycles, such as those assumed for many European species [START_REF] Lascoux | The"New Wave"in plant demographic inference: more loci and more individuals[END_REF][START_REF] Keller | Genomic diversity, population structure, and migration following rapid range expansion in the Balsam Poplar, Populus balsamifera[END_REF][START_REF] Jeffries | Comparing RADseq and microsatellites to infer complex phylogeographic patterns, an empirical perspective in the Crucian carp, Carassius carassius[END_REF][START_REF] Sousa | Impact of range expansions on current human genomic diversity[END_REF]. Our sampling design is likely to have targeted a large part of both species overall diversity and differentiation across the resequenced genic regions. This is first suggested by the small proportion of additional polymorphisms once an initial sample of 8 gametes was included for each species (~10% and decreasing as sample size increases, Fig. 4-A and 7-A). Considering the DiP within each species, each individual brings on average ~166 new variants (~1% of the total). Second, the large variance observed across gene nucleotide diversity estimates (see Table 3) is mostly due to stochastic evolutionary factors rather than to sampling effects so unlikely to be impacted by sample sizes over 10 gametes [START_REF] Tajima | Evolutionary relationship of DNA sequences in finite populations[END_REF]. Also, sampling sites are located in regions which include 4 out of the 5 main cpDNA lineages along white oaks recolonization routes [START_REF] Petit | Chloroplast DNA variation in European white oaks: phylogeography and patterns of diversity based on data from over 2600 populations[END_REF], so only the less frequent D lineage from South-western Spain might not be represented in our DiP. Therefore, if new populations were being sampled within the geographical range considered, they would likely include many of the alleles here observed within species. Indeed, older and more recent reports showed a low genetic structure and high gene flow among distant populations within each species, and a much larger overall differentiation component than the local one [START_REF] Bodénès | General vs. local differentiation between two closely related white oak species[END_REF][START_REF] Mariette | Comparison of levels of genetic diversity detected with AFLP and microsatellite markers within and among mixed Q. petraea (Matt.) Liebl. and Q. robur L. stands[END_REF][START_REF] Petit | Hybridization as a mechanism of invasion in oaks[END_REF][START_REF] Muir | Evidence for shared ancestral polymorphism rather than recurrent gene flow at microsatellite loci differentiating two hybridizing oaks (Quercus spp.)[END_REF][START_REF] Derory | Contrasting relationships between the diversity of candidate genes and variation of bud burst in natural and segregating populations of European oaks[END_REF][START_REF] Guichoux | Outlier loci highlight the direction of introgression in oaks[END_REF][START_REF] Gerber | High rates of gene flow by pollen and seed in oak populations across Europe[END_REF].

We further tested the frequency spectrum representativity of our DiP by comparing genotypes for a subset of 530 SNPs (sanSNP) to those obtained using the Illumina infinium array technology (illuSNP) for ~70 individuals from Southern France local populations [START_REF] Lepoittevin | Single-nucleotide polymorphism discovery and validation in high density SNP array for genetic analysis in European white oaks[END_REF], using their Table S3). The illuSNP excluded SNPs showing compressed clusters (i.e. potential paralogs) and those showing a high number of inconsistencies with control genotypes, as recommended by the authors. For SNPs exclusive to one species in sanSNP, more than 68% of illuSNP showed either the same pattern or one where the alternative allele was at a frequency below 5% in the other species. Less than 8% of those SNPs were common in both species in illuSNP. Similarly, for singletons in sanSNP, more than two-third showed very low to low frequency (<10%) in illuSNP, and only 11% in Q. petraea and 9% in Q.

robur showed a maf above 0.25. This confirmed the singletons reality in our DiP, and also that some may represent more frequent polymorphisms in local populations. The correlations among maf in both datasets were high and significant (0.66 and 0.68 respectively for Q. petraea and Q. robur, both Pr< 0.0001). These results suggest a small risk of SNP ascertainment bias if the new resources were to be used in populations within the geographic distribution surveyed, in contrast to panels with much less individuals (see [START_REF] Lepoittevin | Single-nucleotide polymorphism discovery and validation in high density SNP array for genetic analysis in European white oaks[END_REF] for a discussion on the consequences of such bias in Quercus species).

Overall, we obtained sequence data for 0.072% (~530 kb) of the haploid genome of Q. robur (size of ~740 Mb in [START_REF] Kremer | Genome Mapping and Molecular Breeding in Plants[END_REF]). We also targeted ~3% of the 25808 gene models described in the oak genome sequencing project (www.oakgenome.fr), and around 1% of the gene space in length. Interestingly, both randomly chosen genic regions and those covering different functional categories have been mapped across all linkage groups (columns F and X

in Table S1, Supporting information). Due to the absence of observed background LD, their diversity patterns can be considered independent. The genes studied represent a large number of categories, as illustrated by very similar distributions for level 2 GO terms to those obtained with the larger ocv4 assembly (Lesur et al. 2015, comparing their Figure 2 to Fig. 4-A to 4-C, Supporting information).

Diversity magnitude and heterogeneity highlight species integrity and introgression patterns

Using a detailed polymorphism typology, we characterized for the first time in two oak species a high proportion of variant positions (30%) that included 1 bp to medium-sized indels and sequence repeats, compared to the more common and commonly reported SNP loci (Table 2). The proportions of indels observed (11.5% of all polymorphisms) is in the range of results available in model tree species (e.g. 13.8% across the genome in Prunus avium, [START_REF] Shirasawa | The genome sequence of sweet cherry (Prunus avium) for use in genomics-assisted breeding[END_REF]19% in Prunus persica, Cao et al. 2014; a minimum of 1.4% in Populus trichocarpa, Evans et al. 2014). Although less abundant than SNPs they represent an important component of nucleotide variation, often have high functional impacts when locating within coding sequences and have been proposed as an easy source of markers for natural populations studies [START_REF] Väli | Insertion-deletion polymorphisms (indels) as genetic markers in natural populations[END_REF]. Larger-sized indels are also likely to be relatively frequent in intergenic regions of the Quercus genome and have been linked to transposable elements (TE, see the BAC clones overlapping regions analyses in [START_REF] Plomion | Decoding the oak genome: public release of sequence data, assembly, annotation and publication strategies[END_REF]. Similarly, large indels and copy number variation linked to TE activity were identified as an important component of variation among hybridizing Populus species (Pinosio et al. 2017). Here when considering variant positions involved in complex polymorphisms, we observed 1 variant position per 48 bp on average within species (resp. 1 per 30 bp in both), compared to the 1 SNP per 68 bp statistic (resp. 1 SNP per 42 bp across both species). Also, some of the SNPs observed were located within complex polymorphic regions that would have been classically filtered out, and nucleotide diversity () estimates were higher by 12% when including all polymorphisms (from 0.0038 to 0.0044 if averaging across both species and all genes, Table 3). These nucleotide diversity estimates are provided for the first time in Q. petraea and Q. robur across a large number genic regions (> 800), compared to previous candidate genes studies across much smaller numbers (< 10) of gene fragments [START_REF] Kremer | Genomics of Fagaceae[END_REF] 2003;[START_REF] Arnold | Evolution through genetic exchange[END_REF], despite more recent evidence of strong reproductive barriers [START_REF] Abadie | Strength, diversity and plasticity of postmating reproductive barriers between two hybridizing oak species (Quercus robur L. and Quercus petraea (Matt) Liebl.)[END_REF]. What has been referred to as "strong species integration" seems nevertheless clearer in our Dip for Q. robur than for Q. petraea, according to genetic clustering inference without any a priori. Three individuals (27%) considered as typical morphological Q. petraea adults [START_REF] Kremer | Leaf morphological differentiation between Quercus robur and Quercus petraea is stable across western European mixed oak stands[END_REF] showed significant levels of introgression (Fig. 3). In contrast, only one Q. robur based on morphology was introgressed to a level matching the least introgressed Q. petraea individual. Discussing species delimitation, [START_REF] Guichoux | Outlier loci highlight the direction of introgression in oaks[END_REF] also showed more robustness in assigning morphological Q. robur individuals to their genetic cluster, illustrating an asymmetry in their introgression levels. We note that among our Dip individuals, Qs28, one parent from two mapping pedigrees [START_REF] Bodénès | High-density linkage mapping and distribution of segregation distorsion regions in the oak genome[END_REF]) is a clear F1 hybrid among both species (Fig. 3), making those pedigrees two back-crosses instead of one cross within species and one between species.

Moreover, after excluding the four most introgressed individuals, nucleotide diversity in Q.

petraea was significantly higher (by ~5% on average) than in Q. robur. This effect is small, detectable only with Wilcoxon paired ranked tests, mostly across the same ~200 regions sampled randomly and in the Abiotic stress category, despite the very large diversity variance across regions, and robust to excluding the highest diversity values. We also sequentially removed the three individuals with the lowest Q-values from the Q. petraea cluster (Fig. 3), since they could still harbor residual heterozygosity due to recent back-crossing events and generate the pattern observed. Remarkably, the same significant patterns of higher diversity in Q. petraea were observed. Therefore, with 8 to 10 gametes in Q. petraea instead of 8 to 24 gametes in Q. robur, and with twice less natural stands sampled, the nucleotide diversity in Q.

petraea was still slightly and significantly higher than in Q. robur (Pr<0.011 and Pr<0.026, using all polymorphisms or SNPs only respectively). Although range-wide population structure within species could differentially affect both species global diversity across our Dip, published results show that these are very small (~1% across SNPs) and similar [START_REF] Guichoux | Outlier loci highlight the direction of introgression in oaks[END_REF]. The main hypotheses developed so far to explain this difference in species diversity pattern relate to their disparities in life-history strategies for colonizing new stands and associated predictions [START_REF] Petit | Hybridization as a mechanism of invasion in oaks[END_REF][START_REF] Guichoux | Outlier loci highlight the direction of introgression in oaks[END_REF]. The colonization dynamics model and patterns observed also assumes very similar effective population sizes in both species, which is a reasonable assumption due to their shared past history and the strong introgression impact at the genomic level. However, given increasing and recent evidence of pervasive effects of different types of selection across genic regions with HT data (e.g. [START_REF] Zhang | Recombination rate variation, hitchhiking, and demographic history shape deleterious load in poplar[END_REF]Christe et al. 2016b in Populus; Chen et al 2017 for long-term perennials), alternative (and non-exclusive) hypotheses worth considering are ones of a higher genomewide impact of selective constraints in Q. robur [START_REF] Gillespie | Genetic drift in an infinite population: the pseudohitchhiking model[END_REF][START_REF] Hahn | Toward a selection theory of molecular evolution[END_REF][START_REF] Cutter | Genomic signatures of selection at linked sites: unifying the disparity among species[END_REF][START_REF] Kern | The neutral theory in light of natural selection[END_REF]e.g. Grivet et al. 2017). Since Q. robur is the most pioneering species, it has likely been submitted to very strong environmental pressures at the time of stand establishment. Selection might be efficient given oak tree reproductive capacities across a large number of genes involved in abiotic and biotic responses. This would be consistent with significantly lower levels of He in Q. robur for genes that were specifically enriched for abiotic stress GO terms [START_REF] Guichoux | Outlier loci highlight the direction of introgression in oaks[END_REF], see their Table S5). Redoing here the same tests across a larger number of independent SNPs (> 1000), whether comparing

Nei"s unbiaised locus diversity (He) or the mean He within regions among species, Q. petraea systematically showed the same trend of a slightly higher diversity overall, and significantly so only for the Abiotic stress category (Pr<0.01) and for a similar outlier SNP category (F ST >0.4, mean He>0.15, Pr<0.001) than in [START_REF] Guichoux | Outlier loci highlight the direction of introgression in oaks[END_REF]. In summary, the absence of the same pattern in any other functional categories might suggest that these are too broad in terms of corresponding biological pathways, hence mixing possible selection signals of opposite effects among species, while we still detect an overall effect due to linked selection on a random set of genes, and on genes involved in abiotic stress. More analyses comparing the nucleotide diversity patterns at genes involved in both species relevant biosynthesis pathways for ecological preferences (e.g. [START_REF] Porth | Identification of adaptation-specific differences in mRNA expression of sessile and pedunculate oak based on osmotic-stress induced genes[END_REF][START_REF] Provost | Role of waterlogging-responsive genes in shaping interspecific differentiation between two sympatric oak species[END_REF][START_REF] Cavender-Bares | Diversity, distributions, and ecosystem services of the North-American oaks[END_REF] are clearly needed in replicated populations, for example to estimate the distribution and direction of selection effects and putative fitness impact across polymorphic sites (Stoletzki and EyreWalker 2011), or to study the interplay between different types of selection and variation in local recombination rates on both diversity and differentiation patterns [START_REF] Payseur | A genomic perspective on hybridization and speciation[END_REF].

A large proportion of shared polymorphic sites (~50% in any species) highlights the close proximity of species at genomic level, consistently with a low mean differentiation across polymorphic sites (F ST ~0.13, Fig. 4-C), and despite the very large heterogeneity observed across differentiation estimates. This has now been classically interpreted (and modeled) as reflecting a strong variance in migration and introgression rates, in oaks in particular [START_REF] Leroy | Extensive recent secondary contacts between four European white oak species[END_REF], with islands of differentiation assumed to represent regions resistant to introgression. However, interpretations of such patterns remain controversial and multiple processes might be involved and worth exploring further in oaks, such as the effects of heterogeneous selection (both positive and background) at linked loci [START_REF] Cruickshank | Reanalysis suggests that genomic islands of speciation are due to reduced diversity, not reduced gene flow[END_REF][START_REF] Wolf | Making sense of genomic islands of differentiation in light of speciation[END_REF]. These effects could be particularly visible in lowrecombination regions [START_REF] Ortiz-Barrientos | Recombination rate evolution and the origin of species[END_REF], and would further interact with the mutational and recombination landscapes during the course of speciation (Ortiz-Barrientos and James 2017) and during their complex demographic history.

Applications and usefulness as reference data

During this project, several studies valued part of these resources, hence illustrating their usefulness. For example, good quality homologous sequences were also obtained for ~50 % of the gene fragments in one individual of Quercus ilex. This species is more distant to both Q. petraea and Q. robur, belonging to a different section, so these data guided the choice of nuclear genes for better inferring phylogenetic relationships across 108 oak species [START_REF] Hubert | Multiple nuclear genes stabilize the phylogenetic backbone of the genus Quercus[END_REF]. Bioinformatics tools and candidate genes annotated during the project were also useful to similar genes and SNP discovery approach in Quercus or more distant Fagaceae species [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] Lalagüe | Nucleotide diversity and linkage disequilibrium at 58 stress response and phenology candidate genes in a European beech (Fagus sylvatica) population from southeastern France[END_REF] [START_REF] Porth | Species-specific alleles at a -tubulin gene show significant association with leaf morphological variation within Quercus petraea and Q. robur populations[END_REF] for range-wide studies). Their power and scope would likely be greatly improved by using medium-scale genotyping dataset including a few thousands SNPs such as those described in our study.

The robust patterns described above of differentiation heterogeneity and consistent differences in diversity magnitude among species call for more studies at both spatial and genomic scales for unraveling these species evolutionary history, in particular regarding the timing, tempo, dynamics and genetic basis of divergence and introgression. Practically in oaks, genome complexity reduction methods such as RAD-seq and similar approaches (e.g. [START_REF] Elshire | A Robust, Simple Genotyping-by-Sequencing (GBS) Approach for High Diversity Species[END_REF]) might be fairly limiting for the research questions mentioned above (Arnold et al. 2013;[START_REF] Henning | Genetic mapping of horizontal stripes in Lake Victoria cichlid fishes: benefits and pitfalls of using RAD marker for dense linkage mapping[END_REF][START_REF] Zhou | Targeted enrichment of the black cottonwood (Populus trichocarpa) gene space using sequence capture[END_REF], especially given the large variance in nucleotide diversity and low overall differentiation characterized here. We therefore do not recommend the development of a very large SNP array in oaks since it is likely to be very costly for a minimal return, especially given the very large and range-wide panel that would be needed to significantly limit ascertainment bias (see [START_REF] Lepoittevin | Single-nucleotide polymorphism discovery and validation in high density SNP array for genetic analysis in European white oaks[END_REF]. Also the very high SNP density required for targeting of functional variants given low overall levels of LD (indicating potentially high recombination rates), would be technically constrained when controlling for genotyping error rates (shown previously to be high).

Indeed, these rates would probably be stronger for high diversity, complex, duplicate or multiple copy genic regions (observed in this study Table S1, supporting information, and shown recently to have evolutionary impact on the Q. robur genome structure, [START_REF] Plomion | Oak genome reveals facets of long lifespan[END_REF], preventing them even to be targeted. In contrast, targeted sequence capture (TSC) strategies for resequencing (Jones and Good 2016), although still uncommon in forest tree species evolutionary studies, might be most useful and efficient since they can be oriented towards recovering long genomic fragments and thus allow more powerful site frequency spectrum and haplotype-based inferences to be pursued (e.g. [START_REF] Zhou | Exome resequencing reveals signatures of demographic and adaptive processes across the genome and range of black cottonwood (Populus trichocarpa)[END_REF][START_REF] Wang | Natural selection and recombination rate variation shape nucleotide polymorphism across the genomes of three related Populus species[END_REF], at the same time avoiding most of the above technical issues. TSC approaches will surely be encouraged and tailored to specific evolutionary research questions in oaks in the next decade, given the new Q. robur genome sequence availability [START_REF] Plomion | Oak genome reveals facets of long lifespan[END_REF][START_REF] Lesur | Development of target sequence capture and estimation of genomic relatedness in a mixed oak stand[END_REF] for the first TSC in oaks). However, the bioinformatics pipelines needed for validating haplotype-based or quality data for population genetics inferences also need constant reassessment according to research questions and chosen technology.

We thus propose, in addition to direct applications to landscape genetics (detailed above) and transferability to other Quercus species (see [START_REF] Chen | A novel set of single-copy nuclear gene markers in white oak and implications for species delimitation[END_REF], and primer information in Table S1, Supporting information), that the high-quality data characterized in this study serve as a reference for such validation purposes. They could not only help for adjusting parameters of the chosen pipelines for data outputs, but also allow estimating genotyping error rates for SNP and more complex classes of variants, either from general patterns comparisons (e.g maf distribution from Tables S3, S4 Supporting information) or using the same control individuals maintained in common garden that could be included in larger-scale studies. Such a reference catalog of SNPs and other types of polymorphisms within gene fragments could also be very useful for solid cross-validation of variants identification, allele frequency and other derived summary statistics in alternative strategies such as Pool-Seq, which allow increasing genomic coverage while sampling cost-effectively by pooling individuals [START_REF] Schlötterer | Sequencing pools of individuals-mining genome-wide polymorphism data without big funding[END_REF].

Indeed, the drawback of these approaches, despite dedicated software (PoPoolation2, [START_REF] Kofler | PoPoolation2: identifying differentiation between populations using sequencing of pooled DNA samples (Pool-Seq)[END_REF]) is that they can give strongly biaised estimates, or ones that do not consider evolutionary sampling [START_REF] Hivert | Measuring genetic differentiation from Pool-seq data[END_REF] Latitude and longitude are given in the WGS 84 coordinate system. Coordinates correspond either to a central point in the mixed forest stand, or the mean of individual trees coordinates. *: parents of controlled crosses used for genetic mapping. Quercus species a priori assigned from morphological information by persons who sampled the trees.

Figures Legends

Figure 1 Bioinformatics strategy for sequence data production, amplicon assembly, functional annotation, and polymorphism discovery. Scripts used are in italics (see text for further details). GO: Gene Ontology, EC: Enzyme Commission ID. * A good sequence is defined as having a minimum of 50% of its nucleotides with a Phred score above 30.

  2014). These diversity levels are also within the range estimated for the long-term perennial outcrosser category in Chen et al. (2017, see Fig.1-D with a mean value of silent close to ~0.005) and can be considered relatively high in the plant kingdom if excluding annual outcrosser estimates or intermediate otherwise. In oaks as in many other tree species with similar life history traits, these levels would be consistent with their longevity, large variance in reproductive success and recolonization or introgression histories, which could have maintained deleterious loads of various origins[START_REF] Zhang | Recombination rate variation, hitchhiking, and demographic history shape deleterious load in poplar[END_REF][START_REF] Chen | Genetic diversity and the efficacy of purifying selection across plant and animal species[END_REF], Christe et al. 2016b) Comparing the nucleotide diversity distributions and examining the range of differentiation across genic regions in our Dip reveal several robust patterns that altogether illustrate historical introgression patterns among both Quercus species, which have long been considered as iconic examples of species exhibiting high levels of gene flow (e.g:Petit et al.

Figure 2

 2 Figure 2 Distributions of GO terms across different gene lists (bud, abiotic and biotic) at Biological Level 3, and Fisher exact tests across pairs of sequence clusters with the same GO terms between the random list and other lists. Significance levels *: P<0.05, **: P<0.01.

Figure 3

 3 Figure 3Posterior assignment probabilities of individuals into two optimal clusters based on STRUCTURE analyses. Probabilities are sorted in increasing order of belonging to cluster 2 (here Q. robur (Qr, in blue/dark grey), the alternative cluster 1 matching Q. petraea (Qp, in green/light grey), apart from individuals with higher introgression levels. Each bar represents one individual and includes mean upper and lower bounds of 90% Bayesian confidence intervals around mean Q-values across 10 replicates. Each replicate is a different random sample of 1785 polymorphisms.

Figure 4

 4 Figure 4Mean proportion of new variant sites brought by each new distinct individual added to all possible initial sample size combinations (-A); Mean nucleotide diversity (considering all polymorphisms) in both species across genic regions, and different functional categories (-B) compared between species with Wilcoxon signed-rank tests: significant at Pr<5% (*), nonsignificant (ns); Histogram of Fst estimates across polymorphic gene regions with a minimum of 8 gametes per species, after excluding singletons and grouping negative with null values (-C); Manhattan plot of Fst estimates sorted by mean Fst values across randomly chosen (black dots) and Bud phenology genic regions (grey dots) (-D).

  Fig. 1

Fig. 4

 4 Fig. 4

  [START_REF] Bodénès | High-density linkage mapping and distribution of segregation distorsion regions in the oak genome[END_REF], see Table1). Leaves were sampled, stored in silica gel and sent to INRA (Cestas, France) for DNA

extraction following

[START_REF] Guichoux | Outlier loci highlight the direction of introgression in oaks[END_REF]

. DNA quality and concentration were assessed with a Nanodrop spectrophotometer (NanoDrop Technologies, Wilmington, 152 DE, USA) and by separating samples in 1% agarose gels stained with ethidium bromide. Extractions were repeated until we obtained at least 20 micrograms of genomic DNA per sample, which was needed for a few thousands individual PCRs.

  [START_REF] Zhou | Targeted enrichment of the black cottonwood (Populus trichocarpa) gene space using sequence capture[END_REF][START_REF] Zhou | Exome resequencing reveals signatures of demographic and adaptive processes across the genome and range of black cottonwood (Populus trichocarpa)[END_REF], Evans et al. 2014[START_REF] Wang | Natural selection and recombination rate variation shape nucleotide polymorphism across the genomes of three related Populus species[END_REF]; P. tremula: ~0.008, P. tremuloides:~0.009 across genic regions,[START_REF] Wang | Natural selection and recombination rate variation shape nucleotide polymorphism across the genomes of three related Populus species[END_REF]; ~0.0026 to 0.0045 in a panel including wild Prunus persica accessions,Cao et al.

in Q. petraea; e.g.

[START_REF] Homolka | Insights into drought adaptation of two European oak species revealed by nucleotide diversity of candidate genes[END_REF]

. Our estimates are consistent with those obtained from genome-wide data and range-wide panels in angiosperm tree species, available mostly from the model genus Populus (e.g. P. trichocarpa: 1 SNP per 52 bp and ~0.003 across genic regions,

  in Fagus sylvatica, El Mujtar et al. 2014 in Nothofagus species). Given the low ascertainment bias and good conversion rate expected within the range surveyed, those genomic resources would be directly applicable to landscape genomics studies at various spatial scales (reviewed in[START_REF] Fetter | Landscape Genomics of Angiosperm Trees: From historic Roots to Discovering New Branches of Adaptive Evolution[END_REF] in both Quercus species. Indeed, easy filtering on SNP statistics provided in the catalog would allow distinguishing among different classes of SNPs (e.g. exclusive to each species, common and shared by both, linked to particular GO functional categories), or delimiting and tracing species in parentage analyses and conservation studies (e.g.[START_REF] Guichoux | Outlier loci highlight the direction of introgression in oaks[END_REF][START_REF] Blanc-Jolivet | Tracing the origin and species identity of Quercus robur and Quercus petraea in Europe: a review[END_REF], or improving estimates of lifetime reproductive success and aiming to understand how demographic history and ecological drivers of selection affect spatial patterns of diversity or isolating barriers[START_REF] Andrew | A road map for molecular ecology[END_REF] e.g. Geraldes et al. 2014). This type of spatial studies are rare in the target species, usually included a small number of SSR markers, and all suggested complexity in geographical patterns of genetic variation and importance of the ecological context (e.gNeophytou et al. 2010;[START_REF] Lagache | Distinct male reproductive strategies in two closely related oak species[END_REF][START_REF] Klein | Demographic and spatial determinants of hybridization rate[END_REF][START_REF] Beatty | Pure species in a continuum of genetic and morphological variation: sympatric oaks at the edge of their range[END_REF] for local or regional studies;[START_REF] Muir | Evidence for shared ancestral polymorphism rather than recurrent gene flow at microsatellite loci differentiating two hybridizing oaks (Quercus spp.)[END_REF][START_REF] Gerber | High rates of gene flow by pollen and seed in oak populations across Europe[END_REF]

Table 1

 1 . Therefore, they require further validation methods which usually value previously developed high-quality and lower-scale data (e.g. Pool-Seq versus Sanger and Rad-Seq in Christe et al. 2016b; Illumina GA2 versus Sanger in Cao et al. 2014; EUChip60K versus deep-whole genome resequencing in Silva-Junior et al. 2015).Finally such a reference dataset would help optimizing the amount of available data from either TSC or whole-genome resequencing experiments in future research challenges.Original assembly used for selecting contigs is in Appendix S2 (Supporting information). For Sanger trace files (with data on at least 2 individuals), see the Dryad repository (doi:10.5061/dryad.h380d51). Consensus sequences are respectively in Supporting information appendices S1 (used to design primers), S5 (used for functional annotation), and S6 (genomic sequences obtained). TablesS1 and S2correct and extend and the oak Candidate .evoltree.eu/index.php/erecources/databases/candidate-genes). SNPs, indels and SSRs catalogs and positions within genomic consensus sequences, and ready-to-use format for genotyping essays are provided in Tables S3 to S5 (Supporting information). SNP data correct and largely expand part of the SNP database of the Quercus Portal (www.evoltree.eu/index.php/snp-db/, "Quercus petraea / robur" field). Geographic location of 25 sampled individuals from Quercus petraea, Q. robur and Q. ilex.

	Data Accessibility

Acknowledgments

The authors thank Alexis Ducousso, Jean-Marc Louvet, Guy Roussel, Pablo Goicoecha, Hervé le Bouler, Félix Gugerli, Csaba Matyas, Sandor Bordacs, Hans P. Koelewijn, Joukje Buiteveld, Stephen Cavers, Bernd Degen and Jutta Buschbom for choosing trees and providing dried leaves of individuals from various Intensive Study Populations of previous European projects populations. We are grateful to H. Lalagüe, G. Vendramin, I. Scotti, and L. Brousseau for testing earlier scripts of SeqQual and to I. Lesur for help in using the ocv4 oak resources. The sequencing work was funded by the EVOLTREE network of Excellence (EU contract n°016322). TL post-doc fellowship was funded by the ANR TRANSBIODIV (06-BDIV-003-04) and LINKTREE (contract n°2008-966). TD salary was funded by the ANR REALTIME (N°59000256). Computing facilities of the Mésocentre de calcul Intensif Aquitain des Universités de Bordeaux, de Pau et des Pays de l"Adour are thanked for providing computer time for this study. We also thank Rémy Petit for funding part of TL fellowship and support in developing SeqQual tools. PA received a Ph.D. grant (2009-2011) from the « Ministère de l"Education Nationale, de l"Enseignement Supérieur et de la Recherche » of France, and additional funding from EVOLTREE.

Counts for Q. petraea exclude the 2 most introgressed individuals (Qs28 and S444 in Table 1); SSR: simple sequence repeats; "N50 size" is the size for which the cumulative sum of gene amplicons' size equal or higher than this value corresponds to 50% of the total amplicons' size sum; The number of polymorphisms for Q. ilex equals the number of heterozygotes in the resequenced individual across amplicons; Numbers of monomorphic regions were computed for those with at least 10 gametes in both species; Some detected SSR patterns were not polymorphic in our samples (detailed in Tables S1 andS5, supporting information). 2. Diversity are computed for regions with a minimum of 200 bp overall and at least 8 gametes per species at variant positions. The 509 chosen genes belong to the different functional categories listed in Table S1. Values in the "both species" column for nucleotide diversity estimates are means across each species values.

Supporting Information

Fig. S1 Sampling site locations within the natural geographic distribution of Q. petraea and Q. robur. Vector map is from http://www.naturalearthdata.com and distribution areas from Euforgen (http://www.euforgen.org/distribution-maps/) Fig. S2 Working assembly steps and softwares (A), and bioinformatic strategy for search of candidate genes and amplicon choice (B). -A, -B, -C, -D, -E), and of the mean values of the estimated probability ln (of the data given K) with standard deviations for K ranging from 1 to 5 (S3-F to S3-J), which show support for K=2. Plots are from the STRUCTURE HARVESTER program.

Fig. S4

Distributions of Gene Ontology (GO) terms for the consensus sequences in Appendix S3, at level 2 (-A, -B, -C) and level 3 (-D, -E, -F): A-and D-for Biological Process, B-and Efor Molecular Function, C-and F-for Cellular Component. Annotation rules: E-value<10 -30 , annotation cut-off 70, GO weight 5, HSP coverage cutoff 33%. Filtering applies for at least 5 sequences and a node score of 5 per GO term (but see rare exceptions in Table S2). Table S1 Description of amplicons: primer sequences, original candidate gene list, targeted biological functions (see references), candidate gene type, fragment expected size and position in original assembly, preliminary results based nucleotide quality for obtained sequences and final decision after excluding paralog amplifications.

Table S2 Functional annotation results from Blast2GO (-A), comparison of BlastX best hits results (according to E-values) between consensus sequences of the orict working assembly and the ocv4 assembly (-B), and comparisons of BlastN results of both orict and corresponding amplicon only to ocv4 (-C).

Table S3 Description of all variants single base positions, with sample sizes, alleles, genotypes counts, various statistics, and generic format for genotyping essays input data. Species samples exclude the 2 most introgressed individuals.

Table S4

Description of all polymorphisms as in Table S3, with additional sequence motifs, contiguous base positions and length for complex ones (derived from Table S3, and see also Table S5 for SSR positions).

Table S5 SSR patterns as detected from the mreps software.

Appendix S1 Additional method details.

Appendix S2 Contigs of the original working assembly used for selecting candidate gene regions and design amplicon primers, including consensus sequences and reads where nucleotides with Phred score below 20 have been masked.

Appendix S3

Sequences of chosen contigs consensus and singletons sequences for functional annotation analyses.

Appendix S4 Consensus sequences of 852 genomic regions obtained in this study for Quercus petraea and Q. Robur individuals. "(N) 9 " : represents a low-quality fragment of a length below ~1 kb separating Forward and Reverse amplicons; "n" position with a majority of nucleotides with phd score below 30. "(-) x ": insertion is a minor allele at that position.

Appendix S5 Nucleotide sequence data of 394 gene regions for one Quercus ilex individual, heterozygote sites being indicated by IUPAC codes.

Appendix S6 Outputs from Blast2GO analyses.