C Perrier
email: charles.perrier@cefe.cnrs.fr

Heritability estimates from genome wide relatedness matrices in wild populations: application to a passerine, using a small sample size Running head: GRM in Corsican blue tits

Keywords: heritability, genetic correlation, RAD-sequencing, GRM, pedigree, partitioning, phenotype, blue tit, SNP, Cyanistes caeruleus

    Heritability estimates from genomewide relatedness matrices in wild populations: Application to a passerine, using a small sample size C. Perrier, B. Delahaie, A.

Introduction

Estimating additive genetic variance and heritability of quantitative traits is a major aim in evolutionary biology because they are crucial components of all evolutionary models. Further, fine-scale dissection of the genomic basis of quantitative traits in the wild, using methods such as genome wide association and linkage mapping, requires validation that the traits being mapped are actually heritable. New sequencing technologies have empowered the estimation of genomic relatedness matrices (GRM) between individuals [START_REF] Visscher | Heritability in the genomics era -concepts and misconceptions[END_REF][START_REF] Huisman | Pedigree reconstruction from SNP data: parentage assignment, sibship clustering and beyond[END_REF][START_REF] Gienapp | Genomic quantitative genetics to study evolution in the wild[END_REF], feasibly without any social pedigree (PED-social) recorded in the field. Since in a quantitative framework, common SNPs additively explain a large part of genetic additive variance, it is possible to use such realized relatedness estimated from genomic data to measure traits' heritability in wild populations (Yang et al. 2010). Furthermore, by aligning markers against a reference genome, heritability can then be partitioned among chromosomes (or at other genomic scales) enabling testing quantitative models of increasing heritability with the number of genes per chromosomes (Yang et al. 2011b), or alternatively, identifying chromosomes displaying higher genetic variance than expected. GRM estimates and chromosome partitioning, first applied on human data (Yang et al. 2010;2011b), were then applied in well-known free-ranging animal populations (Great tits, Parus major, [START_REF] Robinson | Partitioning of genetic variation across the genome using multimarker methods in a wild bird population[END_REF]; Soay sheep, Ovis aries, [START_REF] Bérénos | Heterogeneity of genetic architecture of body size traits in a free-living population[END_REF] for which thousands of individuals were genotyped at thousands of loci. There have been fewer empirical applications of these methods to smaller samples (but see [START_REF] Wenzel | Genome-wide association and genome partitioning reveal novel genomic regions underlying variation in gastrointestinal nematode burden in a wild bird[END_REF][START_REF] Silva | Insights into the genetic architecture of morphological traits in two passerine bird species[END_REF], although they have great potential in the current context of genomic tools being increasingly applied to natural populations [START_REF] Gienapp | Genomic quantitative genetics to study evolution in the wild[END_REF]) that often rely on relatively small samples.

