Charles Perrier
email: charles5perrier@icloud.com

Quentin Rougemont

Anne Charmantier

Demographic history and genomics of local adaptation in blue tit populations

Keywords: adaptation, ecological genetics, population genetics | genomics

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

| INTRODUC TI ON

Local adaptation in heterogeneous environments imposing divergent selection on counterparts of a same species has fascinated scientists for decades [START_REF] Blanquart | A practical guide to measuring local adaptation[END_REF][START_REF] Kawecki | Conceptual issues in local adaptation[END_REF]. However, it is only with the recent advance in sequencing technologies that numerous recent empirical studies started discovering the genomic features and processes underlying local adaptation in heterogeneous environments [START_REF] Hoban | Finding the genomic basis of local adaptation: Pitfalls, practical solutions, and future directions[END_REF][START_REF] Manel | Genomic resources and their influence on the detection of the signal of positive selection in genome scans[END_REF][START_REF] Savolainen | Ecological genomics of local adaptation[END_REF][START_REF] Tiffin | Advances and limits of using population genetics to understand local adaptation[END_REF]. Especially, many studies used so-called F ST genome scans between populations [START_REF] Excoffier | Detecting loci under selection in a hierarchically structured population[END_REF][START_REF] Narum | Comparison of FST outlier tests for SNP loci under selection[END_REF][START_REF] De Villemereuil | Genome scan methods against more complex models: When and how much should we trust them?[END_REF] to detect genes under selection, and potentially implicated in local adaptations [START_REF] Hohenlohe | Population genomics of parallel adaptation in threespine stickleback using sequenced RAD tags[END_REF][START_REF] Lamichhaney | A beak size locus in Darwin's finches facilitated character displacement during a drought[END_REF][START_REF] Perrier | Do genetic drift and accumulation of deleterious mutations preclude adaptation? Empirical investigation using RADseq in a northern lacustrine fish[END_REF][START_REF] Reid | The genomic landscape of rapid repeated evolutionary adaptation to toxic pollution in wild fish[END_REF]. Nonetheless, genome scans need improvements, for example, for studying polygenic variation [START_REF] Rockman | The QTN program and the alleles that matter for evolution: All that's gold does not glitter[END_REF][START_REF] Wellenreuther | Detecting polygenic evolution: problems, pitfalls, and promises[END_REF]. Genome scans also need to be interpreted in the light of recombination variation along the genome [START_REF] Booker | Variation in recombination rate affects detection of FST outliers under neutrality[END_REF][START_REF] Stapley | Variation in recombination frequency and distribution across eukaryotes: Patterns and processes[END_REF]) since it shapes the potential extent of diversity and of differentiation in the different regions of the genome [START_REF] Burri | Linked selection and recombination rate variation drive the evolution of the genomic landscape of differentiation across the speciation continuum of Ficedula flycatchers[END_REF][START_REF] Cutter | Genomic signatures of selection at linked sites: Unifying the disparity among species[END_REF][START_REF] Rougemont | Demographic history, linked selection, and recombination shape the genomic landscape of a broadly distributed Pacific salmon[END_REF][START_REF] Tine | European sea bass genome and its variation provide insights into adaptation to euryhalinity and speciation[END_REF]. Genome scans results need to be interpreted in the light of the demographic history of populations since (a) spatio-temporal variation in gene flow and of effective population size affect the possibility of local adaptations [START_REF] Lenormand | Gene flow and the limits to natural selection[END_REF] and the level of false positive due to genetic drift and (b) regions implicated in reproductive isolation accumulated during allopatry showing increased differentiation upon secondary contacts (SC) can be misinterpreted as genuine footprints of recent divergent selection [START_REF] Bierne | Pervasive selection or is it…? why are FST outliers sometimes so frequent?[END_REF][START_REF] Bierne | The coupling hypothesis: Why genome scans may fail to map local adaptation genes[END_REF]. Lastly, structural variants such as inversions also require special attention, since they have the potential to protect adaptive gene sets from recombination [START_REF] Kirkpatrick | Chromosome inversions, local adaptation and speciation[END_REF][START_REF] Kirkpatrick | How and why chromosome inversions evolve[END_REF][START_REF] Stapley | Variation in recombination frequency and distribution across eukaryotes: Patterns and processes[END_REF][START_REF] Wellenreuther | Going beyond SNPs: The role of structural genomic variants in adaptive evolution and species diversification[END_REF]) and hence enable their persistence and rapid redeployment in heterogeneous environments despite high gene flow [START_REF] Lowry | A widespread chromosomal inversion polymorphism contributes to a major life-history transition, local adaptation, and reproductive isolation (NH Barton[END_REF][START_REF] Sinclair-Waters | Ancient chromosomal rearrangement associated with local adaptation of a post-glacially colonized population of Atlantic Cod in the northwest Atlantic[END_REF][START_REF] Todesco | Massive haplotypes underlie ecotypic differentiation in sunflowers[END_REF].

Here, we aimed at investigating various demographic and genomic aspects of the adaptive divergence of a well-studied passerine bird, the Blue tit (Cyanistes caeruleus). Populations of small insectivorous passerines have long been used to study local adaptation, in both quantitative genetics and population genetics frameworks [START_REF] Broggi | Local adaptation to winter conditions in a passerine spreading north: A common-garden approach[END_REF][START_REF] Carbonell | Effects of habitat heterogeneity and local adaptation on the body condition of a forest passerine at the edge of its distributional range[END_REF][START_REF] Laaksonen | Sympatric divergence and clinal variation in multiple coloration traits of Ficedulafly catchers[END_REF]. In particular, several blue tit populations breeding in heterogeneous habitats in Southern France (Figure 1a) offer an ideal context to study local adaptation. Four of them (two deciduous and two evergreen, circled in black in Figure 1a) have been subject to a long-term project spanning more than 40 years [START_REF] Blondel | A thirty-year study of phenotypic and genetic variation of blue tits in Mediterranean habitat mosaics[END_REF][START_REF] Charmantier | Mediterranean blue tits as a case study of local adaptation[END_REF]. These populations show marked quantitative phenotypic differences (Figure 1b,c), notably in morphological (e.g., tarsus length and body mass), life-history (lay date and clutch size), and behavioral traits (e.g., song characteristics and handling aggression) [START_REF] Charmantier | Mediterranean blue tits as a case study of local adaptation[END_REF]. These phenotypic differences were found at two spatial scales. First, F I G U R E 1 Map of the sampling locations of blue tit populations on the mainland and in Corsica (a) and phenotypic differences (b). In red, habitats dominated by deciduous oaks and in green, habitats dominated by evergreen oaks. Sites with long-term monitoring are circled in black on the map (A1D corresponds to D-Rouvière in Charmantier et al. 2016, B3E to D-Pirio, B4D to D-Muro and B4E to E-Muro). On figures b & c are presented illustrative phenotypic differences (with standard deviations) between the four main populations with long-term monitoring. Are shown laying date (days), clutch size and male tarsus length (mm) (data from table 1 in [START_REF] Charmantier | Mediterranean blue tits as a case study of local adaptation[END_REF] Evol Appl, see this reference for detailed information for these measures and for differences in many other traits)

(a) (b) (c)
birds breeding in deciduous forest habitats are taller, more aggressive, and lay larger and earlier broods than birds in evergreen forests (see table 1 in [START_REF] Charmantier | Mediterranean blue tits as a case study of local adaptation[END_REF]. Strikingly, neighboring populations in deciduous and evergreen habitats are weakly genetically differentiated [START_REF] Dubuc-Messier | Gene flow does not prevent personality and morphological differentiation between two blue tit populations[END_REF][START_REF] Porlier | Habitat-linked population genetic differentiation in the blue tit Cyanistes caeruleus[END_REF][START_REF] Szulkin | Population genomic footprints of fine-scale differentiation between habitats in Mediterranean blue tits[END_REF] despite the short spatial scale, which questioned the mechanisms of persistence of the observed phenotypic differentiation against presumably large gene flow. Second, insular blue tits from Corsica, that might have been isolated since the sea level rise after the last glacial maximum (the sea level raised of 120 m from 17,000 to 5,000 years ago [START_REF] Jouet | Shoreface migrations at the shelf edge and sea-level changes around the Last Glacial Maximum (Gulf of Lions, NW Mediterranean)[END_REF][START_REF] Lambeck | Sea-level change along the French Mediterranean coast for the past 30 000 years[END_REF]), are smaller and more colorful than their mainland counterparts (again, see table 1 in [START_REF] Charmantier | Mediterranean blue tits as a case study of local adaptation[END_REF] and are listed as a separate subspecies (Cyanistes caeruleus caeruleus on mainland Europe and

Cyanistes caeruleus ogliastrae mainly in Corsica and Sardinia). Overall, traits displaying these strong phenotypic differences had heritabilities ranging from 0.20-0.43 (e.g., for lay date, [START_REF] Caro | Local adaptation of timing of reproduction: Females are in the driver's seat[END_REF]) to 0.42-0.60 (e.g., for tarsus length, [START_REF] Delahaie | Conserved G-matrices of morphological and life-history traits among continental and island blue tit populations[END_REF][START_REF] Perrier | Heritability estimates from genomewide relatedness matrices in wild populations: Application to a passerine, using a small sample size[END_REF][START_REF] Teplitsky | Assessing multivariate constraints to evolution across ten long-term avian studies[END_REF]), are classically related with fitness, and hence could be involved in a local adaptation process. In particular, the heritability of lay date and the breeding time gap between populations could result in reproductive isolation by breeding time, limiting gene flow, and favoring local adaptation. The studied traits were typically quantitative [START_REF] Charmantier | Mediterranean blue tits as a case study of local adaptation[END_REF], hence probably controlled by a polygenic architecture [START_REF] Perrier | Heritability estimates from genomewide relatedness matrices in wild populations: Application to a passerine, using a small sample size[END_REF] composed of many loci with small individual effects, as found in similar traits for other passerine birds [START_REF] Bosse | Recent natural selection causes adaptive evolution of an avian polygenic trait[END_REF][START_REF] Hansson | Contrasting results from GWAS and QTL mapping on wing length in great reed warblers[END_REF][START_REF] Lundregan | Inferences of genetic architecture of bill morphology in house sparrow using a high-density SNP array point to a polygenic basis[END_REF][START_REF] Santure | Genomic dissection of variation in clutch size and egg mass in a wild great tit (Parus major) population[END_REF]. Overall, given their phenotypic, demographic, and genetic characteristics, these blue tit populations provide an ideal framework to study the genomic architecture of polygenic adaptation in heterogeneous environments.

We investigated genome-wide patterns of genetic diversity and differentiation and the demographic history between several populations of blue tits from Southern France, in heterogeneous forest habitats (deciduous vs. evergreen) and in insular (Corsica island) and mainland areas (Mainland France) in order to better understand the determinants of their local adaptation. The analysis was based on

Resequencing of birds sampled in four sites studied in the context of a long-term project [START_REF] Blondel | A thirty-year study of phenotypic and genetic variation of blue tits in Mediterranean habitat mosaics[END_REF][START_REF] Charmantier | Mediterranean blue tits as a case study of local adaptation[END_REF] together with three additional pairs of deciduous and evergreen forests in order to test for parallel evolution (see Figure 1a). First, we investigated variation in genetic diversity and differentiation in order to verify that habitat type and geographic distance explained a significant proportion of the genetic structure between populations [START_REF] Szulkin | Population genomic footprints of fine-scale differentiation between habitats in Mediterranean blue tits[END_REF]. Second, we investigated the historical and contemporary demography of each deciduous and evergreen population pair in order to better understand the origin of their differentiation, and the demography between the populations on mainland potentially implicated in adaptation in the two habitat types and between Corsica and Mainland France. Rather than identifying the genes under selection, we wanted to test whether the adaptation had a polygenic or oligogenic architecture and whether outlier loci would be more frequent in regions with low recombination rate. We discuss our results in the light of the recent literature in population genomics, to decipher the role of genomic and demographic variations on the genetic and phenotypic divergence between blue tits in deciduous versus evergreen habitats and in mainland versus insular contexts.

| ME THODS

| Sites and sampling

Nine hundred forty seven blue tits (Cyanistes caeruleus) were captured in 10 locations (Figure 1a, Table 1). Three locations are located in the South of mainland France (codes staring with "A": A1D, A2D, A2E) and seven in Corsica (staring with "B": B3E, B4D, B4E, B5D, B5E, B6D, B6E). The numbers from 1 to 6 indicated the sampling area and hence the population pairs. Half of these sites were dominated by the deciduous downy oak Quercus pubescens (indicated by "D" in the code), and the other half are dominated by the evergreen holm oak Q. ilex (indicated by "E"). Four of these locations (A1D, B3E, B4D, B4E) are monitored as part of a long-term research program started in 1976 and have been described in previous studies [START_REF] Blondel | A thirty-year study of phenotypic and genetic variation of blue tits in Mediterranean habitat mosaics[END_REF][START_REF] Charmantier | Mediterranean blue tits as a case study of local adaptation[END_REF]. The six other sites (A2D, A2E, B5D, B5E, B6D, and B6E) were sampled in order to obtain further replicates of neighboring deciduous and evergreen populations.

Capture and handling of the birds was conducted under permits provided by the Centre de Recherches sur la Biologie des Populations d'Oiseaux (CRBPO) and by the Direction Départementale des Services Vétérinaires (DDSV). Birds were captured during the reproductive period, from early April to late June, on their nesting territory. Birds were captured either in nest boxes during the feeding of their young (in the four sites studied on the long term) or using nets in the vicinity of their nest (in the other sites). Birds were banded with a unique metal ring provided by the CRBPO. Five to 20 μl of blood was sampled from a neck or a wing vein from breeding blue tits. Blood was stored at 4°C in Queen's buffer [START_REF] Seutin | Preservation of avian blood and tissue samples for DNA analyses[END_REF].

| Molecular biology and sequencing

DNA extractions were achieved using Qiagen DNeasy Blood & Tissue kits and were randomized across sites. DNA was quantified using first a NanoDrop ND8000 spectrophotometer and then a Qubit 2.0 fluorometer with the DNA HS assay kit (Life Technologies).

DNA quality was examined on agarose gels. Library preparation using RAD-seq (restriction-site-associated DNA sequencing; [START_REF] Baird | Rapid SNP discovery and genetic mapping using sequenced RAD markers (JC Fay[END_REF]) with the enzyme SbfI was done by Montpellier GenomiX (MGX) platform (CNRS, Montpellier). Each individual was identified using a unique six nucleotides tag, individuals were multiplexed in equimolar proportions by libraries of 36 individuals, and each library was sequenced on a lane of an Illumina HiSeq 2000. Single-end sequencing was used to produce 100 bp reads. This design was used to obtain an average read depth of 50x. The DNA of three individuals was extracted twice and processed in different libraries to test for reliability of the genotyping process.

| Bioinformatics and data filtering

Raw sequences were inspected with FastQC [START_REF] Andrews | FastQC: A quality control tool for high throughput sequence data[END_REF] for quality controls. Potential fragments of Illumina adapters were trimmed with Cutadapt (Martin 2011), allowing for a 10% mismatch in the adapter sequence. Reads were filtered for overall quality, demultiplexed, and trimmed to 85bp using process_radtags, from the Stacks software pipeline 1.39 [START_REF] Catchen | Stacks: An analysis tool set for population genomics[END_REF], allowing for one mismatch in the barcode sequence.

Sequencing RAD-tags resulted in a median value of 5,449,564 reads per individual. BWA-MEM 0.7.13 (Li & Durbin, 2009) was used to map individual sequences against the reference genome of the great tit

Parus major [START_REF] Laine | Evolutionary signals of selection on cognition from the great tit genome and methylome[END_REF] and to produce sam files using default options. Indeed, although great tit and blue tit diverged about 7-14 Ma [START_REF] Päckert | Calibration of a molecular clock in tits (Paridae)-Do nucleotide substitution rates of mitochondrial genes deviate from the 2% rule?[END_REF], the use of a reference genome over a de-novo approach is very often highly recommended (see, e.g., [START_REF] Rochette | Deriving genotypes from RADseq short-read data using Stacks[END_REF]. It is also interesting to note that synteny is typically highly conserved in birds [START_REF] Backström | A gene-based genetic linkage map of the collared flycatcher (Ficedula albicollis) reveals extensive synteny and gene-order conservation during 100 million years of avian evolution[END_REF], suggesting that aligning reads on a divergent species would still allow analyses based on SNP position along the genome. On average, 93% of the raw reads were mapped against the genome (Table S1).

Samtools 0.1.19 (Li et al., 2009) was used to build and sort bam files.

We used pstacks to treat bam files, align the reads as assembled loci, and call SNPs in each locus. We used a minimum depth of coverage (m) of 5, the SNP model, and alpha = 0.05 (chi-square significance level required to call a heterozygote or homozygote). cstacks was used to build the catalog of loci, allowing three mismatches between sample loci when building the catalog. sstacks was used to match loci against the catalog. Lastly, populations program in Stacks was used to genotype individuals. In this program, relatively permissive filters were applied, to retain SNP genotyped in at least 50% of individuals (all individuals from all sites grouped) and with a minimum read depth of 4. 350,941 SNP and 947 individuals were obtained in this "GeneralDataset."

We hence applied additional filters to this GeneralDataset using the programs VCFtools [START_REF] Danecek | The variant call format and VCFtools[END_REF] and Plink [START_REF] Purcell | PLINK: A tool set for whole-genome association and population-based linkage analyses[END_REF]. We filtered for a minimum average read depth of 10 and a maximum average read depth of 100 (corresponding approximately to the 5%-95% distribution of read depth). SNP genotyped in <80% of the individuals were removed. We removed SNPs with observed heterozygosity ≥0.65 among individuals on Corsica or on mainland to reduce the potential occurrence of stacked paralogues.

Individuals below 85% genotyping rate were removed. Identity-bystate (IBS) was estimated for three replicated individuals (from the DNA extraction to SNP calling) in order to investigate reliability of the entire genotyping process. The IBS measured between replicates was high, ranging from 0.9989 to 0.9999, indicating very low genotyping error rate. These replicates were then removed from the dataset. Using the R packages gdsfmt and SNPRelate [START_REF] Zheng | A high-performance computing toolset for relatedness and principal component analysis of SNP data[END_REF],

we measured the realized genomic relatedness (GRM) between individuals in order to remove highly related individuals (i.e., full-sibs and parent-offspring). For each pair of individuals with relatedness ≥0.35, we removed one individual. This procedure was applied in order to limit biases due to highly related individuals. Indeed, since we sampled breeding birds at their nests for the four sites studied on the long term and that they tend to disperse relatively close to their nests, we expected a higher percentage of related individuals than there actually is in the population. We then limited the number of individuals to 60, chosen at random, in populations (A1D, B3E, B4D, & B4E) in which a large number of individuals were genotyped.

This limitation was achieved in order to limit differences in analyses precision between populations due to unequal sample sizes (all individuals will be used in another ongoing study of the genomic architecture of quantitative adaptive traits). We hence removed potential monomorphic SNPs and created the dataset "FilteredDataset." 454 individuals (Table 1) and 144,773 SNPs were kept in this final dataset used for all the genomic analyses excepted analyses of demographic history using ABC. The median genotyping rate across all SNPs for these individuals was 0.981. The median read depth across genotypes (SNPs × Individual) was 49.2x. The number of SNPs per chromosome ranged from 143 (LGE22) to 15,738 (chromosome 2).

Further filtering (e.g., MAF) was often operated, depending on the analyses, and thereafter mentioned if it was the case. To perform analyses of demographic history using ABC analyses, we produced a haplotype VCF file for these 454 aforementioned individuals with the populations module of Stacks and we filtered it as explained in the Note S1.

| Analysis of genetic structure and effects of environmental variables

We used PCAs with the function snpgdsPCA from snprelate to depict genetic structure between individuals. PCAs were run for the entire dataset, for mainland and Corsican populations separately and for each pairs of deciduous and evergreen populations. We inferred admixture proportions for each individual using Admixture 1.23 [START_REF] Alexander | Fast model-based estimation of ancestry in unrelated individuals[END_REF] with K-values ranging from 1 to 12 and 1,000 bootstraps. The different clustering solutions were inferred by plotting cross-validations errors and by plotting individual admixture proportion. The effect of environmental variables, forest phenology (E vs. D), geographic distance (latitude and longitude), and altitude, on genomic differentiation, was measured using a redundancy analysis (RDA) [START_REF] Forester | Detecting spatial genetic signatures of local adaptation in heterogeneous landscapes[END_REF][START_REF] Legendre | Comparison of the Mantel test and alternative approaches for detecting complex multivariate relationships in the spatial analysis of genetic data[END_REF] as implemented in the R package Vegan [START_REF] Oksanen | The Vegan package[END_REF]. We investigated the proportion of genetic variability explained by a constraining covariance matrix consisting of phenology, latitude, longitude, and altitude for each individual. We tested the global significance of the model using 1,000 permutations. We ran marginal effects permutations to address the significance of each variable. Then, we focused on the effect of phenology alone, using partial RDA to take into account the effect of latitude, longitude, and altitude. Significance was tested running 1,000 permutations. For the PCA, the admixture analysis and the RDA, we selected from the FilteredDataset the SNPs with more than 95% genotyping rate, MAF >0.05, retaining one SNP per locus and we remove SNP in linkage disequilibrium using the Plink command "indep 50 5 2." Genome-wide differentiation between each sampling location was measured with Weir and Cockerham's F ST estimator [START_REF] Weir | Estimating F-statistics for the analysis of population structure[END_REF] implemented in StAMPP [START_REF] Pembleton | StAMPP: An R package for calculation of genetic differentiation and structure of mixed-ploidy level populations[END_REF]. F ST was estimated for all the SNPs, the ones on the autosomes and the ones on the sex chromosome Z separately.

Significance was assessed using 1,000 bootstraps replicates.

| Analysis of demographic history

Alternative models of divergence history including the effects of selection at linked sites affecting Ne and of differential introgression (m) were compared using an ABC framework modified from [START_REF] Roux | Shedding light on the grey zone of speciation along a continuum of genomic divergence[END_REF]. Linked selection either under the form of hitchhiking of neutral loci linked to a selective sweep (Maynard Haigh & Smith, 1974) or under the form of background selection [START_REF] Charlesworth | The effect of deleterious mutations on neutral molecular variation[END_REF] has strong effects in regions of low recombination and have been shown to influence model choice and parameter estimates [START_REF] Ewing | The consequences of not accounting for background selection in demographic inference[END_REF][START_REF] Schrider | Effects of linked selective sweeps on demographic inference and model selection[END_REF]. The same is true when populations accumulate reproductive incompatibilities during the divergence process: The resulting barrier to gene flow reduces the effective migration rate along the genome [START_REF] Barton | The barrier to genetic exchange between hybridising populations[END_REF] and not accounting for it can affect demographic model choice and parameter estimates [START_REF] Roux | Crossing the species barrier: Genomic hotspots of introgression between two highly divergent Ciona intestinalis species[END_REF][START_REF] Sousa | Identifying loci under selection against gene flow in isolation-with-migration models[END_REF]. Moreover, including selected loci in demographic inferences can reveal the deeper origins of population divergence [START_REF] Bierne | Pervasive selection or is it…? why are FST outliers sometimes so frequent?[END_REF]. Six scenarios were compared for the four pairs of deciduous and evergreen populations in order to test whether the divergence between different habitats was not due to a divergence in different historical refugia but to a contemporary ecological divergence. We included a model of panmixia Under AM model, gene flow occurs from T split to T am and is followed by a period without gene flow. Under IM, gene flow is continuous after T split . Under, the SC model, T split is followed by a period of SI, after which a SC starts T sc generations. The EQ, IM, AM, and SC models included migration as M = 4 N 0 .m, with M 1←2 being the number of migrants from population 2 to population 1 and M 2←1 being the reverse. The effect of linked selection and barriers to gene flow were accounted for, by modeling heterogeneous population size (Ne) and heterogeneous migration (m), respectively. Such heterogeneity was modeled using beta distributions as hyper-prior on each of these two parameters. These resulted in four alternative versions for models with gene flow EQ, AM, IM, and SC (NhomoMhomo, NhomoMhetero, NheteroMhomo, and NheteroMhetero) and two versions for PAN and SI (Nhomo and Nhetero). We used a modified ABC [START_REF] Csilléry | Approximate Bayesian Computation (ABC) in practice[END_REF] pipeline from [START_REF] Rougemont | The demographic history of Atlantic Salmon (Salmo salar) across its distribution range reconstructed from approximate Bayesian computations[END_REF] to perform model selectionand to estimate posterior probabilities of parameters.

We also investigated the historical demography of the populations from Corsica as compared to the ones from the mainland. Gene flow has probably been impossible since the last deglaciation and sea level rise. Therefore, we compared models of ancient migration (AM) and SI. The pipeline described in the above section, integrating linked selection and barriers to gene flow, was run between the A1D samples (chosen on the mainland because it had the largest sample size) and B4D (chosen in Corsica because it had both a large sample size and was from the same habitat as A1D). We used the same ABC pipeline for the model selection procedure and parameter estimation as described in the above section. Finally, we attempted to convert demographic parameter into biological units assuming a mutation rate of 1e-8 mutations/bp/generations.

| Genomic diversity

Genome-wide genetic diversity was inferred for each sampling location in the dataset by measuring observed heterozygosity (Ho), proportion of polymorphic loci, and MAF spectrums. For each population and chromosome, and subsequently on average for the entire genome, linkage disequilibrium (LD) decay was measured with Plink and smoothed in R. To contrast with long-term Ne estimates from coalescent simulation in our ABC modeling (see "Analysis of demographic history"), we also inferred recent Ne for each population using SNeP V1.1 [START_REF] Barbato | SNeP: A tool to estimate trends in recent effective population size trajectories using genome-wide SNP data[END_REF], which uses LD data, with a MAF ≥0.05 filter per population. We investigated the nature of SNP variation, that is, synonymous or nonsynonymous, blasting all RAD sequences on the reference genome and the transcriptome of the great tit (ftp://ftp.ncbi.nih.gov/genom es/Parus_major/; Santure, Gratten, Mossman, Sheldon, & Slate, 2011) using blastx [START_REF] Mcginnis | BLAST: At the core of a powerful and diverse set of sequence analysis tools[END_REF]. We kept hits with at least 90% similarity and a minimum amino acid sequence length alignment of 25. We kept only SNPs for which both the alternative and reference allele yields the same score. Finally, we tested for differences in the distribution of run of homozygosity (ROH) between the mainland and Corsica that may have resulted from smaller Ne and larger inbreeding in Corsica versus the mainland. We used plink 1.9 to estimate the length and number of ROH. We required a window of 500 kb to be homozyguous in order to be considered as a ROH, and with a maximum of 100 SNP.

| Identification of genomic footprints of selection

We used three methods to search for outlier SNPs potentially under divergent selection between blue tit populations. First, we used Bayescan V2.1 [START_REF] Foll | A genome-scan method to identify selected loci appropriate for both dominant and codominant markers: A Bayesian perspective[END_REF] We filtered each of the five datasets for a minimum MAF of 0.05.

We used default parameters except for prior odds that were set at 10,000 in order to limit false positives. We investigated the paral- 2018) to list loci with outlier loading scores on the first RDA axes. We compared the loci found using these different methods.

We reported in which genes these outliers were found (the list of genes can be found together with the genome published by Laine et al., 2016 on NCBI).

| Variation of genomic differentiation with recombination rate

We investigated variation of F ST with local recombination rate and whether SNP outliers were more often found in regions of low recombination than elsewhere in the genome. We estimated local recombination rates using a coalescent method implemented in Ldhat (McVean, 2004) using linkage disequilibrium signal. Following these authors' recommendation, the dataset was split in blocks of 2,000 SNPs with 500 overlapping SNPs. A MAF of 0.05 and a maximum of 5% of missing data were allowed. The local recombination rate, ρ = 4Ne r, was then estimated in each block independently with the Bayesian reversible jump MCMC scheme implemented in interval. We used a block penalty of 5, with 30 million MCMC iterations sampled every 5,000 iterations. The first 250,000 iterations were discarded. To speed up computations, we used the precomputed two locus likelihood table for n = 190 and assuming theta = 0.001.

We estimated ρ in a composite dataset with individuals from every population, habitats, and sex. We tested for correlations between SNP F ST and recombination rate using linear models, for Corsicamainland and for deciduous-evergreen differentiation, independently, and we represented the correlation using a LOESS fit. We tested whether outliers found using Bayescan and the RDA method for both Corsica-mainland and deciduous-evergreen differentiation were more often found in regions of low recombination than elsewhere in the genome using chi-square (χ 2 ) tests, and we represented the pattern using histograms.

| Detection of genomic inversions

We searched for potential genomic inversions using a variety of descriptive statistics. First, we searched for genomic regions having a particularly low recombination rate and large long-distance linkage disequilibrium, nevertheless associated with a high density of SNPs and therefore unlikely to correspond to peri-centromeric regions but rather to local suppression of recombination that may be due to inversions. Second, we implemented a PCA sliding window analysis in order to identify portions of the genome with individuals carrying an inversion at the homozygous or heterozygous state or individuals exempt from the inversion [START_REF] Ma | Investigation of inversion polymorphisms in the human genome using principal components analysis[END_REF]. We first used 10Mb windows sliding by 1 Mb and then 1 Mb windows sliding by 100 kb, in order to use enough SNPs to perform PCAs. We then used Lostruct [START_REF] Li | Local PCA shows how the effect of population structure differs along the genome[END_REF], with k = 2, sliding by 100 SNPs, to identify particular blocks of linked SNPs explaining an abnormally high proportion of variance between two groups of individuals (e.g., inverted and noninverted).

We detected one putative inversion. We verified, using admixture with K = 2 for analyzing SNPs from the inversion, that putative heterozygous individuals had an admixture ratio close to 1:1 of both putative inverted and noninverted homozygous clusters. We then inspected variations of F ST (per SNP) and π (per 10kb window, using vcftools) along the genome between individuals that were inverted homozygous, noninverted homozygous, and heterozygous for the inversion, looking for potentially reduced diversity and increased differentiation at the inversion. We looked for potential salient variations in read depth in and around the putative inversion. We also tested whether the detected inversion was at Hardy-Weinberg equilibrium (using a common χ 2 test) and whether the frequency of the inversion varied geographically and between evergreen and deciduous habitats.

To study the history of the putative inversion identified, we aimed at measuring intra-and interspecific genetic distance and absolute divergence, at the inversion and for the entire genome, for several blue tits and great tits. To do that, we first generated a new SNP dataset by running stacks with the same pipeline and parameters as explained earlier, with two mainland great tits, two insular great tits, four mainland blue tits homozygous for the inversion, four mainland blue tits not carrying the inversion, and four Corsican blue tits not carrying the inversion (the inversion was not found in Corsica). Second, we selected the SNPs from the region of the inversion. Third, Da was measured between the five aforementioned groups of individuals using PopGenome [START_REF] Pfeifer | PopGenome: An efficient Swiss army knife for population genomic analyses in R[END_REF]. Using this measure of Da, we estimated approximately the inversion apparition time using T = Da/2µ Hudson, Kreitma, Aguadé, 1987), assuming a standard mutation rate of 1e-8 and with the simplifying assumption of no gene flow and no introgression. Lastly, we represented the divergence between these individuals using an UPGMA tree of bitwise distance using the R package poppr 1.1.1 [START_REF] Kamvar | Poppr: An R package for genetic analysis of populations with clonal, partially clonal, and/or sexual reproduction[END_REF].

(

| Gene ontology

We used the R package topGO [START_REF] Fe R E N C E S Alexa | Gene set enrichment analysis with topGO[END_REF] to investigate the potential gene ontologies (GO) that were statistically enriched for the sets of genes identified among outliers and the inversion, compared to the entire list of genes in which all the SNPs from the entire dataset were found. For the outlier tests, we used GO analyses for each gene list obtained using the different outlier identification methods but also for aggregated lists of all of the gene identified for deciduous-evergreen tests or for Corsica-mainland tests. We used the GO referenced for the zebra finch, T. Guttata (tguttata_gene_ensembl). We report a Fisher enrichment test pvalue and a p-value after applying a Benjamini and Hochberg correction for multiple testing to control for false discovery rate.

| RE SULTS

| Genetic structure and effects of environmental variables

The admixture analysis suggested the existence of two main distinct genetic groups corresponding to mainland and island populations (Figure 2a). Increasing K values contributed to delineate the populations, showing the existence of a weak structure between the six main sampling areas. However, the coefficient of variation increased with K. There was very little evidence for genetic structure within each of the four pairs of D-and E-populations. PCAs revealed a clear structure between individuals from the mainland and the island but little structure within each of these two groups (Figure S1). The full RDA model was significant (p < .001, Table S2, Figure 2b), as well as the effects of each variable tested in the model: phenology All pairwise F ST comparisons were significantly different from zero (Table 2). Average F ST calculated between island and mainland populations was 0.081. It was on average 0.003 between close populations within the mainland and 0.005 within the island. F ST was low, on average 0.003, for each of the four pairs of D-and E-populations, ranging from 0.0006 between B6D and B6E to 0.0079 between B4D and B4E. The F ST estimated on the Z chromosome was on average 1.52 times higher than on the autosomes; 1.62 times higher when comparing populations from the mainland versus Corsica , and 1.44 times higher when comparing neighboring populations in the island or in the mainland.

| Demographic history

When deciphering the historical demography of deciduous and evergreen population pairs, the hierarchical model choice procedure strongly rejected models of SI, ancient migration, and of PAN, which were associated with a posterior probability of 0 (Table S3). Instead, models with gene flow were highly supported. In three of the four pairs, the equilibrium model (EQ) received the highest posterior probability with P(EQ) = 0.99 and 0.95, 0,94 in the B6D versus B6E, B5D versus B5E and B4D versus B4E comparisons, respectively. In the A2D versus A2E comparisons, the best supported model was the isolation with migration model, P(IM) = 0.74 with the second best model being SC, with P(SC) = 0.23. Comparing the two models against each other while excluding all remaining models yield unambiguous support for IM, with P(IM) = 0.99 (Table S3). Across all models, comparisons with heterogeneous gene flow and heterogeneous effective population size were not supported, indicating that, if genetic barriers or linked selection were at play they could not be detected. Demographic parameters were estimated for each pair of populations under the best model (Table S4). Posterior distributions were well differentiated from their prior indicating that estimated parameters were confidently estimated (Figure S2). Effective population size (N e ) was slightly higher in deciduous than in evergreen habitats (Ne-ABC in Table 2, Figure 3a, Table S4) and tended to be higher on the mainland than on the island. In a number of comparisons, our migration estimates reached the prior upper bounds.

Therefore, we ran a new set of simulations with wider priors (Note S1). Migration rates were not different between deciduous and evergreen habitats: In two instances, the rate of migration was higher from the deciduous to the evergreen, and in two other instances, the reverse was true (Figure S3).

Deciphering the historical demography between Corsica and the mainland, the ancient migration models (AM) strongly outperformed the SI models, with p(AM) = 0.99 (Table S5). showing the influence of latitude, longitude, altitude, and phenology on genetic structure

A2D A2E A1D B3E B4D B4E B5D B5E B6D B6E K = 2 K = 6 K = 10 CV = 0.302 CV = 0.307 CV = 0.314 (a) Admixture (b) RDA -2 -1 0 1 2 3 -2 -1 0 1 2 PHENOLOGY LATITUDE LONGITUDE ALTITUDE A1D A2D A2E B3E B4D B4E B5D B5E B6D B6E Mainland Corsica -2 -1 0 1 2 -2 -1 0 1 2 PHENOLOGY ALTITUDE RDA1 RDA2 RDA3 RDA4 | 1153 PERRIER Et al.
Ne-ABC in Table 2, Table S6). Our analysis also S6).

| Genetic diversity

The mainland populations displayed significantly higher patterns of observed heterozygosity (on average 0.169, S5a). LD decayed rapidly in the first 5kb and was lower in populations from the mainland (especially A1D) compared to island populations (Figure S5b). This pattern of rapid LD decay was similar between chromosomes (see, e.g., chromosomes 1, 2, and Z, Figure S6). Contemporary N e inferred from LD varied from 142 (in A2E) to 355 (in A1D) (Ne-LD in Table 2, Figure S5c). N e values were rather similar within pairs of E-and D-populations, although Ne was on average smaller in evergreen populations (245 on average) than in deciduous ones (261 on average). The four largest N e (with an average of 318) were found for the four populations monitored on the long term in forest in which hundreds of artificial nest boxes have been installed for several decades (the N e in the 6 other populations was on average of 210). In the entire dataset, 2,599 SNPs were identified as nonsynonymous while 6,751 SNPs were identified as synonymous.

We Table S7). The Z chromosome harbored a proportion of nonsynonymous mutations 1.83 times higher than the average for the autosomes. We also observed significant differences in the count of ROH among populations from the mainland versus those from the island but not among populations within the mainland or the island (ANOVA p < 2e-16; Table S8 for Tukey HSD test).

| Footprints of selection

Manhattan plots of F ST per SNP revealed no clear evidence for high peaks of F ST for the five comparisons considered (First panels of (Table S9, Figure S9). None of these SNPs were detected twice among the different tests conducted with Bayescan. Among these 40 SNPs, 27 were found for the Corsica versus mainland test, 2 were found for A2D versus A2E, 9 were found for B4D versus B4E, 0 were found for B5 versus B5E, and 2 were found for B6D versus B6E. Among these 40 SNPs, 19 genes were identified (Table S9). Sliding windows of F ST S11, Figure S10). For the deciduous versus evergreen test, only one of the RDA outliers was also found outlier in Bayescan tests. For the Corsica-mainland test, eight of the RDA outliers were also found outliers in the Bayescan test. Genes in which the outliers were found are reported in the Table S13.

| Variation of genomic differentiation with recombination rate

F ST between mainland and Corsican populations was negatively correlated to recombination rate (Figure 6a, linear model p < 2e-16;

Figure S11). In contrast, F ST between each pair of deciduous and evergreen populations was not correlated to recombination rate. F ST outlier SNPs between Corsica and the mainland populations and identified by Bayescan or by the RDA were more often found in regions of low recombination (χ 2 test p-values <.01, Figure 6c) than observed for the entire SNPs (Figure 6b). Average recombination rate was on average twice lower for these F ST outlier SNPs compared to the rest of the SNPs (t test p-values < 1e-6). In contrast, F ST outlier SNPs between deciduous and evergreen populations were not preferentially found in regions of low recombination (Figure 6d).

| Genomic inversions

We detected one putative inversion on chromosome 3, spanning 7a). However, the SNP density did not appear reduced at the location of this putative inversion (Figure 7a), suggesting a recent drop of recombination (i.e., dissimilar to what is expected in a peri-centromeric location). We did not find evidence for an increase of read depth at specific loci in the putative inversion or its vicinity in inverted individuals compared to normal ones. We did however notice the absence (read depth = 0) of three loci from the inversion, in the four inverted individuals. The sliding PCA detected either two groups (individuals from Corsica and from mainland) outside of the putative inversion, or four groups at the putative inversion location, separating individuals from Corsica and three groups for mainland individuals (Figure 7b). We assumed that the inverted segment was the one with the lower frequency on the mainland and that was absent from Corsica. The Lostruct analysis confirmed for this same region the presence of four successive blocks of 100 SNPs highly discriminating individuals (Figure 7c). The average F ST between the three groups of mainland individuals was much higher (t test: p < 2.2e-16, Figure 7d) for SNPs in the putative inversion (F ST = 0.26) than for SNPs outside of this region (F ST = 5.8e-4). π was lower (t test: p = 1.9e-5, Figure 7e) for SNPs located in the putative inversion for inverted homozygous individuals (π = 2.5e-5) compared to noninverted homozygous individuals (π = 5.3e-5). Both this inversion compared to heterozygous individuals at this region (π = 7.6e-5) and compared to the three types of individuals for the rest of the chromosome 3 (π ranging from 7.4e-5 to 7.6e-5). The PCA (Figure 7f) and the admixture analysis (Figure 7g) using the SNPs from the inversion clearly delineated inverted, noninverted, and heterozygous individuals, with heterozygous individuals falling at intermediate eigenvalue between the two categories of homozygous individuals. LD was higher (Figure 7h) for the region of the putative inversion (little LD was found when only homozygous noninverted individuals were kept in the analysis). For the region of the inversion, Da was 0.0020 between normal blue tits from the mainland and from Corsica as well as between great tits from the mainland and from Corsica. Da was 0.0044 between inverted and normal blue tits from either the mainland or from Corsica, and ranged from 0.0149 to 0.0152 for interspecific comparisons. Using the formula T = Da/2µ, the inversion likely appeared approximately 220,000 generations ago (around 506,000 years ago). The tree of genetic distance illustrated the divergence of inverted blue tits from noninverted blue tits from both mainland and Corsica (Figure 7i). This region contains 21 genes listed in Table S12. Out of 113 individuals from the mainland, we found 4 inverted homozygous, 23 heterozygous, and 86 noninverted homozygous, and hence, we observed no deviation from the Hardy-Weinberg equilibrium (χ 2 test p-value > .1). The percentage of each genotype was similar across all mainland locations and both habitat types (the inverted segment was at 11% in A1D, 16% in A2D, 17% in A2E, 14% across these three locations).

| Gene ontology

None of the gene lists gathered with the different tests (each outlier test among deciduous and evergreen environment and between

Corsica and the mainland and the inversion test) yielded any significantly enriched GO term after correction multiple testing (Table S14). The most promising GO (uncorrected p-value <.0025) in-

| D ISCUSS I ON

In this study, we investigated demographic history and genomewide patterns of genetic diversity and differentiation between several populations of blue tits presenting strong phenotypic differences between deciduous versus evergreen forest habitats as well as insular versus mainland contexts. Between populations in deciduous and evergreen forest habitats, demographic modeling showed large gene flow and large effective population sizes (Ne), explaining the low genetic differentiation between these populations. Demographic modeling also revealed that each pair of deciduous and evergreen populations most likely never diverged and maintained high connectivity through migration. We identified slight and mostly unrepeated footprints of divergent selection across these evergreen and deciduous population pairs, which is consistent with their demography and the likely polygenic nature of several traits implicated in their local adaptation. In both insular and mainland populations, we found large Ne, although smaller for insular populations than mainland ones, resulting in different distributions and lengths of ROH. Gene flow between Corsica and the mainland most likely stopped at the end of the last glaciation. Both large Ne and gene flow until the end of the last glaciation hence explained moderate genome-wide differentiation. We identified several genomic footprints of selection, enriched in regions of low recombination only in the case of the mainland/island divergence. Finally, we identified a putative genomic inversion spanning 2.8 Mb, polymorphic in mainland populations only. We discuss these results in the context of the vast scientific knowledge acquired on these blue tit populations during the last four decades, and more generally in the context of genomics of adaptation with gene flow.

| Divergence between populations in deciduous versus evergreen habitats

Although we found a significant effect of habitat on genetic structure, the genetic differentiation between neighboring deciduous and evergreen populations was low (F ST ranging from 0.0006 to 0.0079). This result is in line with the primary observations realized earlier on a smaller set of populations [START_REF] Porlier | Habitat-linked population genetic differentiation in the blue tit Cyanistes caeruleus[END_REF][START_REF] Szulkin | Population genomic footprints of fine-scale differentiation between habitats in Mediterranean blue tits[END_REF]. Accordingly, we found high gene flow from deciduous to evergreen populations. Yet, this quantification of high gene flow and low genetic structure contrasted with the demographic knowledge collected on the Blue tit. Indeed, demographic studies suggested restricted dispersal between these populations, with four dispersal events observed between B4D and B4E (5.6 km apart) and none between the B3E and either B4D or B4E (24.1 km), among a total of 2,788 males, 2,672 females, and 25,158 nestling ringed in the three main Corsican sites between 1976 and 2018, with a mean recruitment rate of 6% (Charmantier, com pers). Our interpretation of this contrast between gene flow estimations gathered from population genomic versus recapture data is quadruple.

First, dispersal estimation on the field using capture-mark-recapture is very challenging and may require more data than currently collected, to detect rare dispersal events, even though these affect population genetic parameters. Moreover, since natal dispersal in the Blue tit classically ranges between 330 m and 4 km (see [START_REF] Ortego | Genetic structure reflects natal dispersal movements at different spatial scales in the blue tit, Cyanistes caeruleus[END_REF][START_REF] Tufto | A parametric model for estimation of dispersal patterns applied to five passerine spatially structured populations[END_REF], the long-term monitoring sites in Corsica equipped with nest boxes (black circled dots in Figure 1c)

are not ideally spaced to identify the origin of immigrants and the destination of emigrants, and only a small fraction of the landscape favorable to blue tit breeding is covered. Second, only a few migrants are sufficient to decrease the genetic distance between populations, measured through the F ST [START_REF] Cayuela | Demographic and genetic approaches to study dispersal in wild animal populations: A methodological review[END_REF][START_REF] Marko | The complex analytical landscape of gene flow inference[END_REF].

In that regard, our results may be compatible with the few dispersal events recorded throughout the years. Third, it is important to note that the number of migrants estimated using demographic analyses represents an average over historical time scales on the order of Ne generations and may have varied widely during contemporary times.

Fourth, the large effective population size found both using coalescence and a LD method might be explained by the existence of large "meta-populations" connected by high gene flow and such large Ne

might largely contribute to limit the genetic divergence between populations.

Ne was on average slightly larger in deciduous than in evergreen populations. This could be explained by the higher productivity of deciduous forests compared to evergreen forests, resulting in larger clutches and more fledglings in deciduous habitats (table 1 in [START_REF] Charmantier | Mediterranean blue tits as a case study of local adaptation[END_REF]. However, the very high gene flow between deciduous and evergreen populations and the low genetic differentiation between these populations limits further interpretations. In addition, populations monitored on the long term, for which hundreds of artificial nest boxes have been installed, tended to have larger Ne than nonmonitored populations. Breeding density in the nest box areas studied (Figure 1a) was around 1-1.3 pairs per ha [START_REF] Blondel | A thirty-year study of phenotypic and genetic variation of blue tits in Mediterranean habitat mosaics[END_REF], which is most probably 3-5 times higher than natural densities for blue tits when these secondary hole-nesting birds rely on natural cavities only. It is hence possible that the recent did not originate or build up on pre-existing genetic divergence accumulated during an allopatric phase but rather is supported as a case of genuine recent ecological divergence [START_REF] Pinho | Divergence with gene flow: Models and data[END_REF][START_REF] Wang | Isolation by environment[END_REF].

We did not find strong footprints of divergent selection between any of the four pairs of deciduous and evergreen populations (Figure 4). Moreover, the observed outliers were not found repeatedly across the four deciduous-evergreen comparisons but rather each outlier was found only once, and they were not enriched in regions of low recombination. First, this result is in line with the high migration rates observed. High gene flow indeed most likely limits the potential for local adaptation since it limits the accumulation of allelic differentiation, even with high selection coefficients [START_REF] Lenormand | Gene flow and the limits to natural selection[END_REF]. Second, this pattern is consistent with the best demographic model being the EQ and not a SC. Indeed, the later would more often create repeated and strong outliers located in regions implicated in reproductive isolation between divergent populations (e.g., [START_REF] Rougemont | Inferring the demographic history underlying parallel genomic divergence among pairs of parasitic and nonparasitic lamprey ecotypes[END_REF]. Third, this pattern is congruent with a model of local polygenic adaptation involving a transient genetic architecture with multiple alleles of small effects underlying (multiple) quantitative characters [START_REF] Yeaman | Local adaptation by alleles of small effect[END_REF][START_REF] Yeaman | The genetic architecture of adaptation under migration-selection balance[END_REF] and displaying low F ST among loci under divergent selection. Fourth, most documented traits that are involved in the adaptation of blue tits to the deciduous versus evergreen habitat types, such as clutch size or laying date [START_REF] Blondel | Is the small clutch size of a Corsican blue tit population optimal?[END_REF][START_REF] Lambrechts | A single response mechanism is responsible for evolutionary adaptive variation in a bird's laying date[END_REF], are quantitative and most likely rely on a polygenic architecture with alleles of small effects (see, e.g., [START_REF] Santure | Genomic dissection of variation in clutch size and egg mass in a wild great tit (Parus major) population[END_REF]. However, some statistical issues may limit the interpretation of this result.

First, currently available genome scan methods to detect alleles of small effect are statistically limited, especially when applied to relatively small datasets [START_REF] Hoban | Finding the genomic basis of local adaptation: Pitfalls, practical solutions, and future directions[END_REF][START_REF] Rockman | The QTN program and the alleles that matter for evolution: All that's gold does not glitter[END_REF]. Second, a potential lack of power in detecting outliers could arise from the combination of the use of RAD sequencing and a rapid LD decay along the genome, that may result in too low resolution especially in regions with larger recombination rates [START_REF] Lowry | Responsible RAD: Striving for best practices in population genomic studies of adaptation[END_REF]. Further Charly.

analyses of the potential outliers found here, as well as a new analysis with higher marker density and more individuals, are needed to better document and discuss the potential genes and biological processes implicated. Particularly, quantitative genomics will be useful to establish links between the level of divergent selection on these putative targets and their effect size on phenotypic trait variation [START_REF] Gagnaire | Detecting polygenic selection in marine populations by combining population genomics and quantitative genetics approaches[END_REF][START_REF] Stinchcombe | Combining population genomics and quantitative genetics: Finding the genes underlying ecologically important traits[END_REF]. This would notably reveal whether genes under divergent selection are also responsible for the observed phenotypic variation and contribute to assessing of how much of phenotypic variation is adaptive.

Last, we searched for inversions potentially associated with phenotypic variation and/or segregating in both habitats. We report multiple evidences (Figure 6) for a putative inversion on chromosome 3 in mainland populations. However, the proportion of inversions in the genomes of deciduous-versus evergreen-breeding birds did not differ, which most probably precludes a putative role of the inversion in local adaptation in these habitats. Besides, this polymorphism followed the Hardy-Weinberg equilibrium and we found 4 homozygous inverted individuals, suggesting this inversion does not involve an accumulation of lethal recessive deleterious alleles [START_REF] Jay | Mutation accumulation in chromosomal inversions maintains wing pattern polymorphism in a butterfly[END_REF]. The functional consequences of this inversion on variation in life-history traits in the B lue tit call for further investigations, as achieved in other songbirds [START_REF] Kim | A sex-linked supergene controls sperm morphology and swimming speed in a songbird[END_REF][START_REF] Knief | A sex-chromosome inversion causes strong overdominance for sperm traits that affect siring success[END_REF][START_REF] Tuttle | Divergence and functional degradation of a sex chromosome-like supergene[END_REF]. In order to genotype this inversion in more birds and to link it to putative phenotypic variation, we could use a PCR-RFLP approach (da Silva et al., 2019).

| Divergence between blue tit populations in mainland versus insular contexts

Our demographic modeling approach of Corsican and mainland populations revealed that a model with ancient migration was the most probable, which is coherent with the history previously reconstructed for the blue tit complex. Split time of the ancestral populations in two mainland versus island populations connected by gene flow (~3 M years) was relatively coherent with the diversification time found in the literature for the entire blue tit complex (5M [START_REF] Illera | A multi-gene approach reveals a complex evolutionary history in the Cyanistes species group[END_REF]). Gene flow between the island and the mainland likely stopped around 10,000 years, which is compatible with a gene flow break due to the rise of the flooded stretch between Corsica and the mainland during ice melting and the sea level rise after the last glaciation, during from 17,000 to 5,000 years ago [START_REF] Jouet | Shoreface migrations at the shelf edge and sea-level changes around the Last Glacial Maximum (Gulf of Lions, NW Mediterranean)[END_REF][START_REF] Lambeck | Sea-level change along the French Mediterranean coast for the past 30 000 years[END_REF]. Both the signal of population expansion during the diversification period and the gene flow five times larger from the island to the mainland than from the mainland to the island were coherent with a suspected recolonization of the mainland from the islands, as proposed by [START_REF] Illera | A multi-gene approach reveals a complex evolutionary history in the Cyanistes species group[END_REF] based on analyses of nuclear and mitochondrial DNA sequences. The larger effective population size on the mainland than on the island may be due to multiple sources of colonization on the mainland [START_REF] Taberlet | Unusual mitochondrial DNA polymorphism in two local populations of blue tit Parus caeruleus[END_REF] and to the much larger mainland area and hence meta-population size compared to the island. Expanding this study with a sampling from the South-East of Europe, one from Sardinia and one from the Iberian Peninsula, could help determine whether the large size estimated for the mainland meta-population studied here might be explained by a mixed recolonization from these refugium [START_REF] Kvist | The colonization history and present-day population structure of the European great tit (Parus major major)[END_REF][START_REF] Kvist | Glacial history and colonization of Europe by the blue tit Parus caeruleus[END_REF] and whether insular populations consistently have reduced effective population size. It would also be interesting to explore a model allowing for multiple cycles of isolations and contacts that could have resulted from successive glacial cycles [START_REF] Hewitt | Genetic consequences of climatic oscillations in the Quaternary[END_REF].

Given the demographic parameters inferred, it is not surprising that the genetic differentiation between populations from the mainland and the island was moderate (F ST = 0.08). However, this differentiation value contrasted with the much larger divergence at mitochondrial DNA observed for these populations (OST = 0.67 between A1D and B4D [START_REF] Kvist | Glacial history and colonization of Europe by the blue tit Parus caeruleus[END_REF]). This discrepancy between whole nuclear genome and mitochondrial estimates of differentiation could first be explained by the typically four times smaller Ne in mitochondrial DNA compared to nuclear DNA in diploid organisms, resulting in stronger drift and divergence [START_REF] Smith | Examining the role of effective population size on mitochondrial and multilocus divergence time discordance in a songbird[END_REF].

Second, accumulation of cytonuclear incompatibilities could limit its introgression [START_REF] Burton | A disproportionate role for mtDNA in Dobzhansky-Muller incompatibilities[END_REF]. Alternatively, gene flow may have been greater for males than females during the colonization process from Corsica to the mainland, at the end of the last glacial period. However, this hypothesis is not supported by the general observation that females disperse on average further away than males and that rare long-distance dispersals are also achieved by females [START_REF] Tufto | A parametric model for estimation of dispersal patterns applied to five passerine spatially structured populations[END_REF]. Inspecting mitochondrial DNA variation for the birds studied here would be very useful to compare mitochondrial differentiation in our sample to the previous analysis which also included Corsican birds [START_REF] Kvist | Glacial history and colonization of Europe by the blue tit Parus caeruleus[END_REF]. In any case, one should bear in mind that our results stem from integrating coalescent patterns observed across thousands of loci, therefore providing increased resolution to investigate the determinants of demographic divergence compared to approaches based on fewer mitochondrial and nuclear loci.

The contemporary demography inferred between the populations from Corsica and mainland France, especially the absence of gene flow for approximately 4,500 generations and the large effective population size, likely provides a suitable context for the build-up of local adaptation. Moreover, the level of divergence (F ST = 0.08) was probably low enough to detect divergent outlier loci potentially implicated in local adaptation. We identified several outlier SNPs and outlier 200 kb windows showing elevated differentiation, more often in regions of lower recombination than elsewhere in the genome (contrary to the outliers found for the deciduous vs. evergreen populations). Furthermore, F ST was on average higher in these regions of low recombination compared to elsewhere in the genome (which was not the case for deciduous vs. evergreen populations). These patterns of increased divergence in regions with low recombination between divergent populations are commonly observed in other species [START_REF] Cruickshank | Reanalysis suggests that genomic islands of speciation are due to reduced diversity, not reduced gene flow[END_REF][START_REF] Gagnaire | Analysis of genome-wide differentiation between native and introduced populations of the cupped oysters Crassostrea gigas and Crassostrea angulata[END_REF][START_REF] Nachman | Recombination rate variation and speciation: Theoretical predictions and empirical results from rabbits and mice[END_REF], including other bird species [START_REF] Burri | Linked selection and recombination rate variation drive the evolution of the genomic landscape of differentiation across the speciation continuum of Ficedula flycatchers[END_REF][START_REF] Spurgin | The great tit HapMap project: A continental-scale analysis of genomic variation in a songbird[END_REF]. It is well documented that such increased differentiation in regions with low recombination is not necessarily due to positive selection, or at least not alone, and that it is largely influenced by the effect of recombination in interaction with background selection [START_REF] Burri | Linked selection and recombination rate variation drive the evolution of the genomic landscape of differentiation across the speciation continuum of Ficedula flycatchers[END_REF][START_REF] Charlesworth | The effect of deleterious mutations on neutral molecular variation[END_REF][START_REF] Perrier | On the importance of time scales when studying adaptive evolution[END_REF]. Moreover, a uniform F ST sliding window, sized in Kb, is expected to dilute signatures of selection in regions of the genome where LD is fast decaying while it is expected to over-represent regions where LD decays much more slowly, hence increasing even more the potential false positives in regions with low recombination [START_REF] Beissinger | Defining window-boundaries for genomic analyses using smoothing spline techniques[END_REF][START_REF] Perrier | On the importance of time scales when studying adaptive evolution[END_REF]. Although perfectible, our sliding window approach using a SNP unit attempted to fix this issue and successfully captured several new outlier regions outside of deserts of recombination. This second window approach should however be improved, for example, by estimating the local neutral enveloped, which computation should integrate local variations in LD, recombination, diversity, but also the demography of populations. As mentioned for the study of adaptive divergence between deciduous and evergreen populations, complementary analyses integrating both genome scans and quantitative genomics would improve our comprehension of the genomic and phenotypic divergence observed between Corsican and mainland blue tits [START_REF] Gagnaire | Detecting polygenic selection in marine populations by combining population genomics and quantitative genetics approaches[END_REF][START_REF] Stinchcombe | Combining population genomics and quantitative genetics: Finding the genes underlying ecologically important traits[END_REF].

While the putative genomic inversion on chromosome 3 was absent from Corsica and detected in mainland individuals, its level of divergence from the noninverted sequence indicated that it was likely twice older than the beginning of divergence between blue tit populations from mainland and Corsica. This could first suggest that this polymorphism emerged in mainland blue tit populations and then did not introgress the Corsican populations, maybe due to a local disadvantage and/or genetic incompatibilities or simply due to drift coupled to little gene flow. A second hypothesis could be that this inversion was present in Corsican populations but had been purged out due to a local disadvantage. Such a disadvantage can be due to the typical accumulation of deleterious mutations in such nonrecombining inverted sequences (e.g., [START_REF] Jay | Mutation accumulation in chromosomal inversions maintains wing pattern polymorphism in a butterfly[END_REF].

Lastly, it is possible that this inversion has been acquired in mainland populations only recently, after the last period of contact between mainland and Corsican populations, via gene flow from a distinct refugium with which connectivity could have been enhanced after deglaciation, hence after the stop of gene flow with Corsica. The age of this inversion, its origin, its biological effects, and the potential accumulation of deleterious mutations need to be inferred more thoroughly via genotyping full inversion sequences from individuals in diverse locations and using advanced statistical methods [START_REF] Lohse | Genome-wide tests for introgression between cactophilic Drosophila implicate a role of inversions during speciation[END_REF].

| CON CLUS I ON AND PER S PEC TIVE S

Our study demonstrated the usefulness of demographic modeling and of the analysis of the variation of genomic diversity and recombination along the genome to uncover the genetic determinants of local adaptation in a small passerine with a large distribution, and occupying different forest habitats. Especially, demographic modeling rejected the hypothesis of a secondary contact between deciduous and evergreen populations and favored a situation with continuous gene flow. These results support the idea that blue tits have adapted to their habitats despite ongoing gene flow, while contextualizing how large gene flow most probably constrained local adaptation [START_REF] Lenormand | Gene flow and the limits to natural selection[END_REF] and favored its architecture based on alleles of small effects [START_REF] Yeaman | Local adaptation by alleles of small effect[END_REF]. The genomic modeling also refined our knowledge about the divergence between insular and mainland meta-populations, that have been likely unconnected by gene flow for the last ten thousand years. We also verified the relationship between local recombination rate and differentiation, that should probably be integrated in genome scans looking for footprints of selection [START_REF] Beissinger | Defining window-boundaries for genomic analyses using smoothing spline techniques[END_REF][START_REF] Berner | Genomics of adaptive divergence with chromosome-scale heterogeneity in crossover rate[END_REF][START_REF] Booker | Variation in recombination rate affects detection of FST outliers under neutrality[END_REF][START_REF] Burri | Linked selection and recombination rate variation drive the evolution of the genomic landscape of differentiation across the speciation continuum of Ficedula flycatchers[END_REF][START_REF] Perrier | On the importance of time scales when studying adaptive evolution[END_REF].

Future investigations will require increased sample sizes and marker density [START_REF] Hoban | Finding the genomic basis of local adaptation: Pitfalls, practical solutions, and future directions[END_REF][START_REF] Lotterhos | The relative power of genome scans to detect local adaptation depends on sampling design and statistical method[END_REF] in order to better detect loci with small effects that contribute to the quantitative phenotypic variation and local adaptation in blue tits. Lastly, the putative inversion found here would need further analyses since this type of structural variation is often implicated in phenotypic variation [START_REF] Kirkpatrick | Chromosome inversions, local adaptation and speciation[END_REF][START_REF] Kirkpatrick | How and why chromosome inversions evolve[END_REF][START_REF] Stapley | Variation in recombination frequency and distribution across eukaryotes: Patterns and processes[END_REF][START_REF] Wellenreuther | Going beyond SNPs: The role of structural genomic variants in adaptive evolution and species diversification[END_REF].

ACK N OWLED G EM ENTS

We wish to thank Louis Bernatchez for inviting our contribution in this special issue (Box 1). We wish to thank the many people who helped on the blue tit long-term project. We particularly wish to thank Philippe Perret. We also wish to thank Samuel Perret, Christophe De

France

  Abbreviation

(

  PAN), a model of equilibrium corresponding to the island model with two populations (EQ), a model of isolation with migration (IM), a model of SC, a model of divergence with migration during the first generations, that is, ancestral migration (AM), and a model of strict isolation (SI). The prior and details of the simulation pipeline are fully described in the Note S1. The PAN model assumes that the two focal populations descent from a single panmictic population characterized by its effective size (Ne). The EQ model (equivalent to the island model) assumes that the population is subdivided into two discrete populations of sizes N pop1 and N pop2 that are connected by continuous gene flow at a constant rate each generation. In this model, the divergence time is not a parameter. The IM, SI, SC, and AM models all assume that an ancestral population of size N ANC splits at T split into two daughter populations of sizes N pop1 and N pop2. Under SI, no subsequent gene flow occurs.

  to search for SNPs potentially under divergent selection at two geographic levels: (a) between the mainland and Corsican populations and (b) between each local pairs of evergreen (E) and deciduous (D) populations, that is, A2D versus A2E, B4D versus B4E, B5D versus B5E, B6D versus B6E.

  lelism across pairs of D-E environments. Second, we estimated F ST along the genome using either a 200 kb sliding average with VCFtools or the function "snpgdsSlidingWindow" from the package SNPrelate to estimate F ST in blocks of 50 SNPs moving by 5 SNPs. The second window strategy was used in order to compensate for lower SNP density in regions of low recombination that tend to exaggerate the contribution of individual SNPs in these regions and to dilute the individual SNP contribution in regions of high recombination (see Perrier & Charmantier, 2019 for a broader comment on this).Third, we used an RDA as an alternative method to search for SNPs putatively implicated in multilocus adaptation (a) between populations in deciduous and evergreen habitats, and (b) between populations in Corsica and on the mainland. As suggested by[START_REF] Forester | Comparing methods for detecting multilocus adaptation with multivariate genotype-environment associations[END_REF], such a multivariate method may be more suitable than univariate ones to detect weaker footprints of adaptation that are expected in polygenic adaptations in response to complex environmental heterogeneity. Using a similar procedure as described earlier in the methods, we used two RDAs constrained to investigate the effect of phenology (a) or of the geography (b). We then used a three standard deviation cutoff as suggested byForester et al. (

  (i.e., E vs. D, p = .006), latitude (p < .001), longitude (p < .001), and altitude (p < .001). The first axis of the global RDA explained 4.54% of the variance and was correlated mainly with latitude (-0.987) and longitude (+0.999). The fourth axis of this global RDA displayed the strongest correlation with phenology (-0.856) and explained 0.27% of the variance. The partial RDA model conditioning for the effects of latitude, longitude, and altitude was globally significant (p < .001), as well as the effect of phenology (p < .001). The axis of the partial RDA explained 0.29% of the variance and was correlated with phenology (0.922).

  For the AM model, the simplest model without linked selection or barriers to gene flow was the best supported, with p = .76. Parameters estimates revealed posterior generally well differentiated from the prior providing increased confidence (Figure S4). Effective population size of the mainland was 397,685 [CI = 244,375-600,105], around 7.5 times greater than the one of the island, 52,990 [CI = 33,065-78,475] (Figure 3b, F I G U R E 2 Population structure between populations. (a) Admixture for K = 2, 6 and 10 (Vertical lines correspond to individual admixture; colors correspond to distinct genetic groups). (b) RDA

  indicated a strong population size change during the process of population isolation, since the ancestral Ne was estimated at only 240 [CI = 65-685]. Split time was estimated at 1,417,480 [CI = 740,[0][1][START_REF] Booker | Variation in recombination rate affects detection of FST outliers under neutrality[END_REF]394,000] generations ago, hence around 3.2 Ma (assuming a 2.3 years generation time[START_REF] Charmantier | Testing for microevolution in body size in three blue tit populations[END_REF]) and gene flow subsequently stopped 4,540 [CI = 2,500-6360] generations ago, hence around 10,000 years ago, corresponding closely with the end of the last glacial maximum. The AM rate between Corsica and the mainland, from the split in two populations to the end of the latest gene flow between these populations, was inferred asymmetric, with a five times larger migration rate from the island to the mainland (m ~ 1.9e-4) than the other way around (m ~ 3.8e-5). However, this asymmetry did not generate a strong difference in gene flow between the two groups, given the lower effective population size on the island compared to the mainland, (i.e., respective number of migrants of 10 and 15 from Corsica to the mainland and in the reverse direction, with overlapping credible intervals, see details in Table

Finally, we observed

  significantly slightly longer ROH on Corsica compared to mainland (mean island = 44,769 kb vs. mean mainland = 41,586 kb, Wilcoxon rank sum test W = 32,695, p < 2e-16, Figures S7&S8).

Figures 4

 4 Figures 4 & 5). Throughout the five Bayescan tests (Corsica vs. mainland populations and the 4 D-4 population pairs), we identified 40 SNPs with a log10(BF)> 0 among which 18 SNPs had a log10(BF)> 1

F I G U R E 4 F

 4 ST Manhattan plots between the four population pairs in deciduous and evergreen forests. F ST are given either by SNP (four first graphs), by 200 kb windows (graphs 5 & 6), or by 50 SNPs windows (graphs 7 & 8). For graphical simplicity, the three pairs of populations in Corsica were grouped for the sliding window F ST Manhattan plots (BD vs. BE). Dots alternate colors between chromosomes. F ST for unplaced loci are shown at the end of each of the first four Manhattan plots inverted and noninverted homozygous individuals had lower π at

  cluded the following: (a) for the deciduous-evergreen outlier tests, GO:0060385, axonogenesis involved in innervation; GO:2001013, epithelial cell proliferation; GO:0007194, negative regulation of adenylate cyclase; GO:0090647, modulation of age-related behavioral; GO:1901379, regulation of potassium ion transmembrane transport, (b) for the Corsica-mainland outlier tests, GO:0008344, adult locomotory behavior; GO:2001224, positive regulation of neuron migration, (c) for the inversion, GO:0016446, somatic hypermutation of immunoglobulin; GO:0045910, negative regulation of DNA recombination; GO:0006298, mismatch repair; GO:0008340, determination of adult lifespan.

F I G U R E 5 F

 5 ST Manhattan plots between Corsica and the mainland. F ST are given either by SNP, by 200kb windows or by 50 SNPs windows. Dots alternate colors between chromosomes. F ST for unplaced loci are shown at the end of the first Manhattan plot

F I G U R E 6

 6 Relationship between recombination rate and divergence. (a) Correlation between local recombination rate and F ST between populations on the mainland and in Corsica; (b) SNP frequency per class of recombination rate; (c) outlier SNP frequency for the mainland-Corsica comparison, per class of recombination rate; (d) outlier SNP frequency for the deciduous-evergreen comparisons, per class of recombination rate. In panels c & d, outliers are given for the Bayescan and the RDA tests

F

  Detection and characteristics of a putative inversion on chromosome 3. (a) Recombination rate and SNP density along chromosome 3; (b) individual values for eigenvector 1 of a PCA of all the blue tit individuals from Corsica and the mainland; (c) MDS 1 values from Lostruct analysis;(d) F ST between inverted, noninverted, and heterozygous individuals from the mainland; (e) nucleotide diversity for inverted (in red), noninverted (in green), and heterozygous individuals (in blue) from the mainland; (f) PCA results for the region of the putative inversion for all the blue tit individuals; (g) mainland blue tits' admixture proportions for K = 2 for the region of the putative inversion (individual ordered by decreasing ancestry value for the inverted cluster); (h) linkage disequilibrium for SNPs found in the putative inversion (plus 1 Mb each side of the inversion), for mainland blue tits; (i) genetic distances at the inversion, between four inverted blue tits, eight noninverted blue tits, and four great tits (sites of capture are indicated for each individual) availability of nest boxes locally increased productivity, effective population size, and heterozygosity (Box 1).Demographic modeling indicated that secondary gene flow between previously isolated populations, which could have eroded most genome-wide differentiation outside regions implicated in local adaptation and/or acting as reproductive barriers in particular in low recombining areas[START_REF] Bierne | Pervasive selection or is it…? why are FST outliers sometimes so frequent?[END_REF], was very unlikely. Secondary contacts were strongly rejected in favor of an island model supporting ongoing gene flow without divergence (EQ model) in three out of four population pairs, the last one supporting an isolation with migration model (IM). As is often the case, our demographic estimates should be interpreted cautiously, particularly in the absence of a well-documented mutation rate in the Blue tit. The relatively narrow credible intervals around parameter estimates nevertheless provide confidence in their biological relevance. These results suggest that the focal population pairs occupying deciduous/evergreen habitats have been continuously connected by relatively high gene flow. Hence, genetic divergence between populations in each pair

Franceschi,

  Jacques Blondel, Marta Szulkin, Monica Arias, Patricia Sourouille, Annick Lucas & Boris Delahaie. We thank the associate editor, two reviewers, and the production team for their constructive comments and edits on previous versions of this manuscript. This project was funded by the European Research Council (starting grant ERC-2013-StG-337365-SHE to AC) and the OSU-OREME. study, conducted fieldwork, laboratory analyses, bioinformatics and statistical analyses, interpreted the results and wrote the manuscript. QR conducted bioinformatics and statistical analyses, interpreted the results and edited the manuscript. AC obtained the ERC grant, conceived the study, performed fieldwork, interpreted the results and edited the manuscript. DATA AVA I L A B I L I T Y S TAT E M E N T Individual sequences with metadata are available on NCBI with the SRA accession PRJNA630135. A filtered VCF file with metadata are available on Dryad with the https://doi.org/10.5061/dryad.x69p8 czfg.

Table 2

 2

) than Island populations (on average 0.147), all t tests' p-values <.001. On the contrary, there was no significant difference of observed heterozygosity at smaller spatial scale among populations within the island (either same or different tree phenology, E-and D-) or within the mainland (all t tests' p-value >.05). More SNPs were polymorphic in mainland populations

(77,616

SNPs for a MAF > 0.05) than in island populations (62,728-65,876 SNPs for a MAF > 0.05). The MAF spectrum showed enrichment of variants with smaller frequencies in mainland populations versus island populations (Figure

  Diversity and differentiation across populations. Fixation index between populations (F ST ), effective population size (Ne) estimated with ABC (Ne-ABC) using models for deciduous-evergreen population pairs or for Corsica-mainland (A1D-B4D), Ne estimated with a LD method (Ne-LD), and observed heterozygosity (Ho). Color gradients indicate differentiation intensity (Darker red for higher differentiation and darker blue for lowest differentiation). Underlined values indicate deciduous and evergreen pairs

			Ne-LD Ho	355 0.171	161 0.167	142 0.168	334 0.149	310 0.149	271 0.148	210 0.145	216 0.146	268 0.146	260 0.146
		Corsica-	Mainland	397,685	na	na	Na	52,990	na	na	na	na	na
	Ne-ABC	Deciduous-	Evergreen	Na	400,005	100,055	Na	43,070	37,100	46,520	42,800	46,885	162,940
			B6E	0.0784	0.0792	0.0787	0.0058	0.0062	0.0085	0.0032	0.0033	0.0006
			B6D	0.0785	0.0793	0.0790	0.0056	0.0061	0.0083	0.0035	0.0036		0.0006
			B5E	0.0804	0.0810	0.0807	0.0071	0.0067	0.0097	0.0025		0.0035	0.0033
			B5D	0.0796	0.0799	0.0796	0.0067	0.0063	0.0093		0.0025	0.0035	0.0032
			B4E	0.0853	0.0872	0.0871	0.0089	0.0079		0.0093	0.0097	0.0083	0.0085
			B4D	0.0830	0.0849	0.0846	0.0058		0.0079	0.0063	0.0067	0.0061	0.0062
			B3E	0.0812	0.0832	0.0828		0.0058	0.0089	0.0067	0.0071	0.0056	0.0058
			A2E	0.0043	0.0029		0.0828	0.0846	0.0871	0.0796	0.0807	0.0790	0.0787
			A2D	0.0045		0.0029	0.0832	0.0849	0.0872	0.0799	0.0810	0.0793	0.0792
	F ST		A1D		0.0045	0.0043	0.0812	0.0830	0.0853	0.0796	0.0804	0.0785	0.0784
				A1D	A2D	A2E	B3E	B4D	B4E	B5D	B5E	B6D	B6E

observed 1.18 times lower allele frequencies for nonsynonymous variants (average MAF = 0.09) than for synonymous ones (average MAF = 0.11) (t test p-value < .000001). This overall lower frequency of nonsynonymous compared to synonymous mutations was similar between populations from Corsica (1.19) and the mainland (1.16), as well as between deciduous (1.18) and evergreen (1.19) habitats (see TA B L E 2

Box 1 Personal (Charles Perrier) reflections on my career and my collaboration with Louis Bernatchez

	These are personal reflections regarding my early career at
	Louis Bernatchez lab. This box is hence not related to the
	population genomics project on blue tits.
	I have many ideas to write a very scientific comment about
	my collaboration with Louis and how it influenced my
	young carrier… but in those times of confinements due to
	COVID-19 I am rather willing to discuss more personal and
	social aspects that also matter a lot in academia. During
	my stay in Louis laboratory, I admired especially two of
	Louis talents. The first one, maybe unsurprising, the out-
	standing leader that Louis is and the positive effect that
	this has on getting the job done on his laboratory. He is
	always enthusiastic about new ideas, permanently push-
	ing you to pursue your scientific objectives, and constantly
	finding the bright side of every hard work done to get it
	published. I will not detail here how it helped me to grow
	as a scientist but it surely did in many ways! To my opin-
	ion, Louis leadership transcends his indisputable academic
	performance and reaches other important values, namely
	equity and diversity. And very importantly for me, Louis
	stands up and fight when necessary. I will give a small ex-
	ample (hoping I am authorized to): While receiving a price
	for his hard work from the Quebecois government, Louis
	had to give a speech. Instead of seeking to cast even more
	light on his own glory, he gave an engaged speech defend-
	ing the long standing equitable and open university prin-
	ciples that the government of Quebec was threatening at
	that time by liberal reforms including the rising of tuition
	fees. Louis basically claimed that the previous system was
	the one that allowed many guys like him, coming from
	relatively modest families to obtain a university degree,
	to pursue an academic career in science, and to bring the
	diversity of perspectives the university and the society
	need… Hence arguing that such price was not only reward-
	ing him but also a relatively open access to the university. I
	found that speech extremely appropriate and inspiring. In
	such good times as that ceremony for Louis' price, as well
	as in bad times as the current covid19 epidemic, I find par-
	ticularly important to remember that societal and scientific
	"advances" or "progress" could or should also (or only?) be
	measured by their contribution to human equity, diversity,
	and fraternity. Louis, keep up doing good science, motivat-
	ing young scientists, and fighting for what you believe is

right for fish, science, and people! See you soon for a fishing trip in Camargue and "un café double"! Happy birthday!

Funding information

European Research Council, Grant/Award Number: ERC-2013-StG-337365-SHE