Three important limitations for estimating GRM and GRM based heritability are the number of SNPs, the number of individuals and the sampling variance [START_REF] Visscher | A General Unified Framework to Assess the Sampling Variance of Heritability Estimates Using Pedigree or Marker-Based Relationships[END_REF]. The limitation from the number (as well as diversity and quality) of markers and their ability to resolve relatedness essentially vanished when we entered the genomic era [START_REF] Csillery | Performance of Marker-Based Relatedness Estimators in Natural Populations of Outbred Vertebrates[END_REF][START_REF] Stanton-Geddes | Estimating heritability using genomic data (J Hadfield[END_REF][START_REF] Gienapp | Genomic quantitative genetics to study evolution in the wild[END_REF]. The sample size of individuals genotyped and phenotyped, however, remains an issue. For a given dataset, the number of individuals included will influence both error and median values of heritability estimates, as tested in two studies [START_REF] Stanton-Geddes | Estimating heritability using genomic data (J Hadfield[END_REF][START_REF] Visscher | A General Unified Framework to Assess the Sampling Variance of Heritability Estimates Using Pedigree or Marker-Based Relationships[END_REF]. Nevertheless, reasonable heritability estimates can be reached using 150 to 200 individuals, as showed in a study on barrelclover (Medicago trunculata, [START_REF] Stanton-Geddes | Estimating heritability using genomic data (J Hadfield[END_REF]. However, different species, demographic contexts, and sampling strategies will likely perform differently for a similar number of individuals, number of markers and traits' architecture. [START_REF] Visscher | A General Unified Framework to Assess the Sampling Variance of Heritability Estimates Using Pedigree or Marker-Based Relationships[END_REF] showed that for designs that use genetic markers to estimate relatedness among randomly sampled individuals from a population (which may be a common situation for non-model species in natural populations), the error in heritability estimation is inversely proportional to the squared sample size and is proportional to the effective population size. Therefore, heritability estimates based from genetic marker relatedness in extremely large populations will likely require thousands of samples (e.g. marine pelagic fish [START_REF] Gagnaire | Detecting polygenic selection in marine populations by combining population genomics and quantitative genetics approaches[END_REF]), but also thousands of loci (given the small linkage disequilibrium among loci). Nevertheless, for populations of smaller size, a sampling strategy aiming at capturing variance in realized relatedness offers possibilities to obtain robust heritability estimates bounded by reasonable standard errors.

While genomics offers promising potential to obtain quantitative genetic parameters for virtually any species or population [START_REF] Gienapp | Genomic quantitative genetics to study evolution in the wild[END_REF], it has only been used in a few studies to-date and more applications in smaller samples and non-model species are needed. The first studies to estimate genomic heritability and partition variation across the genome elsewhere than in humans were, understandably, focusing on already well-known populations and species. In particular, analysis of 2644 individuals at 7203 SNPs from a long term population study of great tits was the first application of quantitative genomics to a wild population [START_REF] Robinson | Partitioning of genetic variation across the genome using multimarker methods in a wild bird population[END_REF]. The Soay sheep (Ovis aries) also benefited from such a genomic quantitative genetic study, with 5805 individuals genotyped at 37,037 autosomal SNPs [START_REF] Bérénos | Estimating quantitative genetic parameters in wild populations: a comparison of pedigree and genomic approaches[END_REF]). This study revealed that most of the additive genetic (co)variances in sheep body size traits were captured by half of the SNPs. These pioneering studies enabled the verification of the power of quantitative genomics on known populations with long-term pedigrees and characterized quantitative genetic parameters [START_REF] Edwards | Next-generation QTL mapping: crowdsourcing SNPs, without pedigrees[END_REF]. Other authors also estimated heritability from GRM in datasets with much less individuals, e. These examples confirmed the usefulness of the method in new cases and with fewer individuals. Although we did not perform a thorough meta-analysis, there was a trend among these studies for an increase, among and within studies, of heritability estimates' standard errors with decreasing sample size and SNP amount. Nevertheless, and surprisingly considering the number of studies performing genomic analyses since the beginning of the genomic era, there are relatively few applications of genomic data to estimate heritability, genetic correlations and chromosome partitioning, using smaller datasets.

Here, we assessed heritability, genetic correlations, and chromosome partitioning for four phenotypic traits, namely tarsus length, wing length, bill length, and body mass, in a collection of 494 individuals of Corsican blue tits (Cyanistes caeruleus ogliastrae) from three closely located sites in Corsica (Figure 1). Based on mitochondrial genetic and phenotypic data, blue tits in Corsica and Sardinia have been qualified as a subspecies, Cyanistes caeruleus ogliastrae [START_REF] Kvist | Glacial history and colonization of Europe by the blue tit Parus caeruleus[END_REF]. The focal populations are ideal for quantitative genetic study because of the availability of pedigree, phenotypic and genetic data gathered through a long-term study [START_REF] Charmantier | Mediterranean blue tits as a case study of local adaptation[END_REF]. Based on microsatellite and SNP genetic data, important gene flow together with significantly low genetic differentiation have been suggested between Corsican 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]. Phenotypic variance for the four aforementioned traits has been characterized as quantitative and moderately to highly heritable (Charmantier et al. 2004b;[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 latest published results based on 17 to 27 years of monitoring, provided population-specific significant heritability ranging from 0.43 to 0.57 for tarsus length, 0.20-0.33 for wing length, 0.18-0.34 for bill length, and 0.22-0.32 for body mass. All trait combinations showed significant additive genetic covariance (COVA) in at least one population; in particular COVAs between wing length and tarsus length and between wing length and body mass were significantly positive [START_REF] Delahaie | Conserved G-matrices of morphological and life-history traits among continental and island blue tit populations[END_REF]. Lastly, from 18 to 25% of young have been identified as extra-pair offspring in these Corsican blue tit populations (Charmantier et al. 2004a), which may lead to slight underestimation of trait heritability from a fieldbased social pedigree [START_REF] Charmantier | How do misassigned paternities affect the estimation of heritability in the wild?[END_REF].

The general objective of this study was to test the effectiveness of GRM to estimate heritability and genetic correlations among traits and to perform heritability partitioning while using a relatively small sample size of Corsican blue tits. We produced a dataset of 494 Corsican blue tits, from three sites (sample sizes of 110, 185 and 199), genotyped at nearly 50,000 filtered SNPs derived from RAD-seq and phenotyped for four quantitative phenotypic traits (tarsus length, wing length, bill length, and body mass). We looked for erroneous links in the social pedigree (PED-social) using the genetic data, in order to build a corrected social pedigree (PED-corrected). We estimated heritability using both GCTA (Yang et al. 2011a) and MCMCglmm [START_REF] Hadfield | MCMC methods for multi-response generalized linear mixed models: the MCMCglmm R package[END_REF] for the 494 genotyped individuals (with the three sites pooled or separately), based on GRM, PED-social and PED-corrected. Heritability differences between PEDcorrected and PED-social were computed in order to infer the possible effects of pedigree errors on estimates of traits heritability (notably underestimation originating from extra-pair copulations [START_REF] Charmantier | How do misassigned paternities affect the estimation of heritability in the wild?[END_REF]. Comparisons between heritability estimated using PED-corrected and GRM aimed at defining the effect of the greater precision and information content of the GRM. We estimated the effect of the number of SNPs on tarsus length heritability estimates, in order to discuss the SNP density needed to recover much of the additive genetic variance. Finally, we estimated genetic correlations between traits and we partitioned heritability between chromosomes, although these quantitative genomic tests may challenge our sample size and probably require much larger sample size to produce accurate results. We discuss the usefulness and limitations of these heritability inferences in small samples from wild populations of non-model species.

Methods

Study sites, bird monitoring and phenotypic data

Data were collected in three locations (E-PIRIO, D-MURO and E-MURO) in Corsica (France, Figure 1, Table 1). Approximately six kilometers separate D-MURO from E-MURO and 27 km separate these sites from E-PIRIO. The landscape of these sites is dominated either by the evergreen holm oak Quercus ilex (E-PIRIO and E-MURO) or by the deciduous downy oak Q. pubescens (D-MURO). All sites were monitored as part of a long-term research programme for 20 (E-Muro), 24 (D-Muro) and 42 years (E-Pirio) until 2017. A total of 360 nest-boxes were monitored and 26,650 birds were ringed. Each year, nest boxes were visited at least once a week during the reproductive period (from early April to late June). Breeding blue tits were captured in nest boxes during the feeding of their young, and banded (if not already banded earlier) with a unique metal ring provided by the CRBPO. Nestlings were also banded before fledging at 9-15 days old.

Four phenotypic traits were measured on male and female breeders: tarsus length (from the intertarsal joint to the most distal undivided scute on the tarsometarsus), flattened wing length, bill length (from the anterior end of the nares to the tip of the upper mandible), and body mass. Phenotypic differences among individuals from the different sites and used in genomic heritability estimates are presented in Table 1 and Supplementary material 1 (see [START_REF] Charmantier | Mediterranean blue tits as a case study of local adaptation[END_REF] for a review of phenotypic divergence in these sites). Five to 20µl of blood were sampled from a small neck vein or from a wing vein from adult breeders for later DNA extraction. 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].

DNA extraction, RAD-sequencing, SNP calling and data filtering

We selected a set of 494 individuals captured between the years 2010-2016, from which we extracted DNA from collected blood. The individuals were captured from three sites (E-PIRIO, n = 185; D-MURO, n = 199; E-MURO, n = 110; Table 1). Individuals were chosen according to the presence of several phenotypic measurements and a blood sample. DNA extraction was achieved using Qiagen DNeasy Blood & Tissue kits. DNA extractions 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 checked on agarose gels.

Library preparation using restriction-site-associated DNA sequencing (RAD-seq; [START_REF] Baird | Rapid SNP Discovery and Genetic Mapping Using Sequenced RAD Markers[END_REF] with the enzyme SbfI was done by MGX (CNRS, Montpellier). Each individual was identified using a unique 6 nucleotide tag and individuals were multiplexed in equimolar proportions by groups of 36 individuals. Each library was sequenced on one of 20 lanes of an Illumina HiSeq 2000 (libraries also included blue tit individuals from another site that was not analyzed in this study).

Raw sequences were inspected with FastQC (Andrews 2010) for quality controls.

Potential fragments of Illumina adapters were trimmed with Cutadapt [START_REF] Martin | Cutadapt removes adapter sequences from high-throughput sequencing reads[END_REF], 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 V1.39 [START_REF] Catchen | Stacks: an analysis tool set for population genomics[END_REF], , allowing for one mismatch in the barcode sequence. BWA-MEM 0.7.13 [START_REF] Li | Fast and accurate short read alignment with Burrows-Wheeler transform[END_REF] was used to map individual sequences against the reference genome of the great tit [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. Samtools 0.1.19 (Li et al. 2009) was used to build and sort bam files. Back in Stacks V1.39, we used pstacks to treat bam files, align the reads into matching stacks, infer loci, and detect SNPs at 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 catalogue of loci using n = 3 (number of mismatches allowed between sample loci when build the catalog). sstacks was used to match loci against the catalog. populations, the last Stacks program used here, genotyped individuals. Loci were retained if genotyped in at least 90% of individuals (all individuals from all sites grouped), with heterozygosity per site ≤ 0.60, and with individual minimal read depth of 10 ("na" replaced genotypes below a read depth of 10). VCFtools [START_REF] Danecek | The variant call format and VCFtools[END_REF] was used for further filtering of loci for a minimum average read depth of 20 across all genotypes, and a maximum average read depth of 100 across all genotypes. Individuals were genotyped for at least 90% of all loci; a filter would have otherwise been implemented to remove individuals with low genotyping rate. The dataset was filtered to retain loci with minimum allele frequency larger than 5% (MAF ≥ 0.05). Subsequently the dataset was pruned for linkage disequilibrium with the plink command LD indep 50 5 2.

Population genetic structure and genome wide relatedness matrix

Data were converted to several formats using the R package radiator [START_REF] Gosselin | Radiator: RADseq Data Exploration[END_REF].

Expected heterozygosity per site and Fst between sites were estimated using Genodive [START_REF] Meirmans | genotype and genodive: two programs for the analysis of genetic diversity of asexual organisms[END_REF]. A PCA resolving genetic distance between individuals was calculated using the R packages gdsfmt and SNPRelate (Zheng et al.

2012).

Genome wide relatedness matrices (GRM) were computed with GCTA (Yang et al. 2010;2011a) using the markers found on autosomes (i.e. excluding markers from sex chromosomes). We estimated GRMs for each of the three sites and for the three sites pooled. The GRM based on the three sites pooled was represented using a heatmap and histograms in order to depict relatedness diversity in the different sites (Figure 2). Finally, we computed several other GRMs for the unique purpose of chromosome partitioning (detailed in the corresponding section).

Social pedigree

Social pedigrees (PED-social) were constructed using the pedantics R-package [START_REF] Morrissey | pedantics: an r package for pedigree-based genetic simulation and pedigree manipulation, characterization and viewing[END_REF],on the basis of the complete pedigree generated from the whole long-term monitoring period. We first included all ringed individuals and assigned their mother and father based on observational data at capture. Unknown parents of a given nest were coded using a dummy identity in the PED-social to preserve sibship information. PED-social were pruned for each dataset studied (each population separately, and the combined pedigree of the three sites pooled) to retain only genotyped individuals and their ancestors (see supplementary material 2 for detailed characteristics of the pedigree).

Comparison between social pedigree and genome wide relatedness matrix and creation of a corrected social pedigree

We compared the PED-social and GRM for the 494 genotyped individuals mainly in order to infer the potential presence of erroneous pedigree relationship (particularly parent-offspring and sibling relationship). Indeed, errors may occur in social pedigrees, essentially due to extra-pair paternities, but also observational errors such as misreading of bird ring, successive egg-laying by different females in the same nest box, and brood parasitism. Such errors, notably extra-pair paternities, can downwardly bias estimates of heritability [START_REF] Charmantier | How do misassigned paternities affect the estimation of heritability in the wild?[END_REF]. We illustrated the concordances and discrepancies between PED-social and GRM using a biplot (Figure 2D). We then created corrected social pedigrees (PED-corrected) on the basis of discrepancies found between the PED-social and the GRM. For pairs of individuals having PED-social relatedness values of 0.5 (full-sibs or parent-offspring relationships) but GRM relatedness estimates below 0.1, we assumed that these relationships corresponded to erroneous parent-offspring relationship due to extra-pair paternity and assigned them a 0 in PED-corrected. For PED-social relatedness values of 0.5 but GRM estimates between 0.15 and 0.35, we assumed they corresponded to full-sibs being half-sibs and assigned them a value of 0.25 in PED-corrected. We also corrected false half-sibs in the pedigree (pairs of individuals having a 0.25 relatedness value in the PED-social but a GRM estimates smaller than 0.1) and assigned them a 0 relatedness values. We acknowledge that probably some more relationships could have been corrected based on discrepancies between PED-social and GRM, but the distribution of GRM values below 0.15 were largely overlapping and it would have been difficult to assign correct PED-corrected values. Using a method such as the one implemented in the R package sequoia [START_REF] Huisman | Pedigree reconstruction from SNP data: parentage assignment, sibship clustering and beyond[END_REF] would certainly have allowed the correction of these potential errors as well. Another potential limitation comes from the fact that ancestors of the genotyped individuals in the pedigree were not genotyped and therefore erroneous relationships between these individuals were impossible to correct.

Inference of heritability

We used both a frequentist and a Bayesian method to estimate heritability of the four traits. The frequentist method was implemented in GCTA. The Bayesian framework was implemented in the MCMCglmm R-package [START_REF] Hadfield | MCMC methods for multi-response generalized linear mixed models: the MCMCglmm R package[END_REF]. Bayesian inference is renowned to have a clear advantage over the other existing methods since the use of posterior distributions propagates the errors in estimates derived from animal models [START_REF] Morrissey | Bayesian approaches to the quantitative genetic analysis of natural populations[END_REF].

GCTA was used to estimate heritability for each trait, separately for each of the four GRMs. We here used best linear unbiased estimates (BLUPs) for each individual and phenotype since we could not implement more complex models in GCTA. BLUPs were estimated for the entire dataset of 494 genotyped individuals using a generalized linear mixed model with the MCMCglmm function, integrating a random effect for the observer and permanent environment effects accounting for multiple measurements of the same individual. Differences between sites and sex were also accounted for by adding these two factors as fixed effects. We then extracted the posterior mode of the BLUPS for each individual, and ran GCTA for each of the four GRMs and the four traits.

MCMCglmm was used to estimate heritability for each trait, separately for each of the four GRMs, the four PED-social and the four PED-corrected. Sex was integrated as a fixed effect in all models. Site was integrated as a fixed effect in the models that included all the individuals from the different sites. Measurer, identity and animal (corresponding either to the GRM, PED-social or PED-corrected matrix) effects were incorporated as random factors in order to partition the phenotypic variance into its observer, permanent environment and additive genetic components. Running models based on the GRMs, we used both BLUPs (to compare with GCTA) and alternatively all phenotypic measurements. For the models using GRMs and BLUPs, the only random factor was the animal effect (ie. the GRM) as permanent environmental and observer effects were already taken into account in the BLUPs. For the models based on PEDsocial and PED-corrected, we did not use BLUPs (since it could bias the analysis and/or give anticonservative results, [START_REF] Hadfield | The Misuse of BLUP in Ecology and Evolution[END_REF] but instead used all phenotypic measurements. We used identical parameters, priors and iterations for each estimate.

Random effects included additive genetic effects (V A ) estimated through the inclusion of pedigree data, permanent environmental effects (V PE ) accounting for repeated measurements of the same individual, measurer identity controlling for any potential confounding measurer effect (V OBS ) and residual variance (V R ). The models used can be described as follow:

Y = ! + !" + ! ! ! + ! !" !" + ! ! !"# + ! Equation 1
Equation 1 describes the animal models run on phenotypic traits with Y the vector of phenotypic observations for all individuals and the vector of mean phenotypes. Xb stands for the fixed effects (containing sex and site for the models on all individuals, and sex only for within site models). Z A , Z PE and Z O correspond to the random factors:

additive genetic (a), permanent environment (pe), and measurer (obs) random effects, respectively. e is the vector of residual errors. Posterior distributions were composed of 1,000 values per parameter. We used 120,000 iterations per model with sampling every 100 steps and with 20,000 discarded burn-in iterations. We used slightly informative priors to facilitate convergence, with V = V P / (r + 1), nu = 1, V P being the phenotypic variance, and r the number of random factors (results were quantitatively and qualitatively similar using uninformative priors, V = 1, nu = 0.002, but convergence takes a bit longer). We checked the models graphically. We verified that autocorrelations were less than 0.05. We finally reported posterior median, posterior mode, and 95% credible interval.

Inference of the effect of the number of SNPs on heritability estimates

We inspected whether the number of SNPs included resulted in i) an increase in heritability estimation precision (approximated via standard deviation, and LRT (Likelihood Ratio Test) in the case of GCTA) and ii) an increase in heritability (median) via saturation of the GRMs by SNPs in linkage disequilibrium with loci most likely causative of the phenotypic variation [START_REF] Gienapp | Genomic quantitative genetics to study evolution in the wild[END_REF]. We reported the median and standard deviation values of tarsus length heritability inferred using GCTA and MCMCglmm for several GRMs produced using a variable number of SNPs. Using GCTA, we analyzed 991 GRMs made up of a decreasing number of randomly chosen SNPs from the entire dataset, by step of 50 SNPs from the total number of SNPs. To confirm GCTA results, we analyzed a smaller number of 25 GRMs with MCMCglmm (much less than for GCTA since the Bayesian analysis was much more time consuming than the frequentist one) made from randomly chosen SNPs, concentrated mainly between 1000 to 10,000 SNPs (since GCTA indicated that the rate of improvement in estimates was particularly concentrated between these numbers).

Inference of genetic correlations

We used bivariate models in GCTA and in MCMCglmm in order to estimate genetic correlations between each of the four traits. These bivariate models were achieved using the same data, same fixed and random effects, and same number of iterations in the case of MCMCglmm, as for the univariate models. For MCMCglmm, we used slightly

informative priors with V = ! ! ! 0 0 ! ! !
and nu = 2 with Vp being the phenotypic variance, and r the number of random factors.

Chromosome partitioning

Finally, we partitioned heritability across the chromosomes using GCTA, using two methods: i) fitting the univariate GCTA model simultaneously on each GRM corresponding to each autosome; ii) fitting n times, where n is the number of autosomes, the univariate GCTA model simultaneously on two GRMs, one GRM computed for the focal chromosome and the other GRM computed using the other autosomes pooled. In both cases (i) and (ii), several microchromosomes (22, 27, 28, 25LG1, 25LG2, LGE22)

were pooled into one artificial autosome because they had too few SNPs. When a model did not converge, we discarded the smallest autosomes and ran the model again. We did so until the model would converge.

Results

Phenotypes

Among the 494 genotyped individuals, 494, 493, 470, and 488 were measured for tarsus length, wing length, bill length, and body mass, respectively (Table 1, Supplementary material 1). The sex ratio of genotyped individuals was near 0.5, with 246 females out of 494 individuals. Phenotypic traits varied slightly between sexes and between sites, justifying the inclusion of sex and site effects in the models partitioning phenotypic variance.

SNP calling and population genetic structure

The median number of reads per individuals was 6,066,514. The heterozygosity was 0.205, 0.205 and 0.204 in E-PIRIO, D-MURO and E-MURO, respectively (Supplementary material 3). Genetic differentiation between sites estimated using an Fst index was low yet significant, ranging from 0.006 to 0.008 (pvalue < 0.001; Supplementary material 3). The two first axes of the PCA (Figure 2A) explaining each 0.98% and 0.93% of the genotypic variance, and the heatmap of the GRM (Figure 2B), depicted the low genetic structure between the three sites.

Genome wide relatedness matrix and social pedigree

The GRM (Figure 2C) mostly included non-related (or distantly related) individuals.

Zooming in the histogram (Figure 2C) shows the presence of related individuals with parent-offspring or full sib-like links (grey), and half-sib-like links (blue). Comparing the relatedness values from the GRM to the one from the PED-social showed great consistency between the two matrices (Figure 2D). However, as expected (and shown elsewhere, Charmantier and Réale 2005), 112 links out of 858 known links presented PED-social relatedness higher than GRM relatedness by at least 0.10. Among these, 39 links exhibited PED-social relatedness of 0.5 but GRM relatedness of 0, and 36 links showed PED-social relatedness of 0.25 but GRM relatedness of 0. In addition, the GRM allowed the identification of many links that cannot be documented in the PED-social because of immigration, incomplete sampling or extra-pair paternity [START_REF] Pemberton | Wild pedigrees: the way forward[END_REF].

Notably, we found 460 links with PED-social relatedness was equal to 0 but GRM relatedness larger than 0.1. In total, we corrected 83 erroneous PED-social relatedness values (corresponding to ca. 10% of the total number of links known) thanks to the GRM analysis to create the PED-corrected.

Heritability estimates

In general, heritability measures based on PED-social, PED-corrected and GRM, using

BLUPs or all the measurements, and with MCMCglmm or GCTA, were highly consistent for the same trait (Table 2, Figure 3; see supplementary material 4 for variance estimates associated with the other random factors). The variations were very high among the three localities, for which sample sizes were very small. When pooling the three localities, variations between methods and credible intervals decreased.

Considering estimates from MCMCglmm and models fitting all phenotypic measurements, there was a consistent trend for lower heritability when using PED-social compared to GRM (Figure 3). For example, tarsus length and body mass heritability were lower by 22% and 31%, respectively, when using PED-social compared GRM.

Specifically, there was a trend for lower heritability when using PED-social than PEDcorrected. For example, tarsus length and body mass heritability were lower by 7% and 11%, respectively, when using PED-social compared to PED-corrected. There was also a trend for lower heritability when using PED-corrected compared to GRM. For example, tarsus length and body mass heritability were lower by 15% and 22%, respectively, when using PED-social compared to PED-corrected.

Heritability estimates were slightly higher when using BLUPs instead of all the phenotypic measures, together with the GRM, in MCMCglmm (e.g. only 3% higher for tarsus length, but 24% higher for body mass). Furthermore, credible intervals were slightly higher when using BLUPs. Using BLUPs and the GRM in GCTA, heritability estimates were relatively similar to the one obtained using MCMCglmm on all the phenotypic measures and either the GRM or the PED-corrected.

The number of SNPs included to compute GRMs influenced the tarsus length heritability estimates using either GCTA on BLUPs or MCMCglmm on all of the phenotypic measures (Figure 4). For this trait, heritability sharply increased from a handful of SNPs to approximately 5,000 SNPs. Heritability increased slowly from approximately 5,000 to 15,000 SNPs and then reached a plateau along which a very slow increase was however apparent. It is also worth mentioning a similarly shaped increase of likelihood ratio tests from GCTA.

Genetic correlations between traits

In general, genetic correlations between traits calculated with the different methods, were consistent for the same pair of trait (Table 3). However, the credible intervals were generally large, with 13 of the 30 intervals including 0. The largest estimates and smallest credible intervals were obtained for the genetic correlation between tarsus length and body mass, the most heritable traits. For each pair of traits there was a trend for a higher genetic correlation when using all measures than BLUPs.

Chromosome partitioning

When fitting the univariate GCTA model simultaneously on each autosome, the correlation between heritability and chromosome length was positive for tarsus length (r2 = 0.26; p-value = 0.009), being the most heritable traits of this study. The correlation was non-significant for body mass (r2 = 0.08; p-value = 0.20), wing length (r2 = 0.03; p-value = 0.41), and bill length (r2 = 0.00; p-value = 0.85). Standard error intervals for chromosome scaled heritability did not include zero only for chromosomes 1A and 3 for tarsus length (Table 4, Figure 5), and for chromosome 6 for body mass. In the case of tarsus length, no autosome had to be removed to enable the convergence of the model fitting all the autosomes simultaneously. In contrast, for wing length, bill length and body mass, the microchromosomes group had to be removed and respectively one, 12 and three additional autosomes (by order of increasing size) had to be removed to enable convergence of the models. When fitting the univariate GCTA model separately for each autosome versus the rest of the genome, several heritabilities were likely overestimated since the total heritability was larger than with the previous method (h2 = 1.47 vs 0.84 for tarsus length, 0.85 vs 0.36 for wing length, 1.11 vs 0.31 for bill length, and 0.95 vs 0.62 for body mass).

Discussion

Genomic analyses offer great opportunities for measuring relatedness precisely among individuals from natural populations and hence estimating heritability and genetic correlations among traits in the wild [START_REF] Edwards | Next-generation QTL mapping: crowdsourcing SNPs, without pedigrees[END_REF][START_REF] Gienapp | Genomic quantitative genetics to study evolution in the wild[END_REF]. Heritability partitioning along the genome can also help in identifying regions explaining the additive genetic variation for polygenic traits (Yang et al. 2011b). Determining trait heritability and its partitioning along the genome, and understanding trait genetic covariances, are also essential first steps to detect individual loci in the genome that contribute to trait differences between individuals, using for example genome wide association or linkage mapping [START_REF] Schielzeth | Challenges and prospects in genome-wide quantitative trait loci mapping of standing genetic variation in natural populations[END_REF]). Although such approaches have been recently

tested in a few datasets of well-studied wild organisms [START_REF] Robinson | Partitioning of genetic variation across the genome using multimarker methods in a wild bird population[END_REF][START_REF] Stanton-Geddes | Estimating heritability using genomic data (J Hadfield[END_REF][START_REF] Bérénos | Estimating quantitative genetic parameters in wild populations: a comparison of pedigree and genomic approaches[END_REF], we require more explorations on whether smaller datasets could be used for such procedures and on how they compare with the classic animal model approach based on long-term pedigrees. In this context, we used both social pedigrees (PED-social) and genome wide relatedness matrices (GRM), to compute heritability estimates, genetic correlations, and chromosome partitioning for several phenotypic quantitative traits (tarsus length, wing length, bill length, and body mass) in a rather small dataset of individuals sampled in wild populations of Corsican blue tits and genotyped at nearly 50,000 RADseq derived SNPs. We discuss the power of such relatively small GRM, which may be typically produced while studying wild populations of non-model species, to examine the aforementioned quantitative genetic parameters.

Our main findings were i) a high congruence between heritability estimated using GRM, PED-corrected and PED-social. However, the median heritability obtained using the GRM was generally slightly higher than using PED-corrected and PED-social, which may be attributed to absence of erroneous links, greater precision and higher information content of the GRM. ii) The number of markers was not an issue for computing the GRM and estimating heritability, when approximately 15,000 SNPs were reached. iii) Genetic correlations among traits and chromosome partitioning were less likely than heritability estimates to be very informative or robust given the large credible intervals and the number of non-significant estimates, most likely due to power issues (insufficient number of SNPs per chromosome in the case of the micro-chromosomes and insufficient number of individuals). We thereafter discuss the usefulness and limitations of inferences of quantitative genomic parameters, including heritability, in relatively small sample sizes, with a medium number of genetic markers, for wild populations of non-model species.

Comparing PED-social and GRM

The GRM obtained differed in two aspects from the PED-social (Fig 2D ): higher density and relatedness variance on one hand, and no pedigree errors on the other hand. The first major difference was that the GRM enabled estimating relatedness among individuals that were not connected in PED-social (Figure 2D; 460 GRM relatedness values were higher than PED-social relatedness values by at least 0.1). This concerned individuals at every degree of relatedness but obviously primarily among individuals from disconnected families. Therefore, GRM increased the depth and also the size of the pedigree. This illustrates the usefulness of GRM in cases for which a dense and deep pedigree is difficult to obtain due to for example dispersal, large population size compared to the fieldwork capacity, and family structure [START_REF] Pemberton | Wild pedigrees: the way forward[END_REF].

Furthermore, GRM allows capturing much more variance in realized relatedness among individuals, both close and more distantly related (Figure 2D), which may also increase the power of such GRM-based estimates in large and open populations.

The second major difference between PED-social and GRM was that several individuals had much higher relatedness in PED-social than in the GRM, most likely indicating errors in PED-social. The most common cases consisted of wrongly assigned parentages (relatedness of 0.5 in PED-social and 0 in GRM), half sibs most likely originating from extra-pair paternities (relatedness of 0.5 in PED-social and 0.25 in GRM), and the spread of such errors in the pedigree (other smaller relatedness in GRM than in PED-social). Using this comparison between GRM and PED-social, we created a PED-corrected in which these aforementioned erroneous links, accounting for approximately 10% of the PED-social links, were corrected. This proportion of erroneous linked corrected in our dataset should be expected since 18 to 25% of young have been previously identified as extra-pair offspring in these Corsican blue tit populations (Charmantier et al. 2004a).

Heritability estimates from PED-social, PED-corrected and GRM

Heritability estimated for tarsus length (0.62 to 0.81 depending on the method, for three sites pooled, see Table 2) and body mass (0.34-0.61) were on the upper limit compared to previous estimates in this species and for other close species [START_REF] Jensen | Evolutionary dynamics of a sexual ornament in the house sparrow (Passer domesticus): the role of indirect selection within and between sexes[END_REF][START_REF] Husby | Testing Mechanisms of Bergmann's Rule: Phenotypic Decline but No Genetic Change in Body Size in Three Passerine Bird Populations[END_REF][START_REF] Postma | Four decades of estimating heritabilities in wild vertebrate populations: improved methods, more data, better estimates?[END_REF][START_REF] Delahaie | Conserved G-matrices of morphological and life-history traits among continental and island blue tit populations[END_REF]) while heritability for wing length (0.26-0.34) and bill length (0.14-0.21) were rather on the lower ranges usually found.

Here again it should be noted that we did not optimize the animal models by integrating factors that may have contributed to explain trait variation, for example to account for daily and seasonal variation in body mass. Moreover, it has been observed that the number of individuals included can affect the heritability estimates [START_REF] Silva | Insights into the genetic architecture of morphological traits in two passerine bird species[END_REF].

Credible intervals and standard errors for heritability estimates were much larger than in previous studies based on larger sample size [START_REF] Robinson | Partitioning of genetic variation across the genome using multimarker methods in a wild bird population[END_REF][START_REF] Stanton-Geddes | Estimating heritability using genomic data (J Hadfield[END_REF][START_REF] Bérénos | Estimating quantitative genetic parameters in wild populations: a comparison of pedigree and genomic approaches[END_REF][START_REF] Silva | Insights into the genetic architecture of morphological traits in two passerine bird species[END_REF]. The credible intervals per site were larger than for the three sites pooled, suggesting the decrease of power with decreasing sample size. Similarly, [START_REF] Silva | Insights into the genetic architecture of morphological traits in two passerine bird species[END_REF] obtained increased standard errors when decreasing sample size. This most likely suggests that the number of individuals used here conferred limited power.

Overall, when using MCMCglmm and models using all of the phenotypic measurements, heritability estimates were relatively similar for GRM, PED-corrected and PED-social for the same trait (Figure 3). These results are in line with previous comparisons realized with much larger sample size [START_REF] Robinson | Partitioning of genetic variation across the genome using multimarker methods in a wild bird population[END_REF][START_REF] Bérénos | Estimating quantitative genetic parameters in wild populations: a comparison of pedigree and genomic approaches[END_REF]). However, unlike these previous studies showing very small differences in heritability based on PED-social and GRM, our estimates based on GRM were larger than using PED-social and to a lesser extent, than using PED-corrected, for tarsus length and for body mass, but not for wing length and bill length. Here, wing length heritability was 0.28 using PED-social, 0.26 using PED-corrected and 0.26 using GRM (Table 2). In turn, tarsus length heritability obtained from GRM was 18% higher than from PED-corrected and the one from PED-corrected 8% higher than from PED-social.

Similarly, body mass heritability obtained from GRM was 29% higher than from PEDcorrected and the one from PED-corrected 12% higher than from PED-social. The difference between PED-corrected and PED-social for these two traits may be attributable to the erroneous links in PED-social. Indeed, it has been estimated that PED-social containing 5 to 20% of extra-pair paternities result in underestimated heritability by up to 17%, in blue tits [START_REF] Charmantier | How do misassigned paternities affect the estimation of heritability in the wild?[END_REF]. Then, the differences observed between heritability estimated using GRM and PED-corrected are likely to originate from higher density of the GRM and the fact that GRM incorporate more variance in relatedness for a given class of kinship. Further changing the pedigree by not only correcting erroneous existing links but also informing previously unidentified links (eg. linking previously unknown extra-pair fathers to their offspring) based on the genetic data, could also decrease the heritability difference between PED-corrected and GRM.

Heritability estimates, and their credible intervals, based on BLUPs were slightly higher compared to estimates based on models using all phenotypic measurements when using MCMCglmm (this was not true for GCTA). This is not surprising given that BLUPs are known to artificially increase the precision of the measure and therefore increase statistical significance [START_REF] Hadfield | The Misuse of BLUP in Ecology and Evolution[END_REF]. While the difference was low for tarsus and bill length (3% and 5% higher, respectively), the difference was particularly pronounced for wing length and body mass (31% and 24%, respectively). One speculation is that wing length and body mass are less repeatable than tarsus length and bill length and that our BLUPs did not adequately take into account individual variance.

Effect of the number of SNPs on the heritability estimates based on GRM

The number of SNPs had a relatively minor effect on heritability estimates, from approximately 15,000 SNPs (Figure 4). This is concordant with what has been shown in previous studies [START_REF] Stanton-Geddes | Estimating heritability using genomic data (J Hadfield[END_REF][START_REF] Bérénos | Estimating quantitative genetic parameters in wild populations: a comparison of pedigree and genomic approaches[END_REF]. Of course, the effect of the number of SNPs on the heritability estimation depends on several parameters (e.g. genome size, LD across the genome, LD among markers, genetic architecture of the focal trait, polymorphism of the markers, genotyping quality, missing data) and is not transposable from one study to another. These days, the number of markers may no longer be a problem, given the rise of the accessibility of many genomic tools. However, the genotyping quality and missing data occurrence may be challenging for studies using RADseq and GBS methods, as noted by [START_REF] Gienapp | Genomic quantitative genetics to study evolution in the wild[END_REF]). Here, we prioritized the quantity and the quality of the RAD markers over the number of individuals, as well as genotyping only individuals with high quality DNA, ending up with a relatively high read depth and low missing rate. Future studies aiming to use RADseq derived markers to estimate heritability should also find a tradeoff between the number of genomic markers and the number of genotyped individuals. A formal simulation study determining such a tradeoff between numbers of individuals and of markers and markers read-depth (directly linked to quality) would be highly valuable for guiding the design of RADseq analyses in the context of quantitative genomics.

Genetic correlations using PED-social, PED-corrected and GRM

Genetic correlations obtained using GRM, PED-social and PED-corrected had large credible intervals, preventing us to make robust interpretations, except that our sample size was probably too small to obtain acceptable credible intervals. This was partly expected since bivariate models are known to be data hungry. While such lack of power is a recurrent issue when estimating genetic correlations even with much larger sample sizes [START_REF] Vattikuti | Heritability and Genetic Correlations Explained by Common SNPs for Metabolic Syndrome Traits (PM Visscher[END_REF][START_REF] Visscher | Statistical Power to Detect Genetic (Co)Variance of Complex Traits Using SNP Data in Unrelated Samples (GS Barsh[END_REF][START_REF] Ni | Estimation of genetic correlation using linkage disequilibrium score regression and genomic restricted maximum likelihood[END_REF], it was nevertheless a disappointment while comparing these results to the one from [START_REF] Bérénos | Estimating quantitative genetic parameters in wild populations: a comparison of pedigree and genomic approaches[END_REF], in which the credible intervals were relatively small. The median values of genetic correlations for these traits were in line with previous studies in other bird species (Teplitsky et al. 2014).

Chromosome partitioning based on GRM

Chromosome partitioning of heritability of tarsus length was on par with what has been observed in previous studies, i.e. increasing heritability of quantitative traits with increasing number of genes per chromosome (Yang et al. 2011b;[START_REF] Robinson | Partitioning of genetic variation across the genome using multimarker methods in a wild bird population[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][START_REF] Santure | Replicated analysis of the genetic architecture of quantitative traits in two wild great tit populations[END_REF][START_REF] Bérénos | Heterogeneity of genetic architecture of body size traits in a free-living population[END_REF][START_REF] Wenzel | Genome-wide association and genome partitioning reveal novel genomic regions underlying variation in gastrointestinal nematode burden in a wild bird[END_REF][START_REF] Silva | Insights into the genetic architecture of morphological traits in two passerine bird species[END_REF].

Similarly to several of these previous studies, there were few chromosomes for which standard error did not include zero. Of particular interest, chromosomes 3 and 1A

showed high heritability for tarsus length. These chromosomes may be of particular interest for future studies investigating the genomic bases of this quantitative trait variation in blue tits. Chromosome 3 was already identified in a recent study on morphology of Dutch and British great tits [START_REF] Santure | Replicated analysis of the genetic architecture of quantitative traits in two wild great tit populations[END_REF] as explaining high heritability for tarsus length. In contrast, chromosome 1A explained a low proportion of genetic variance in this last study. These kind of discrepancies might be congruent with the fact that this same study [START_REF] Santure | Replicated analysis of the genetic architecture of quantitative traits in two wild great tit populations[END_REF]) also showed very little consistency between the variance explained by each chromosome for great tit populations from the Netherlands and from the United-Kingdom.

Regarding the three other traits analyzed, the chromosome partitioning results appeared much less robust since several chromosomes had to be removed (7, 18, and 29 9 autosomes removed, for wing length, bill length and body mass, respectively) to enable models' convergence. The need for such chromosomes removals has previously been reported [START_REF] Wenzel | Genome-wide association and genome partitioning reveal novel genomic regions underlying variation in gastrointestinal nematode burden in a wild bird[END_REF][START_REF] Silva | Insights into the genetic architecture of morphological traits in two passerine bird species[END_REF]) and has received detailed consideration by [START_REF] Kemppainen | Inference of genetic architecture from chromosome partitioning analyses is sensitive to genome variation, sample size, heritability and effect size distribution[END_REF]. In addition, the relationship between the number of genes and the heritability explained by chromosomes were deviating from theoretical expectations of increasing heritability with increasing chromosome size for these polygenic traits. For comparison, [START_REF] Bosse | Recent natural selection causes adaptive evolution of an avian polygenic trait[END_REF] showed a significant increase of the proportion of additive genetic variance of bill length explained by a chromosome in relation to chromosome size in great tits. Such deviation may not have been affected by removing chromosomes [START_REF] Kemppainen | Inference of genetic architecture from chromosome partitioning analyses is sensitive to genome variation, sample size, heritability and effect size distribution[END_REF]) but most likely originated from a limited power of the small sample size, alongside relatively small heritabilities [START_REF] Kemppainen | Inference of genetic architecture from chromosome partitioning analyses is sensitive to genome variation, sample size, heritability and effect size distribution[END_REF].

Conclusion on the usefulness of GRM based heritability measures using small samples from wild populations

Overall, our study reveals that RADseq data on around 50k SNPs (as advised by Berenos et al 2014) for around 500 phenotyped individuals can provide estimates of heritability that are close to, and probably more accurate than, estimates based on social pedigree data from seven years monitoring for more than 1600 individuals (as included in the PED-social used here). This opens very interesting avenues in the field of quantitative genetics since estimating heritability and genetic correlations in the wild have long been restricted to study systems where long-term monitoring is feasible.

RADseq data could allow the estimation of quantitative genetic parameters for a much lower cost in terms of time (years) spent in the field, although individual capture and phenotyping will most likely remain time consuming. Using GRM on few individuals gives the power to estimate such parameters on virtually any species even though no longterm pedigree is available (Gienapp 2017). This will notably allow moving on from individual study to more comparative approaches and answer important questions such as understanding the spatial variation of evolutionary potential or the role of evolutionary constraints in phenotypic evolution. Such a genomic approach also provides the possibility to explore genetic covariance between traits as well as chromosome partitioning to confirm the polygenic architecture for phenotypic traits classically measured in birds. Less than 500 individuals appears however to provide insufficient power to correctly estimate either genetic covariances between traits or contribution of individual chromosomes to overall heritability. Table 1. Geographic coordinates, number of genotyped individuals (N), and trait mean values (±SE) for this sample, for 868 each of the three sites and for the entire dataset. The phenotypic variance (Vp) is also given for each trait.

Table 3. Genetic correlations (rG), and confidence (GCTA) / credible (MCMCglmm) intervals (95%CI), between traits estimated for the individuals from the three site pooled, using GRM, the PED-social or the PED-corrected, based on either BLUPS or all the measures, and implementing models in GCTA or MCMCglmm.

  g. 200 barrelclover (Medicago trunculata) individuals genotyped at more than 5 million SNPs. Wenzel et al. (2015) estimated genome-wide heritability in 695 red grouse (Lagopus lagopus scotica) genotyped at 384 SNPs. Silva et al. (2017) estimated quantitative genomic parameters in 1,898 house sparrows (Passer domesticus) genotyped at 6,348 SNPs and in 825 collared flycatchers (Ficedula albicolis) genotyped at 38,689 SNPs.

  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).

Figure captions Figure 1

 captions1 Figure captions

Figure 2 .

 2 Figure 2. Genetic structure and relatedness among individuals. A-Biplot of the two first

Figure 3 .

 3 Figure3. A-Heritability estimates using the GRM, the PED-social or the PED-corrected,

Figure 4 .

 4 Figure 4. Effect of the number of SNPs used to build the GRM on the tarsus length

Figure 5 .

 5 Figure 5. Chromosome partitioning. The univariate GCTA model was fitted

  Figure 4.

  The median read depth per individual and per locus was 56. The Stacks program population outputted 52,783 loci totaling 96,009 SNPs. After MAF pruning, we retained 41,986 loci totaling 68,114 SNPs. After LD pruning, we retained 38,030 loci totaling 49,682 SNPs. 47,865 of these SNPs were on autosomes. The number of filtered SNPs per chromosome ranged from 5523

	(chromosome 2) to 56 (LGE22).

Table 4 .

 4 Chromosome partitioning of the heritability estimated for the four phenotypic traits using GCTA. Partitioning was

			Tarsus	Tarsus	Tarsus	Wing	Wing	Bill
			length	length	length	length	length	length
			Wing	Bill	Body	Bill	Body	Body
			length	length	mass	length	mass	mass
	GCTA	rG	0.40	0.37	0.48	-0.06	0.22	0.13
	BLUPs	0.05	0.05	-0.14	0.23	-0.86	-0.022	-0.49
	GRM	0.95	0.75	0.88	0.73	0.86	0.67	0.76
	MCMCglmm	rG	0.41	0.32	0.52	0.10	0.21	0.03
	BLUPs	0.05	-0.08	-0.09	0.20	-0.57	-0.28	-0.50
	GRM	0.95	0.67	0.73	0.75	0.48	0.66	0.59
	MCMCglmm	rG	0.48	0.45	0.56	0.50	0.50	0.33
	all measures	0.05	0.16	0.06	0.37	-0.08	0.05	0.05
	GRM	0.95	0.81	0.80	0.76	0.87	0.71	0.88
	MCMCglmm	rG	0.48	0.43	0.57	0.46	0.49	0.52
	all measures	0.05	0.13	0.00	0.40	-0.12	0.06	0.05
	PED-corrected	0.95	0.79	0.81	0.78	0.89	0.74	0.84
	MCMCglmm	rG	0.46	0.52	0.61	0.42	0.41	0.37
	all measures	0.05	0.17	0.16	0.38	-0.13	0.09	0.08
	PED-social	0.95	0.77	0.79	0.82	0.79	0.72	0.81

B 48 Figure 3.

Acknowledgments

We thank all the people who helped to maintain the study sites and to conduct the blue tit monitoring, in particular Christophe de Franceschi, Claire Doutrelant, Arnaud Grégoire, Marcel Lambrechts, Jacques Blondel and Philippe Perret. We thank Pierre Alexandre Gagnaire for useful advices on genomics. We thank three anonymous reviewers and Anna Santure for precious comments on earlier versions of this manuscript. This project was funded by the European Research Council (Starting grant ERC-2013-StG-337365-SHE to AC) and a long-term funding from OSU-OREME.

Data archiving

Data (Plink files, PED-social, and repeated and BLUPs phenotypic traits values) are available from the Dryad Digital Repository: doi:10.5061/dryad.k6r1mk8

Conflict of interest

The authors declare no conflict of interest.

Author contributions

Conception

Supplementary material

Supplementary material 1. Box plots showing sex-and site-specific median and variation in tarsus length, wing length, bill length and body mass for the 494 genotyped individuals. The first four letters code for the site and the last letter codes for the sex (F: female; M: male). Supplementary material 2. Social pedigree characteristics. Supplementary material 3. Heterozygosity (He) in E-Piro, D-Muro and E-Muro and FST among the three sites (FST above diagonal and p-value below the diagonal). Supplementary material 4. Variance estimates associated to the other random factors.

Supplementary material 1.


Adeline Manuel
email: adeline.manuel@map.cnrs.fr

Anas Alaoui
email: anas.alaoui@map.cnrs.fr

M ' Darhri

Violette Abergel
email: violette.abergel@lyon.archi.fr

Fabien Rozar
email: fabien.rozar@gmail.com

Livio De Luca
email: livio.deluca@map.cnrs.fr

A semi-automatic 2D/3D annotation framework for the geometric analysis of heritage artefacts

Keywords: CH documentation, semantic annotation, geometrical and visual descriptors, photogrammetry, segmentation

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

INTRODUCTION

Heritage asset monitoring and management are multidisciplinary fields by nature. They involve a plurality of actors, methods, and tools related to a wide variety of scientific issues and objectives. This context naturally leads to a phenomenon of dispersion of data, information and knowledge, which poses major problems in terms of sustainability. The heterogeneity of this context also causes difficulties in linking different levels of observation and analysis of the built environment by the various actors involved. Finally, the sensitive nature of the studied sites and their fragility, as well as their evolution over time, imply the need for regular monitoring that requires a constant supply of information.

Gathering of analysis strategies

Based on these findings, our main objective is to merge all the different aspects of the patrimonial survey (geometry, visual appearance, and semantics) within an integrated documentary approach built on an informative continuity that merges the acquisition, analysis and interpretation phases in a collaborative framework. This approach raises three main problematics :

• Semantic: To bring concrete meaning to geometric data

• Collaborative : To allow the various actors of the heritage to better communicate their results

• Interoperability: To provide a common platform, independent of the device and the operating system, so that collaborators can exchange their analyses in a transparent way.

I. STATE OF THE ART

Heritage assets study involves the collection of numerous data in addition to visual representations. This data must be structured, based on theoretical work on the naming of concepts and their interconnection [START_REF] Doerr | Ontologies for cultural heritage[END_REF]. The general goal is to associate semantic data with documents, which can take the form of simple vocabulary terms, or through ontologies, structuring terms and concepts related to a specific field of knowledge ( [START_REF] Messaoudi | An ontological model for the reality-based 3d annotation of heritage building conservation state[END_REF] [3] [START_REF] Wefers | Development of a platform recommending 3d and spectral digitisation strategies[END_REF] ). Several studies have been conducted on this association, and more particularly on the semantic enrichment of 3D digitization ( [START_REF] Havemann | The Arrigo Showcase Reloaded -towards a sustainable link between 3D and semantics[END_REF] [6] [7] [8] ). [START_REF] Shi | CHER-Ob: A Tool for Shared Analysis in Cultural Heritage[END_REF] also offers the ability to integrate and annotate different types of 2D media in addition to the digital model. In terms of web visualization, the new web technologies like the WebGL API have enabled the development of many viewers dedicated to the manipulation of 3D data ( [10] [11] [START_REF] Schuetz | Potree: Rendering large point clouds in web browser[END_REF] ), possibly annotated. Nevertheless, most of them consider 2D and 3D visualization as two distinct subjects. Let us however note [START_REF] Graf | smapshot -geolocate our images on a virtual 3d globe![END_REF], whose objective is to allow users to geo-reference archive photographs, then visualize them by transparency on a 3D representation of the globe, or [START_REF] Brivio | Photocloud: Interactive remote exploration of joint 2d and 3d datasets[END_REF] and [START_REF] Snavely | Photo tourism: Exploring photo collections in 3d[END_REF] which also make it possible to visualize 3D and oriented images within a single scene. However, if these viewers allow to spatialize 2D resources, they do not offer specific interaction modalities or annotation visualization. Finally, [START_REF] De Luca | Semantic 3d media and content: A semantic-based platform for the digital analysis of architectural heritage[END_REF] offers all of these features, but uses meshes instead of raw data, on which annotations are overlaid as textures [START_REF] Stefani | Developing a toolkit for mapping and displaying stone alteration on a web-based documentation platform[END_REF]. This data structuring does not allow a direct link to be established between the real object and its digital model and is therefore not sufficient to construct the expected digital epidermis around the studied object.

II. MAIN APPROACH

Our approach -the Aïoli platform [START_REF] Aioli | A reality-based 3d annotation cloud platform for the collaborative documentation of cultural heritage artefacts[END_REF], currently under development -is based on the observation that the element most likely to federate all the above-mentionned actors is the "physical" entity of the heritage object under study. The ambition of this platform is to establish a link between the real object and its "digital ecosystem", in order to guarantee a continuum of information ranging from data acquisition to the construction of semantically enriched 3D digital representations. This information continuum is established, from photographs of the studied object, thanks to a processing pipeline integrated into the Aïoli platform. The latter responds to contextual issues by considering the study of an object (artifact or heritage building) as a "project", and by providing users with a working environment -a workspace -integrating all the above approaches.

A. Workspace Preparation

The first step is based on photogrammetry, which allows us, from an initial set of images, to construct a point cloud in three dimensions, a representation of the object acting as a digital clone [START_REF] Pierrot-Deseilligny | Automated image-based procedures for accurate artifacts 3d modeling and orthoimage generation[END_REF]. The second phase of this implementation is the point cloud indexing, in order to establish a correspondence between the points of the cloud and the pixels of the images [START_REF] Manuel | 2D/3D Semantic Annotation of Spatialized Images for the Documentation and Analysis of Cultural Heritage[END_REF]. The bijective relationship thus established allows us to know, for each 3D point, the list of images on which the point appears, and for each of them, the corresponding pixels. Conversely, starting from one or more point(s) in an image, we can get a list of the corresponding 3D points (Figure 1).

B. Spatializing semantic data

This principle is more specifically used to allow users to draw annotations (observation areas) on an image and automatically reproduce them not only on the point cloud but also on all other images on which this area is visible. By allowing the various actors to associate to these annotations semantic descriptions coherent with their field of study, we finally spatialize information, by storing it as close as possible to the real object.

C. Annotation Containers

Once the working environment is ready, the users only have to define "annotation layers", named data structures. These structures allow to group in a semantic way the annotations which would have common aspects. These layers give us a way to gather similar regions, and offer them to the users' browsers for display, with a familliar way of storing data, since the layers mechanism is already implemented in a lot of different softwares. They also are a good way to organize multi-user works, as they're holding different rights (read-only, read/write) for each user. In this way, the project owner can control more precisely who can see his work, and who can edit and enhance his analysis. Finally, they also are a way to hold formal description sheets, which normalize the annotations for each layer, as described later.

Fig. 2. Example of annotated layers on the same image [START_REF] Messaoudi | An ontological model for the reality-based 3d annotation of heritage building conservation state[END_REF] For example, for the monitoring of architectural heritage [START_REF] Messaoudi | An ontological model for the reality-based 3d annotation of heritage building conservation state[END_REF], the user can create 3 layers, one layer named "Degradations", supporting all the annotations describing degradations of the material (cracks, detachment of matter, biological colonization...), one named "Material", created to describe the materials used for such constructions, and finally one named "Architectural elements", which will hold all the regions describing classical architectural elements (walls, columns, vaults...) (Figure 2). The comparison of these different layers then permits to identify the materials most exposed to particular alterations, or to target -with more accuracy -the architectural elements to be restored.

D. Formal Description sheets

The user-defined descriptors add a way to store information within an annotation, allowing the layer owner to define flexibly their own data structure, by adding fields: text, date, URL, concept from a thesaurus, and even attached documents (images, pdfs, sound, video). In this way, the user makes sure that their working tool will be adapted to their observations and to the needs of their study. Thus, all the annotations of a layer inherit the parent structure, which ensures a certain semantic uniformity. In the "Degradations" layer of the example above, 3 description fields could be defined: The type of degradation, the date of the survey and different observations related to the degradations. When a user characterizes a degradation, he must then fill these fields, fixed for each annotation of the layer. This structure allows the user to visualize, collect and organize data in a semantic way. These are mainly qualitative data, related to the users' fields of study; but what about quantitative data describing the object's geometry, taking advantage of the three-dimensional resources?

III. AUTOMATIC DATA ENHANCEMENT: GEOMETRICAL & VISUAL SUPPLY
The characterization of a point cloud involves extracting data representative of all the points or the geometry at each point of the cloud.

A. Computed Quantitative Descriptors

For each point cloud, whether it is the complete point cloud of the object, or the 3D representation of an annotation, a set of "basic" descriptors can be extracted to bring out information related to the region involved. We extract numerous data such as the number of points composing the 3D region, the barycenter, information related to the bounding box (position, center, magnitude, volume), the average orientation (average normal) and the average color automatically for each point cloud representation. In addition, for normal vectors and point colors, a range of representative values for the region is extracted to nuance the information provided by the simple average. Thus, to describe each annotation, the user has, on the one hand, semantic data associated by himself through the description sheets and on the other hand, numerous data automatically generated by the system. All this data can then be used to define queries, for research, classification, and comparaison purposes.

B. Visual Geometrical Descriptors

Images of the object in the visible spectrum do not always allow a good appreciation of the different characteristics of the object, such as shapes or surface appearance, and can therefore sometimes make it difficult to select a region. Thus, to facilitate analysis and interpretation operations, we propose the generation of images highlighting the geometric information extracted automatically on the point cloud. The estimation of the values of curvature, roughness and ambient occlusion [START_REF] Duguet | A Point-Based Approach for Capture, Display and Illustration of Very Complex Archeological Artefacts[END_REF] for each point of the cloud are made comprehensible on a point cloud by the use of a color gradient. Thanks to the 2D/3D relationship, this information is reprojected in the view of each of the images and thus generate images color-coding this geometric information of the point cloud (Figure 3). These images can then be used as support for the selection in the same way as photographies. This supposes therefore to set up a dynamic viewer allowing at the same time to display the point cloud, the source images, but also the reprojection of the various "textures" of the images (normal, curvature, roughness, ambient occlusion, etc...).

C. Interaction with the information continuum

The convergence of these analysis tools resulting from this new informational continuum implies the definition of a new approach concerning user interactions. More specifically, this led us to develop a hybrid 2D/3D viewer allowing users to take advantage of the bijective relationship by offering interaction tools that can be adapted to different types of device. For this, the first phase was to literally translate within a dynamic viewer the relationship linking 2D and 3D, positioning and orienting the different iconographic resources relative to the point cloud, in order to create a coherent representative space (Figure 4). In this way, the user can position himself on each view and visualize in transparency the underlying 3D data. For this purpose, he benefits from all the available visual supports: both the original images and the images created by projection of the calculated geometric descriptors.

A second phase was to propose interaction tools adapted to the practice of annotation on 3D models. The 2D annotation tools were therefore completed by a real-time projection of the plot control points on the point cloud, in order to offer a 3D preview of the current annotation. As long as it is not validated, the trace remains editable. It is now also possible to define an annotation directly on the point cloud, and to edit it, before propagating it to all the project resources. Finally, the multiplicity of profiles of actors that can intervene on the same project, and the complexity of some studied entities led us to propose more varied representation tools, allowing users to define their own graphic charter. Thus, users can now create annotation typologies, for which they can define the type of line (dotted or solid line) and thickness, spacing, color, opacity, filling, etc. Thus, during an annotation, the user selects in his personal list the relevant type of trace in order to obtain a coherence in the graphic representation. Such a space of visualization and interactions presents then a double interest: it allows not only to facilitate the study of a heritage object by facilitating the observation and the collaboration, but it also increases considerably the traceability and lasting quality, thanks to the safeguard of the various handling leading to the elaboration of a data. In fact, through this approach we aim to assign as much value to processes as to their results.

IV. POINT CLOUD SEGMENTATION

Segmentation methods offer automatic solutions for subdividing a digital model. The principle is based on the clustering of points in sets that share a common property (geometric or colorimetric). Many segmentation methods exist [START_REF] Grilli | A review of point clouds segmentation and classification algorithms[END_REF] and all face difficulties related to the precision of the point cloud generated (number of points, irregularity, noise, outlier points...). In order to offer users a simple way to create annotations, the use of segmentations is currently being implemented on our platform. The idea is to automatically detect on the point cloud sets of points sharing a common characteristic and then provide the user with regions to annotate predefined by segmentation. Thus, for the time being, only a plane segmentation tool has been integrated in order to validate our approach.

A. Plane segmentation

A plane is the simplest surface to identify from a mathematical perspective in regards to its geometric characteristics. In order to overcome the difficulties aforementioned, our planesegmentation is based on an octree-based region growing method initially created in [START_REF] Wang | Lidar data segmentation and classification based on octree structure[END_REF]. This method, described in [START_REF] Vo | Octreebased region growing for point cloud segmentation[END_REF] works in two iterative steps:

• Split : This step divides the input point cloud into sets of points in the form of an octree i.e. in the form of cubes containing points satisfying a same refinement criterion.

In our case, this criterion is satisfied for a cube if the average distance of the points from their best fit plane is lower than the user-defined threshold. Starting initially from a cube containing all the points, each iteration subdivides the different cubes by testing, for each of them, the refinement criterion until all the cubes satisfy it. This process generally builds an irregular octree where each leaf contains the points closest to their best fit plane.

• Merge: This step groups the points identified as coplanar into regions. A region is built with the following process: it is initialized with an octree leaf and then iteratively constructed by exploring and adding neighbouring octree leaves as long as the merge criterion is met. In our case, this criterion is satisfied if the angle between the normal vector of the current region and the one of the considered octree leaf is lower than a user-defined threshold. The set of different regions is therefore obtained when all leaves of the octree have been tested and assigned.

B. Integration with Aïoli platform

The result of this segmentation gives a clustering of the 3D points, forming different regions which can then be reproduced on each image thanks to the 2D/3D indexation mechanism (Figure 5). This feature allows automatic generation of an annotation layer grouping the regions detected by segmentation. The user only have to associate the semantic descriptions to each of the newly created layer regions. The plane segmentation algorithm above-mentionned is one of many examples we could use. The Aïoli system is built on Docker© [START_REF] Docker | Docker : Build, ship and run apps anywhere[END_REF] to encapsulate the different software used to process the projects. We benefit from this flexible architecture by using all external dedicated software (photogrammetry engine MicMac, point cloud indexation, annotation propagation, and also segmentation plugins) as Docker images, so we can use all software components the same way. This allows us to not only run efficiently by preserving computing resources, but also be able to easily develop and import all kind of needed software, and more specifically segmentation algorithms as plugins processes for Aïoli. Finally, if our toolset (geometric descriptors and segmentation) uses 3D data as a support for the information generation and the visualization of the results on the images, the inverse process can easily be envisaged thanks to the 2D/3D relationship. All processing performed on images can be easily projected onto the 3D model. Thus the image segmentation methods [START_REF] Vantaram | Survey of contemporary trends in color image segmentation[END_REF] could help in the 3D model segmentation and the machine learning methods generally applied to images rather than 3d data [START_REF] Ioannidou | Deep learning advances in computer vision with 3d data: A survey[END_REF] could allow to classify in a semi-automatic way the elements of the 3D model by detection in images (visible or geometric descriptors).

V. SEMI-ASSISSTED SELECTION

If one of Aïoli's specificities is to establish a link between semantic data and geometric data, it is essential to consider different annotation modalities adapted to these two fields: user-centered, when the annotation depends entirely on the user's field of expertise, automatic or assisted, when the annotation results from calculations. Thus, the existing semantic annotation tools are essentially manual since they aim to bring semantically qualified information to the 2D and 3D data of the project. For now, three tools are available on the platform: point selection, polygon drawing, or free form. These three modalities are mainly carried out on photographs. The next part introduces new tools for assisted annotation, based on the exploitation of segmentations and the implementation of graphical selection tools.

A. From multisegmentation to annotation layers

Today, the work performed on the geometric analysis of project data allows us to consider extending these methods with an assisted annotation tool. The latter will be based on the implementation of multi-segmentation plugins, like the one described above for plane detection. The principle is simple: a preliminary segmentation operated on a point cloud can automatically suggest annotations to the user in order to select regions with similar geometric characteristics. If the annotation proposal is not sufficient for the user, he can modify it in order to fit his objectives, thanks to existing graphic tools. Once the regions have been selected, they can be enriched via their layer description sheets.

B. Graphical selection

Aïoli's main annotation functionality is based on the observation of objects by experts. Starting from this principle, the most suitable media to collect annotations are the images (photos) of the object. These annotations, made using the graphic tools mentioned earlier, are not automatic. However, it is possible to rely on geometric descriptors in order to make precise annotations according to particular features. Some automatic selection tools available on large image editing software (Photoshop and others) would facilitate the detection and annotation of features on the Aïoli platform. Among these tools we are mainly interested in magnetic lassos and magic wands, allowing, starting from a simple parameter (threshold), to select in a more or less precise and effective way regions of interest, in particular those which are not visible on the photographs, but stand out on the geometric descriptors reproduced. A magic wand tool would not only allow to locate areas with close colors on the photos, but also to locate areas with close features on different textures. For example, this tool automatically selects a plane by applying it to the normal map, with the correct threshold parameter (Figure 6). But we can go further, by complexifying a bit the parameters. If we provide a function to the wand, we could find not only constant (± a threshold) data, but also predictable evolution, described by mathematical profiles. For instance, if we provide a function of the form f (x) = T hreshold + x, by applying this function to the X-axis, we can allow the user to automatically select curves with a single click. These tools offer the user a range of interactions according to his description objectives. These different tools are flexible enough to allow the ergonomic creation of regions for less experienced users, but also to offer precise annotation support to users eager to make the most of the tool's capabilities. These tools are designed to support and facilitate the user's annotation work in an interactive way. These two methods of automatic region extraction are not disjointed; we can indeed pre-establish annotations using a preliminary segmentation, and then complete and correct it manually or using graphical tools.

CONCLUSION

This paper describes an approach for the semantic annotation of heritage assets using a bijective relationship established by indexation, between a 3D representation of an object and all the associated oriented 2D images. Our objective is to establish a digital continuum linking the acquisition, interpretation and representation phases, by providing relevant annotation tools to help cultural heritage experts to enrich data collected on the study object. While this goal is reachable, a few questions still raise our attention. DISCUSSION AND FUTURE WORKS Our workflow must be simple and flexible, and must allow to compare the annotations made within collaborative and multidisciplinary frameworks. Our web architecture might render collaborative work on multiple devices easy, but it also raises the obvious problem of Internet access: without a connection (within archaeological sites for example), the application is not accessible, thus preventing in situ uses. Although the regular deployment of 4G antennas is gradually improving network coverage in remote areas, it is still worth considering alternatives. For instance, we are studying the feasibility of implementing an offline mode that allows users to pre-load their projects in order to perform their annotations locally. The data entered could then be updated at a later Internet connection. Unfortunately, this solution is not fully adapted for collaborative use, as it would require a postponedsynchronisation algorithm to allow multiple users to work simultaneously on the same offline project.

Furthermore, the core of Aioli is based on a multi-stereo correlation engine. Its limits are therefore intrinsically linked to the photogrammetric pipeline, both for calculation and visualization. These drawbacks in mind, users are required to take important precautions when acquiring raw data and respect the good practices related to this metrological method: choice of sensor, camera settings, covering, etc. On the other hand, any image or data compatible with the photogrammetry process is theoretically directly usable in Aioli. A long-term perspective would therefore be to work on these computation steps in order to extend our approach to different scientific imaging techniques (e.g. multispectral imagery).

Last but not least, the question of ergonomics and acceptability remains open: the research carried out so far demonstrates the possibility of creating new tools for the study of Cultural Heritage. However, a paradigm shift in this area will only occur if these tools are truly adapted to the needs and practices of experts, and sufficiently intuitive and flexible. The choice of semi-automatic methods for annotation was made for this reason: the results are certainly less direct, but this allows users to configure and correct the data obtained, whether by the segmentation parameters, or by refining the temporary regions made by magic wand and magnetic lasso tools. The Aioli collaborative cloud platform is already open for beta-testing, and all the tools detailed in this paper have been developed with the aim of being integrated soon. Following the various expert feedback, further developments are focused on modalities for visualization and data capture, in order to facilitate observation and annotation of barely perceptible elements, and to establish a more tangible link between the digital clone and the real object. Indeed, despite the presence of ever more precise digital data, the study of built heritage still requires many round trips with the real object, to contextualize the data and compare it with what can be perceived in situ. From this perspective, Augmented and Mixed Reality appears to be a high-potential medium, on which we are currently working, but which requires requires solving the problem of the alignment of a project's data with the real heritage object.

Fig. 1 .

 1 Fig. 1. Bijective relationship between images and point cloud

Fig. 3 .

 3 Fig. 3. Examples of reprojected geometrical descriptors

Fig. 4 .

 4 Fig. 4. User interface for the Aïoli platform

Fig. 5 .

 5 Fig. 5. Example of segmentation answer, and its reprojection on the images

Fig. 6 .

 6 Fig. 6. Magic wand possibilities on normal map and colored textures