Matthew L Spangler

Dominique Hazard
email: dominique.hazard@inrae.fr

Florence Plisson-Petit

Carole Moreno-Romieux

Stéphane Fabre

Laurence Drouilhet

Genetic Determinism Exists for the Global DNA Methylation Rate in Sheep

Keywords: epigenetics, blood, heritability, genetic correlation, QTL

Recent studies showed that epigenetic marks, including DNA methylation, influence production and adaptive traits in plants and animals. So far, most studies dealing with genetics and epigenetics considered DNA methylation sites independently. However, the genetic basis of the global DNA methylation rate (GDMR) remains unknown. The main objective of the present study was to investigate genetic determinism of GDMR in sheep. The experiment was conducted on 1,047 Romane sheep allocated into 10 halfsib families. After weaning, all the lambs were phenotyped for global GDMR in blood as well as for production and adaptive traits. GDMR was measured by LUminometric Methylation Analysis (LUMA) using a pyrosequencing approach. Association analyses were conducted on some of the lambs (n = 775) genotyped by using the Illumina OvineSNP50 BeadChip. Blood GDMR varied among the animals (average 70.7 ± 6.0%). Female lambs had significantly higher GDMR than male lambs. Inter-individual variability of blood GDMR had an additive genetic component and heritability was moderate (h 2 = 0.20 ± 0.05). No significant genetic correlation was found between GDMR and growth or carcass traits, birthcoat, or social behaviors. Association analyses revealed 28 QTLs associated with blood GDMR. Seven genomic regions on chromosomes 1, 5, 11, 17, 24, and 26 were of most interest due to either high significant associations with GDMR or to the relevance of genes located close to the QTLs. QTL effects were moderate. Genomic regions associated with GDMR harbored several genes not yet described as being involved in DNA methylation, but some are already known to play an active role in gene expression. In addition, some candidate genes, CHD1, NCO3A, KDM8, KAT7, and KAT6A have previously been described to be involved in epigenetic modifications. In conclusion, the results of the present study indicate that blood GDMR in domestic sheep is under polygenic influence and provide new insights into DNA methylation genetic determinism.

INTRODUCTION

Genetic selection in several livestock species has recently progressed to include heritable DNA polymorphisms (i.e., genomic selection) for the improvement of production traits [START_REF] Goddard | Genomic selection in livestock populations[END_REF]. However, recent studies also showed that non-genetic information responsible for phenotypic differences between animals can also be inherited across generations [START_REF] Johannes | Epigenome dynamics: a quantitative genetics perspective[END_REF][START_REF] Danchin | Beyond DNA: integrating inclusive inheritance into an extended theory of evolution[END_REF]. One such non-genetic inherited piece of information concerns epigenetic marks. Epigenetic modifications include biochemical modifications of DNA (methylation) or proteins tied to DNA (methylation and acetylation of histone proteins). These epigenetics processes act in an interrelated way to influence gene expression and hence phenotype in response to environmental conditions, thereby demonstrating that the eukaryotic genome can respond dynamically to changes in the environment to which every individual is exposed [START_REF] Geiman | Chromatin remodeling, histone modifications, and DNA methylation-how does it all fit together?[END_REF].

One of the key mechanisms in the regulation of gene expression in mammals is cytosine methylation at CpG dinucleotides, the most common DNA modification in eukaryotes. The importance of epigenetics in mammals and plants has been demonstrated in several studies and epigenetics contribute to the inheritance of traits of interest. In plants, DNA methylation in several differentially methylated regions was shown to explain 60-90% of the heritability of complex traits [START_REF] Cortijo | Mapping the epigenetic basis of complex traits[END_REF]. In humans, DNA methylation explained 7-11% of the body mass index while SNP explained 5-8% of the same trait [START_REF] Dick | DNA methylation and body-mass index: a genome-wide analysis[END_REF]. In livestock, current studies aim to identify epi-loci underlying traits of agronomical interest. For instance, in sheep, genome-wide analysis of DNA methylation profiles has already revealed numerous epi-loci associated with prolificacy or body size [START_REF] Cao | DNA methylation Landscape of body size variation in sheep[END_REF][START_REF] Zhang | Genome-wide analysis of DNA Methylation profiles on sheep ovaries associated with prolificacy using whole-genome Bisulfite sequencing[END_REF]. Epigenetic sources of inheritance are not currently taken into consideration in selection strategies in livestock production but these studies could help identify the origin of variability of production traits and later on, improve genetic progress.

Taking into account variations in epigenetic marks may enable better modeling of phenotypic variation and improving the precision of genetic evaluation of traits and breeds. In addition, taking epigenetic events into account may help not only understand but also improve adaptive processes in both plants and mammals. Particular epigenetic recombinant inbred lines of Arabidopsis thaliana showed highly variable DNA methylation, while being genetically very similar [START_REF] Latzel | Epigenetic diversity increases the productivity and stability of plant populations[END_REF]. Using these plant model lines, it was demonstrated an important role of epigenetics in adaptive process facing saline conditions for example (Kooke et al., 2015). Moreover, evolutionary studies showed that epigenetic inheritance may be important for natural selection and transmission of advantageous traits, for a review (see [START_REF] Donohue | THE epigenetics of adaptation: focusing on epigenetic stability as an evolving trait[END_REF]. Recent advances in epigenetic studies in natural populations addressed how epigenetic inheritance may influence adaptive evolution by focusing on epigenetic stability and inheritance itself as a potentially evolving trait [START_REF] Herman | How stable 'should' epigenetic modifications be? Insights from adaptive plasticity and bet hedging[END_REF][START_REF] Schlichting | Phenotypic plasticity and epigenetic marking: an assessment of evidence for genetic accommodation[END_REF]. These articles brought into focus the genetic and ecological basis of epigenetic stability and raised a number of questions including "is the stability of epigenetic marks heritable?" [START_REF] Turck | Natural variation in epigenetic gene regulation and its effects on plant developmental traits[END_REF] showed that DNA methylation has the potential to be mitotically and meiotically stable, whereas histone modification is involved in environmentally induced epigenetic regulation that may be reversible. These authors hypothesized that both the degree of DNA methylation and the stability of environmentally induced changes in gene expression via histone modification are associated with changes in the DNA sequence that can cause heritable variation in epigenetic regulation. New scientific questions arise concerning the relationship between genetics and epigenetics as part of the improvement of livestock species and their adaptation to changing production systems. In the present study, we used a complementary approach to those used in existing studies by considering the global methylation of genomic DNA as a novel quantitative phenotype to provide essential information on the processes involved in epigenetic inheritance. DNA methylation in differentially methylated regions representing a very small proportion of genomic DNA global methylation, we hypothesized that the stable fraction of DNA methylation is genetically determined and can be inherited across generations. To test this hypothesis, we quantified the global DNA methylation rate in blood samples collected from domestic sheep. Sheep not only has the advantage of being a species of agronomical interest, but also represents an interesting animal model to study adaptive processes. The aim of the present study was to estimate the genetic component of the global DNA methylation rate and to identify the genomic loci underlying the phenotype.

MATERIALS AND METHODS

Animals and Management

The experimental animals were Romane lambs, a fixed crossbreed between Romanov x Berrichon du Cher [START_REF] Ricordeau | La race ovine INRA 401: un exemple de souche synthétique[END_REF]. Data collected from a total of 1,047 male and female lambs over a period of 5 years were used in this study. The lambs originated from 10 rams, each ram was used for only 1 year. All the animals were born in April, identified at birth using electronic ear tags, and reared outdoors with their dams under extensive conditions in the experimental "la Fage" farm on 280 ha of rangeland (Causse du Larzac, Roquefort sur Soulzon, South of France). The farming system, management and environment characteristics are detailed in [START_REF] Gonzalez-Garcia | Circannual body reserve dynamics and metabolic profile changes in Romane ewes grazing on rangelands[END_REF]. The lambs were weighed regularly from lambing to weaning to estimate average daily gain. All the lambs were weaned at 85 ± 4 days of age. Two weeks after weaning, experimental lambs were individually exposed to two behavioral tests. Tests of each animal were all performed on the same day, with a total of around 35 lambs tested per day on 10-11 consecutive test-days per year. After the behavioral tests (i.e., approximately 1 month after weaning), lambs were blood sampled for genotyping and determination of the global DNA methylation rate.

Global DNA Methylation Rate

The global DNA methylation rate (GDMR) was measured from whole blood samples taken from sheep. Genomic DNA was extracted from the blood samples using a high-salt extraction method [START_REF] Roussot | AFLP linkage map of the Japanese quail Coturnix japonica[END_REF]. Methylation analyses were performed using the LUMA assay (Karimi et al., 2006a,b). Briefly, this method relies on the use of methylation-sensitive and -insensitive restriction endonucleases: HpaII and MspI, respectively. The target sequence for both enzymes is CCGG, HpaII is not able to cut if the internal cytosine is methylated (C m CGG), whereas MspI cuts the restriction site whatever the methylation status. Moreover, EcoRI (recognition site: GAATTC) was used as an internal control for normalization of DNA amount. For each sample, DNA was independently digested by EcoRI + HpaII and EcoRI + MspI restriction enzymes (New England Biolabs) and then analyzed using a Q24 Pyromark sequencer and Q24 1.0.10 software (Qiagen). The dispensation order for nucleotides was GTGTCACATGTGTG. Methylation rate was calculated from peak heights as [1 -[(HpaII(G)/ EcoRI_Hpa(T))/(MspI(G)/ EcoRI_Msp(T))] × 100]. The same internal control DNA sample was used for each pyrosequencing run. Finally, considering complete digestions, GDMR integrates data from nearly six million CCGG sites detected in silico on sheep reference genome OAR v3.1.

Behavioral Traits

Two behavioral tests were used to evaluate sociability for conspecifics and reaction to a human. A complete description of the devices and test procedure used for behavioral measurements are given in [START_REF] Hazard | Genotype by environment interactions for behavioral reactivity in sheep[END_REF]. Briefly, the arena test (AT) aimed to evaluate sociability for conspecifics in three successive phases: (1) attraction for conspecifics with visual contact (Arena test phase 1, AT1), (2) reactivity to social separation from conspecifics by preventing visual contact between the tested lamb and conspecifics (Arena test phase 2, AT2), and (3) attraction for conspecifics in presence of a motionless human standing in front of conspecifics (Arena test phase 3, AT3). The variables used were the frequency of low-pitched bleats (i.e., the lamb bleated with its mouth closed) recorded during AT1 (AT1-LBLEAT), and the frequency of high-pitched bleats (i.e., the lamb bleated with its mouth open) recorded during AT2 (AT2-HBLEAT). Locomotor activity was assessed by measuring the number of zones crossed during AT2 (AT2-LOCOM). Vigilance postures (i.e., animal motionless, head in an upright position, and ears perpendicular to the head) were measured during AT2 (AT2-VIGIL). The time spent in each virtual zone was recorded and the ewe's proximity score to conspecifics during AT3 (AT3-PROX) was calculated by weighting the time spent in each zone.

The corridor test (CT) aimed to evaluate the reaction to a human. The second phase of the CT was used to evaluate the lamb's reaction to a walking human. The mean distance between the tested lamb and the walking human was recorded (CT2-DIST).

Zootechnical Traits

At birth, the type of lamb birthcoat was graded on a scale of one to nine based on composition and structure of the coat, i.e., grade 1, single wooly coat with no halo-or coarse hair, hair length < 10 mm; grade 9, double hairy coat (coarse hair mixed with fine down) > 25 mm in length. Further details on this trait are provided in [START_REF] Allain | Importance of birthcoat for lamb survival and growth in the romane sheep breed extensively managed on rangelands[END_REF]. Lambs were weighed at birth (birth weight, kg) and weaning (weaning weight, kg). Lambs were also weighed regularly from lambing to weaning to estimate average daily gain. Growth rates (i.e., average daily gain, ADG, in g) were estimated in all lambs from birth to 30, from 30 to 60 and from 60 to 90 days of age (called ADG0-30, ADG30-60, ADG60-90, respectively). In addition, some carcass traits were measured at slaughter only for males, including dressing yield (DY, percentage), conformation score (CONF, score 0-6), compactness (COMP) (i.e., width/length ratio of the carcass, percentage), external fat score (FAT Score, score 0-9), and back fat depth at the 12th rib (FAT depth, mm).

Statistical Analysis

Analyses of variance were performed to assess the significant fixed effects affecting each measured trait. The tested effects were sex (male or female), born and reared litter size, age of the dam, age and weight of the lambs and year of measurement. The litter size effect was classified according to the number of lambs born and suckled (class 1, ewes lambing and suckling singletons; class 2, ewes lambing twins but suckling only one lamb; class 3, ewes lambing and suckling twins; and class 4, ewes lambing and suckling more than two lambs). The age of the dam effect included ewes that were one, two, and three or more than 3 years old (classes 1, 2, and 3, respectively). The year of measurement effect covered the 5 years over which data were collected. The GLM procedures of the SAS R software (SAS Institute Inc, 1999) were applied to the variables to test the fixed effects and first order interactions. The factors of variation showing significant effect (i.e., P< 0.05) on the considered trait were included in subsequent genetic analyses. The fixed effects considered depended on the trait. In addition, for the GDMR trait, we investigated differences between lamb families by testing a family effect (10 classes), nested in the year effect, according to the GLM procedure describe above. This was only used for GDMR because the genetic variability of the other traits included in the present study has already been described [START_REF] Bibé | Genetic parameters of growth and carcass quality of lambs at the French progeny-test station Berrytest[END_REF][START_REF] Allain | Importance of birthcoat for lamb survival and growth in the romane sheep breed extensively managed on rangelands[END_REF][START_REF] Hazard | Identification of QTLs for behavioral reactivity to social separation and humans in sheep using the OvineSNP50 BeadChip[END_REF].

Genetic Analyses

The (co)variance components of each trait were estimated by restricted maximum likelihood applied to an animal model using ASREML 3.0 software [START_REF] Gilmour | ASReml User Guide Release 3[END_REF]. For each trait, the model included the appropriate fixed effects (i.e., sex, born and reared litter size, age of the dam, year of measurement) and a direct additive genetic effect of the animal considered as random effect. For some traits analyzed, a litter or a dam permanent environmental effect was added as random effect in the model based on present results or the results of previous studies that described an appropriate model for the trait concerned [START_REF] Bibé | Genetic parameters of growth and carcass quality of lambs at the French progeny-test station Berrytest[END_REF][START_REF] Allain | Importance of birthcoat for lamb survival and growth in the romane sheep breed extensively managed on rangelands[END_REF][START_REF] Hazard | Identification of QTLs for behavioral reactivity to social separation and humans in sheep using the OvineSNP50 BeadChip[END_REF]. For GDMR, weight and growth traits, the dam was considered as a permanent environmental effect. For birthcoat, the animal's litter was considered as a permanent environmental effect. The following complete animal mixed model was fitted:

y = Xβ + Z a a + W c c + e (1)
where y is the vector of observations for the trait(s) being analyzed, β is the vector of appropriate fixed effects (sex of lambs, litter size born and suckled, age of the dam and year of measurement), a is the vector of random genetic effects and c is the vector of permanent environmental effects, when appropriate, with incidence matrices X, Z a , and W c , respectively, and e is the vector of residual effects. Univariate analyses were performed to estimate the variance of each trait. Bivariate analyses were performed to estimate the genetic correlations between the GDMR and the other traits using the same model as that used in univariate analyses.

Genomic Analyses

Genomic analyses were only done for the GDMR trait. Nine of the 10 families were used (i.e., family 10 was not genotyped). Among the 1,047 lambs phenotyped for GDMR, 800 lambs were genotyped (i.e., after filtering an individual call rate ≥ 98%) as well as their nine respective sires using the Illumina OvineSNP50 BeadChip (i.e., 54,241 SNPs). Outlier individuals were removed, and genomic analyses were performed using 775 phenotyped and genotyped lambs. SNP quality was checked as described by [START_REF] Hazard | Identification of QTLs for behavioral reactivity to social separation and humans in sheep using the OvineSNP50 BeadChip[END_REF] (i.e., SNP call rates < 97%, a minor allele frequency < 1%, inconsistent for Hardy-Weinberg disequilibrium). Sex chromosomes were not included for analysis. Finally, 40,725 autosomal SNPs were retained for QTL analyses. Genome wide association analysis was performed of the whole population genotyped using joint analysis considering simultaneously linkage association and linkage disequilibrium (LDLA) to take advantage of both family structure and linkage disequilibrium [START_REF] Legarra | Linear models for joint association and linkage QTL mapping[END_REF] using QTLMAP software [START_REF] Gilbert | QTLMAP, a software for QTL detection in outbred populations[END_REF]. The LDLA method consisted in interval mapping and LDLA model considered the sire haplotypes effects of the LD model in addition to sire QTL effect. A haplotype size of four SNP was used, and when haplotype frequency was less than 1%, haplotypes were assigned to a rare haplotype group. The chromosome-wise p-values were estimated assuming that, conditional on the QTL position, the likelihood ratio test statistics followed a χ2 -distribution with k degrees of freedom, k being the number of genetic effects [START_REF] Piepho | A quick method for computing approximate thresholds for quantitative trait loci detection[END_REF]. In our study, k was equal to the number of haplotypes plus the number of families for LDLA. Genome-wise p-values were obtained by correcting for multiple testing assuming that the number of independent tests was equal to the number of chromosomes analyzed (i.e., 26 independent tests corresponding to the 26 autosomes) [START_REF] Benjamini | Controlling the false discovery ratea pratical and powerful approach to multiple testing[END_REF][START_REF] Knott | Multiple marker mapping of quantitative trait loci in a cross between outbred wild boar and large white pigs[END_REF]. Confidence intervals were determined using the "2-LOD drop-off " criterion and assuming 1 LOD = 4.61 LRT [START_REF] Lynch | Genetics and Analysis of Quantitative Traits[END_REF].

For each QTL, the confidence interval expanded by 2 Mb on each side was used to extract gene annotation by using the Biomart tool from Ensembl release 101 of the sheep reference genome OAR v3.11 . Functional annotations for epigenetic modifiers were extracted from dbEM 2 [START_REF] Nanda | DbEM: a database of epigenetic modifiers curated from cancerous and normal genomes[END_REF] and from the review about epigenetic modifiers by [START_REF] Feinberg | Epigenetic modulators, modifiers and mediators in cancer aetiology and progression[END_REF].

RESULTS

Phenotypes and Genetic Analyses

Descriptive statistics of GDMR are summarized in Table 1. The mean value of GDMR was 70.71% and the coefficient of variation was 8.40%. Analysis of the gender effect indicated that female lambs had higher GDMR than male lambs. Marked differences were observed between extreme values within a gender. Phenotypic variability of male lambs was higher than female lambs. No significant effects of litter size, age of the dam, age of the lamb and weight (at birth or at weaning) were found, but a significant effect of family (nested in the year effect) was found on GDMR. One family had a low GDMR at 66.60% (family 3) contrasting with four families having a higher GDMR ranging from 71.95 to 74.11% (families 1, 2, 7, and 8). Four families had intermediate GDMR ranging from 68.86 to 70.85% (families 4, 5, 6, 9, and 10). The difference between extreme families reached 1.25 SD. Differences in phenotypic variability were observed between families.

The estimated variance components of GDMR are listed in Table 2. The estimated heritability for GDMR was moderate (0.20 ± 0.05) and the permanent environment effect of the dam was very low. The additive genetic coefficient of variation for GDMR was also low (3.5%).

Descriptive statistics and genetic parameters for additional traits recorded in the lambs used in the present study are reported in Table 3. Heritability estimates for behavioral, zootechnical and carcass traits were generally moderate to high (0.20 ± 0.05 to 0.74 ± 0.07). Only low-pitched bleats and conformation score traits showed low heritability. No genetic correlation was found between GDMR and behavioral, birthcoat or zootechnical traits.

Genomic Analyses

The significant QTLs found using LDLA and reaching genomewide (GW) or chromosome-wide (CW) significance thresholds are reported in Table 4 and Figure 1. Four haplotype-trait associations reached the GW significance thresholds and 16 haplotype-trait associations reached the 1% CW significance threshold. Eight significant associations reaching the 5% CW threshold were also detected. Associations found in LDLA mapped to 14 chromosomes. The four associations that reached the GW threshold of significance localized on chromosomes 1 (OAR1, 13.20 and 179.50 Mb) and 17 (OAR17,40.14 and 55.54 Mb). Associations reaching the 1% CW significance threshold were mapped on chromosomes OAR3, 5, 13, 14, 16, 17, 18, 21, 23, 24, and 26. Effects of QTLs found in LDLA ranged from 0.68 to 3.20% of the total phenotypic variance of GDMR.

By crossing information from sheep reference genome annotation (OAR v3.1), the database of epigenetic modifiers (dbEM, [START_REF] Nanda | DbEM: a database of epigenetic modifiers curated from cancerous and normal genomes[END_REF] and those reported by [START_REF] Feinberg | Epigenetic modulators, modifiers and mediators in cancer aetiology and progression[END_REF], we highlighted CHD1, NCO3A, KDM8, KAT7, and KATA6 as the most obvious candidate genes from each QTL (+2 Mb each side of the confidence interval), likely to play a role in the GDMR variation (Table 4).

DISCUSSION

The global DNA methylation rate measured by LUMA reached 70% in sheep blood. This result is consistent with global DNA methylation percentages reported in human blood analyzed by LUMA (for a review, see [START_REF] Soriano-Tarraga | DNA isolation method is a source of global DNA methylation variability measured with LUMA. experimental analysis and a systematic review[END_REF]. The value of GDMR in the present study is also in the range of calculated DNA methylation either using the global CpG methylation rate (50-55%) measured by reduced representation bisulfite sequencing (RRBS) in sheep muscle [START_REF] Couldrey | Genome-wide DNA methylation patterns and transcription analysis in sheep muscle[END_REF] or the total methyl content measured by high performance liquid chromatography reported in somatic tissues (80%) in Significant QTLs reaching the chromosome-wide or the genome-wide thresholds are listed. a *p < 5%; **p < 1%; ***p < 0.1%. b Percentage of phenotypic variance explained by the haplotype with the most effect (absolute value). c SNPs flanking the QTL with significant association. CW, chromosome wide; GW, genome wide; Position and genes are based on the reference genome Oar_v3.1.

Frontiers in Genetics | www.frontiersin.org , 2, 3, 5, 11, 13, 14, 16, 17, 18, 21, 23, 24, and 26. The horizontal lines indicate the 5, 1, and 0.1% genome wide thresholds. One percent chromosome wide thresholds were OAR3,47.6;OAR5,37.4;OAR13,42.6;OAR14,35.9;OAR16,38.7;OAR17,42.3;OAR18,40.6;OAR21,37.3;OAR23,45.7;OAR24,42.8,and OAR26,47.1. 5% chromosome wide thresholds were OAR1, 45.5; OAR2, 45.0; OAR11,41.3;OAR13,38.3;OAR24,38.4. humans [START_REF] Ehrlich | Amount and distribution of 5-methylcytosine in human DNA from different types of tissues or cells[END_REF]. The large range observed for DNA methylation rate across the genome is likely due to the methods used for quantification of DNA methylation and/or the tissues studied. Indeed, in their review, Soriano-Tarraga et al. ( 2013) reported that variation in global DNA methylation in human blood analyzed by LUMA ranged from 52 to 78% depending on the DNA isolation method used. In addition, [START_REF] Doherty | Exploring genorne wide bisulfite sequencing for DNA rnethylation analysis in livestock: a technical assessment[END_REF] reported that RRBS libraries are biased because they contain promoter regions which have high CG content known to be largely devoid of DNA methylation. Consequently, RRBS libraries are expected to display lower methylation on average across genome than unbiased libraries.

Variability was observed for the GDMR phenotype recorded in this study. Part of this variability originated from differences between male and female lambs. Sex differences in the methylome have already been reported in humans. In human liver, females displayed higher average DNA methylation in the X-chromosome whereas males presented higher DNA methylation in autosomes in this tissue [START_REF] García-Calzón | Sex differences in the methylome and transcriptome of the human liver and circulating HDL-cholesterol levels[END_REF]. Genome-wide sex differences in locus-specific DNA methylation across autosomes have also been reported in other human tissues such as blood [START_REF] Singmann | Characterization of whole-genome autosomal differences of DNA methylation between men and women[END_REF]. In the present study, the positions of DNA methylation across genome were unknown, so we are unable to conclude whether gender difference in DNA methylation found in lambs was of similar magnitude on the X-chromosome and autosomes. In addition, it should be kept in mind that we investigated DNA methylation in young lambs (4 months of age), not in adult sheep. A study investigating gender differences in DNA methylation in human blood at birth reported that CpG located on autosomes, as well as DMR, were hypermethylated in girls compared to boys [START_REF] Yousefi | Sex differences in DNA methylation assessed by 450 K BeadChip in newborns[END_REF]. Because gender differences in DNA methylation across autosomes in human blood varied between newborns and adults, we cannot exclude that a sex by age interaction could affect DNA methylation in sheep blood.

We hypothesized that GDMR has a genetic component in sheep. The variance component decomposition analysis of GDMR revealed moderate heritability (20%). To our knowledge, this is the first heritability estimate for GDMR in livestock species. Estimation of heritability of DNA methylation levels in human whole blood at thousands of sites has shown that heritability estimates vary across the genome but is on average 20%, consistently with our estimate in sheep blood [START_REF] Van Dongen | Genetic and environmental influences interact with age and sex in shaping the human methylome[END_REF]. Additional studies in humans also indicated that the average twin based heritability of DNA methylation across genome-wide CpGs varied between 5 and 19% depending on the tissue (i.e., h 2 = 5, 7, and 12% in placenta, human umbilical vascular endothelial cells and cord blood mononuclear cells from newborns twins, respectively, h 2 = 18% in peripheral blood lymphocytes from adolescent twins [START_REF] Gordon | Neonatal DNA methylation profile in human twins is specified by a complex interplay between intrauterine environmental and genetic factors, subject to tissue-specific influence[END_REF][START_REF] Mcrae | Contribution of genetic variation to transgenerational inheritance of DNA methylation[END_REF]). Interestingly, [START_REF] Shah | Genetic and environmental exposures constrain epigenetic drift over the human life course[END_REF] reported that DNA methylation in whole blood remained stable over the human lifetime, as indicated by the high correlation (0.68) between DNA methylation repeatability in older people (median age 70-90) and heritability estimated in teenagers.

Genetic relationships between GDMR and production or adaptive traits were also investigated in the present study. We did not find any relationship between GDMR and production or adaptive traits in domestic sheep. Our results contrasted with the negative correlation between total DNA methylation level and daily growth rate in fish embryos [START_REF] Ou | The DNA methylation level is associated with the superior growth of the hybrid fry in snakehead fish (Channa argus x Channa maculata)[END_REF]. To our knowledge, no published works provide evidence of a relationship between GDMR in blood and the specific adaptive traits used in the present study (social behaviors and birthcoat). Some relationships would be expected, since several studies suggested epigenetic influence on adaptive processes in mammals [START_REF] Massicotte | DNA methylation: a source of random variation in natural populations[END_REF][START_REF] Lea | Resource base influences genome-wide DNA methylation levels in wild baboons (Papio cynocephalus)[END_REF]. Nevertheless, only a few adaptive traits were explored in this study and we cannot exclude the possibility that such relationships exist for other adaptive traits or in older animals. Additionally, regardless of the absence of relationships between GDMR and adaptive traits in our study, we cannot exclude, by considering individually DNA methylation sites, to find epi-loci associated with production or adaptive traits investigated in the present study. Indeed, epiloci have previously been reported to be associated with animal behaviors [START_REF] Champagne | Epigenetic mechanisms mediating the long-term effects of maternal care on development[END_REF][START_REF] Roth | Epigenetic mechanisms in the development of behavior: advances, challenges, and future promises of a new field[END_REF][START_REF] Verhulst | Evidence from pyrosequencing indicates that natural variation in animal personality is associated with DRD4 DNA methylation[END_REF][START_REF] Seebacher | Epigenetics of social behaviour[END_REF].

Even if the greater attention paid to epigenetics in farm animals has led to an increasing number of GWA epigenetic studies, genomic studies of DNA methylation rate had remained to be performed. In our study, we were particularly interested in detecting QTLs for GDMR in domestic sheep. The QTL analysis for GDMR resulted in mapping 28 QTLs on 14 chromosomes, with four regions (on OAR1 and OAR17) reaching the genome wide level of significance. Each detected QTL explained less than 3% of the phenotypic variance. Our experimental design enabled to reach a power of 90% for the detection of a QTL explaining 8% of the phenotypic variance while the power was 20% for a QTL that explains 1% of phenotypic variance [START_REF] Luo | Detecting linkage disequilibrium between a polymorphic marker locus and a trait locus in natural populations[END_REF]. Consequently, the undetected QTLs likely had a limited effect on GDMR. Thus, the low proportion of variance explained by each QTL supports the hypothesis that GDMR in domestic sheep is under polygenic influence and unlikely under the control of a major gene. Nevertheless, several QTLs mapped in the present study could probably act in combination to account for substantial genetic variation of GDMR.

To better understand the underlying molecular mechanisms associated with these QTLs, we looked for possible overlapping between the location of the identified QTLs (confidence interval enlarged by 2 Mb each side) and the location of annotated genes relevant to epigenetics. Despite the large number of genes located close to the QTL regions identified, only few genes seem relevant as epigenetic modifiers through their action on methylation or acetylation processes on DNA or histones. Indeed, by crossing the list of positional genes in enlarged QTL regions with that of known epigenetic modifiers either from the dbEM database [START_REF] Nanda | DbEM: a database of epigenetic modifiers curated from cancerous and normal genomes[END_REF] or those listed in [START_REF] Feinberg | Epigenetic modulators, modifiers and mediators in cancer aetiology and progression[END_REF], only CHD1, KDM8, NCO3A, KAT7, and KATA6 appeared as the most obvious candidate genes. CHD1 gene encodes the chromodomain-helicase-DNA-binding protein 1 (OAR5: 95.4 Mb), that functions as substrate recognition component of the transcription regulatory histone acetylation complex SAGA. Mutations in CHD1 are frequently associated with prostate cancers [START_REF] Berger | The genomic complexity of primary human prostate cancer[END_REF]. Functional analyses have shown that CHD1 was involved in transient DNA methylation and several loci were hypermethylated in CHD1 deleted strains of Neurospora crassa [START_REF] Belden | CHD1 remodels chromatin and influences transient DNA methylation at the clock gene frequency[END_REF]. KDM8 (Lysin demethylase 8 or Jumonji C domain-containing demethylase 5, JDMD5; OAR24: 24.8 Mb) is a histone H3 demethylase with specificity for Lys-36. This demethylase promotes homologous recombination [START_REF] Amendola | JMJD-5/KDM8 regulates H3K36me2 and is required for late steps of homologous recombination and genome integrity[END_REF], cell proliferation [START_REF] Hsia | KDM8, a H3K36me2 histone demethylase that acts in the cyclin A1 coding region to regulate cancer cell proliferation[END_REF] and protects from nerve demyelination [START_REF] Fuhrmann | Miz1 controls schwann cell proliferation via H3K36(me2) demethylase Kdm8 to prevent peripheral nerve demyelination[END_REF]. Interestingly, a second gene encoding an important histone lysine demethylase, KDM2B (Lysine-specific demethylase 2B, OAR17: 53.23 Mb), was found close to one of the enlarged QTL regions (OAR17: 55.4-55.6 Mb). In this region, we also identified SETD1B gene (Histone-lysine N-methyltransferase SETD1B, 53.06 Mb) which encodes a component of a histone methyltransferase complex that specifically methylates Lys-4 of histone H3 and is responsible for the epigenetic control of chromatin structure and gene expression. A specific hypermethylation signature was associated with loss of function mutations in the SETD1B gene [START_REF] Krzyzewska | A genome-wide DNA methylation signature for SETD1Brelated syndrome[END_REF]. The NCOA3 gene (nuclear receptor coactivator 3; OAR13: 32.0 Mb) encodes a nuclear receptor coactivator with histone acetyltransferase activity. It particularly interacts with other transcriptional activators such p300/CBP-associated factor and CREB binding protein (CREBBP) as part of a multi-subunit coactivation complex. The CREBBP gene, also located near a significant QTL position (OAR 24: 3.1 Mb) has been recently reported to be involved in the crosstalk between DNA methylation and histone acetylation [START_REF] Zhang | Crosstalk between DNA methylation and histone acetylation triggers GDNF high transcription in glioblastoma cells[END_REF]. Finally, two lysine acetyltransferase coding genes of the MYST family, KAT6A (lysine acetyltransferase 6A, also named MOZ, OAR26: 35.3 Mb) and KAT7 (lysine acetyltransferase 7, also named HBO1, OAR11: 36.4 Mb) has been reported to play a key role in acetylation of lysine residues in histone H3 and/or H4 [START_REF] Kitabayashi | Activation of AML1-mediated transcription by MOZ and inhibition by the MOZ-CBP fusion protein[END_REF]. Interestingly, KDM8, NCOA3 and KAT7 act as androgen or estrogen receptor co-regulators [START_REF] Sharma | Androgen receptor interacts with a novel MYST protein, HBO1[END_REF][START_REF] Vija | Expression and characterization of androgen receptor coregulators, SRC-2 and HBO1, during human testis ontogenesis and in androgen signaling deficient patients[END_REF][START_REF] Wagner | NCOA3 is a selective co-activator of estrogen receptor alpha-mediated transactivation of PLAC1 in MCF-7 breast cancer cells[END_REF][START_REF] Wang | KDM8/JMJD5 as a dual coactivator of AR and PKM2 integrates AR/EZH2 network and tumor metabolism in CRPC[END_REF]. Polymorphisms that alter the function of these genes could explain that female lambs had higher GDMR than male lambs.

The level and state of histone acetylation regulated by histone acetyl transferase enzymes are widely reported to be an epigenetic regulation of gene expression [START_REF] Kouzarides | Chromatin modifications and their function[END_REF]. Therefore, we cannot exclude the possibility that global DNA methylation in sheep is associated with histone acetylation, since relationships between DNA methylation and histone deacetylation have already been described [START_REF] Geiman | Chromatin remodeling, histone modifications, and DNA methylation-how does it all fit together?[END_REF][START_REF] Vaissière | Epigenetic interplay between histone modifications and DNA methylation in gene silencing[END_REF]. Several other coding genes, located close to QTLs we identified, have been reported to be involved in the regulation of transcription or regulation of chromatin architecture. Although no scientific evidence has suggested their involvement in DNA methylation to date, these genes could mediate this process through their active role in the coregulatory complexes that associate epigenetic modifiers to modulate transcription.

CONCLUSION

Inter-individual variability in global DNA methylation rate in blood has an additive genetic component in sheep and heritability was moderate (h 2 = 0.20 ± 0.05). This opens the way for a possible genetic selection of this trait and create experimental divergent lines to investigate further the relationship between GDMR and adaptive traits as already evidenced in plants [START_REF] Johannes | Assessing the impact of transgenerational epigenetic variation on complex traits[END_REF][START_REF] Kooke | Epigenetic variation contributes to environmental adaptation of Arabidopsis thaliana[END_REF]. Moreover, this work reported the first SNP-based QTL detection study for GDMR as a quantitative trait in livestock species. The evidence of 28 QTLs associated with blood GDMR, each explaining a small proportion of the phenotypic variance (between 1 and 3%), most likely indicated a polygenic determinism of this trait. The further identification of genes and their polymorphisms underlying the global DNA methylation rate assessed in our study will paves the way for a deeper understanding the genetic component of such a trait.

FIGURE 1 |

 1 FIGURE 1 | Chromosome plots of the likelihood ratio test values obtained for global DNA methylation rate by LDLA analysis. The likelihood ratio test (LRT) is plotted against SNP haplotype positions (four consecutive SNPs) for chromosomes 1, 2, 3, 5, 11, 13, 14, 16, 17, 18, 21, 23, 24, and 26. The horizontal lines indicate the 5, 1, and 0.1% genome wide thresholds. One percent chromosome wide thresholds wereOAR3, 47.6; OAR5, 37.4; OAR13, 42.6; OAR14, 35.9; OAR16, 38.7; OAR17, 42.3; OAR18, 40.6; OAR21, 37.3; OAR23, 45.7; OAR24, 42.8, and OAR26, 47.1. 5% chromosome wide thresholds were OAR1, 45.5; OAR2, 45.0; OAR11, 41.3; OAR13, 38.3; OAR24, 38.4.

TABLE 1 |

 1 Number of records, mean (± SD), minimum, maximum, quartiles, and least squares means for global DNA methylation rate in the whole population and for each sex and each family.

			Number of lambs	Mean (±SD)	Min	Max	Q1	Q3	P-value	LSMeans
	Population		1047	70.71 (5.97)	23.08	87.94	68.72	74.24
	Sex	Male	426	68.75 (7.21)	23.08	80.63	66.54	72.60	***	69.12
		Female	621	72.05 (4.47)	40.61	87.94	70.01	74.85		72.17
	Family (year)	1	93	73.40 (4.37)	55.82	81.28	71.17	76.81	***	73.05
		2	77	74.73 (3.55)	60.56	80.79	72.82	77.21		74.11
		3	116	66.94 (8.15)	27.04	75.34	65.34	71.46		66.63
		4	91	71.27 (6.7)	23.08	79.85	69.83	75.05		70.85
		5	115	68.95 (4.95)	37.96	77.06	66.24	72.05		68.86
		6	117	69.40 (5.62)	43.39	87.94	66.82	72.18		69.29
		7	101	73.00 (3.89)	55.01	77.89	71.30	75.87		72.57
		8	108	72.20 (3.31)	61.02	78.83	70.52	74.43		71.95
		9	123	69.92 (5.60)	41.62	79.56	68.65	73.03		69.81
		10	106	69.64 (6.73)	28.03	77.22	68.19	73.32		69.33
	***P-value < 0.001. Min, minimum; max, maximum; Q1, first quartile; Q3, third quartile; LSMeans, least squares means; Family effect was nested in the year effect.

TABLE 2 |

 2 Estimates of heritability, repeatability, permanent and residual effects (±S.E.) for the global DNA methylation rate.

				Component		Repeatability (r 2 )	Total σ p	2
		n	Genetic (h 2 )	Permanent (m 2 )	Residual (e 2 )
	GDMR	1,047	0.20 ± 0.05	0.02 ± 0.03	0.78 ± 0.05	0.22 ± 0.05	31.3 ± 1.5
	h 2 , m 2 , e 2 , proportion of total phenotypic variance attributed to additive genetic, permanent, and residual effects, respectively, Total σ p	2 , total phenotypic variance;
	r 2 = h 2	+ m 2 ; n, number of animals.

TABLE 3 |

 3 Number of records, mean (±SD), estimates of heritability, and permanent environmental effects for behavioral, zootechnical and carcass traits, and genetic correlations between these traits and the global DNA methylation rate.

		n	Mean (±SD)	h 2	p 2	r g
	AT1-LBLEAT	1,047	0.28 (0.28)	0.17 (0.07)	ne	<0.05 (<0.01)
	AT2-HBLEAT	1,047	2.91 (1.39)	0.46 (0.07)	ne	0.05 (0.03)
	AT2-LOCOM	1,017	19.44 (8.5)	0.26 (0.07)	ne	<0.05 (<0.01)
	AT2-VIGIL	1,017	20.68 (9.81)	0.32 (0.07)	ne	<0.05 (<0.01)
	AT3-PROX	990	26.26 (16.89)	0.27 (0.06)	ne	<0.05 (<0.01)
	CT2-DIST	1,047	5.37 (1.14)	0.20 (0.06)	ne	<0.05 (<0.01)
	Birthcoat (score 0-9)	1,032	6.06 (3.38)	0.74 (0.07)	0.15 (0.04)	<0.05 (<0.01)
	Birth weight (kg)	1,047	3. 9 (0.8)	0.26 (0.09)	0.27 (0.05)	<0.05 (<0.01)
	Weaning weight (kg)	1,043	23.3 (4.2)	0.36 (0.09)	0.20 (0.04)	<0.05 (<0.01)
	ADG0-30 (g)	1,047	360.3 (73.04)	0.30 (0.08)	0.24 (0.04)	<0.05 (<0.01)
	ADG30-60 (g)	1,047	227.4 (50.5)	0.34 (0.08)	0.06 (0.04)	<0.05 (<0.01)
	ADG60-90 (g)	1,046	216.0 (58.0)	0.32 (0.07)	0.10 (0.04)	<0.05 (<0.01)
	DY (%)	386	42.07 (2.19)	0.41 (0.15)	ne	<0.05 (<0.01)
	CONF (score 0-6)	390	3.44 (0.8)	0.09 (0.08)	ne	<0.05 (<0.01)
	COMP (%)	388	32.26 (1.46)	0.58 (0.20)	ne	<0.05 (<0.01)
	FAT score (score 0-9)	390	5.07 (1.39)	0.22 (0.12)	ne	<0.05 (<0.01)
	FAT depth (mm)	389	1.87 (1.36)	0.36 (0.13)	ne	<0.05 (<0.01)
	AT1/2/3, arena test phase 1/2/3; CT2, corridor test phase2; LBLEAT, low bleats; HBLEAT, high bleats; LOCOM, locomotion; VIGIL, vigilance behavior; PROX, proximity
	score; DIST, mean distance between tested animal and a human; ADG, average daily gain; DY, dressing yield; CONF, conformation score; COMP, compactness. h 2 ,
	p

2

, proportion of total phenotypic variance attributed to additive genetic or permanent environmental effects, respectively, r g , genetic correlations; standard errors are in parentheses; ne, not estimated.

TABLE 4 |

 4 List of QTLs detected in joint analysis (linkage disequilibrium linkage analysis) associated with the global DNA methylation rate.

	Chr.	Significance	level a	QTL effect (%) b	Confidence interval (Mb)	Flanking markers c		Annotated protein coding genes	Epigenetic modifier
								within confidence interval	candidate
	1	GW	*	2.48	13.1-13.4	OAR1_13033760.1	s62110.1	RRAGC, GJA9, RHDB
	1	CW	*	2.75	69.3-69.5	OAR1_74015424.1	OAR1_74065628.1	MTF2
	1	CW	*	3.0	138.0-138.3	OAR1_149183890.1	OAR1_149346875.1	-
	1	GW	*	2.4	179.4-179.6	OAR1_193623565.1	OAR1_193664187.1	-
	2	CW	*	1.44	199.7-199.9	OAR2_211501156.1	OAR2_211573792.1	-
	3	CW	**	1.0	114.5-114.7	OAR3_122117752.1	OAR3_122207054.1	-
	5	CW	**	0.83	96.6-96.9	OAR5_105186382.1	OAR5_105374983.1	-	CHD1
	11	CW	*	2.95	10.5-10.9	s74938.1	OAR11_10495137.1	MAD13, INTS2, BRIP1
	11	CW	*	0.85	37.4-37.6	OAR11_39854191.1	OAR11_39948154.1	HOXB1/2/3/4/5/6/7/9	KAT7
	13	CW	*	2.9	17.8-18.5	OAR13_20537279.1	OAR13_20578061.1	PARD3
	13	CW	*	2.7	32.1-32.3	OAR13_35499997.1	OAR13_35648372.1	NSUN6, EPC1
	13	CW	**	1.17	74.5-74.7	s65442.1	s62908.1	CDH22, SLC35C2, ELMO2	NCOA3
	14	CW	**	2.45	7.7-7.9	s16689.1	OAR14_8081413.1	CMIP, PLCG2, SNORA70
	14	CW	**	2.06	10.7-10.9	s26487.1	s43948.1	-
	14	CW	**	2.77	14.8-15.2	OAR14_15219778.1	OAR14_15268863.1	GTP2, DNAJA2, NETO2, ITFG1
	14	CW	**	2.85	17.5-18.3	s22016.1	s31525.1	ZNF423, CNEP1R1, HEATR3, TENT4B,
								ADCY7, BRD7
	16	CW	**	2.38	4.5-4.8	OAR16_4698169.1	OAR16_4737110_X.1	DUSP1, ERGIC1, RPL26L1, ATP6V0E1
	17	CW	**	2.7	26.7-26.9	OAR17_29456634.1	s52728.1	-
	17	GW	***	2.71	40.0-40.2	OAR17_43310387.1	OAR17_43356226.1	PPID, ETFDH, C4orf46, RXFP1
	17	GW	*	2.8	55.4-55.6	OAR17_60522910.1	s57641.1	CCDC60
	18	CW	**	1.68	6.6-6.8	OAR18_6256173.1	s37275.1	TTC23, SYNM
	21	CW	**	0.80	3.7-4.1	OAR21_4678890.1	s05427.1	-
	23	CW	**	0.68	8.3-8.5	OAR23_9159961.1	OAR23_9204273.1	-
	23	CW	**	2.7	34.4-34.6	OAR23_36433223.1	s07796.1	GATA6
	24	CW	*	2.5	7.4-7.7	OAR24_8719574.1	s03253.1	TMEM114, METTL22, ABAT,
								TMEM186, PMM2, CARHSP1
	24	CW	**	3.2	23.6-24.0	s12637.1	OAR24_26210500.1	-	KDM8
	24	CW	**	1.08	32.6-33.0	s64882.1	s17643.1	RCC1L, NCF1, GTF2IRD1, CLIP2
	26	CW	**	2.0	33.4-33.6	OAR26_38152631.1	OAR26_38251830.1	ADAM18, IDO1, IDO2	KAT6A

Frontiers in Genetics | www.frontiersin.org 2 December 2020 | Volume 11 | Article 616960

Frontiers in Genetics | www.frontiersin.org 3 December 2020 | Volume 11 | Article 616960

http://www.ensembl.org

http://crdd.osdd.net/raghava/dbem

Frontiers in Genetics | www.frontiersin.org

December 2020 | Volume 11 | Article 616960

ACKNOWLEDGMENTS

We are indebted to Paul Autran for the management of the experimental farm, Didier Foulquié and Fabien Carrière for managing the experimental flock, for animal care as well as for their active role in collecting experimental data. We also thank Alain Boissy and Eric Delval (INRAE UMR Herbivores) for their involvement in behavioral experiments, and Dominique François for his involvement in the management of animal breeding.

FUNDING

This study was granted by the INRAE metaprogram SELGEN (GENEPIGEN project), particularly for the phenotyping of GDMR. Genotyping was funded by ANR and APIS-GENE organizations in the frame of SheepSNPQTL project (ANR-08-GENM-039). The INRAE Division of Animal Genetics supported computational facilities and experimental costs for the production of the animals and the collection of behavioral, zootechnical and carcass traits.

DATA AVAILABILITY STATEMENT

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

ETHICS STATEMENT

The animal study was reviewed and approved by the Local ethics committee SCIENCE ET SANTE ANIMALES N • 115 (Ecole Nationale Vétérinaire de Toulouse).

AUTHOR CONTRIBUTIONS

DH, LD, CM-R, and SF conceived and designed the genetic experimental design and were involved in planning the study. FP-P was in charge of the phenotyping of GDMR. DH performed statistical and genetic analyses. DH, LD, SF, FP-P, and CM-R analyzed and interpreted the results. DH drafted the manuscript. LD, SF, and CM-R helped in writing, reviewing, and editing the final manuscript. All authors have read and approved the final manuscript.

Conflict of Interest:

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.


SAMSEL Pierre Guyomarc'h
email: pierre.guyomarch@u-bordeaux.fr

Mathilde Samsel

Patrice Courtaud

Pascal Mora

Bruno Dutailly

Sébastien

DUTAILLY Mora Pascal

Villotte Sébastien
email: sebastien.villotte@u-bordeaux.fr

Sébastien Villotte

Pierre Guyomarc

New data on the paleobiology of the Gravettian individual L2A from Cussac cave (Dordogne, France) through a virtual approach

Keywords: Virtual Anthropology, Photogrammetry, Sexual dimorphism, Probabilistic sex diagnosis, Coxal bone, Skull, Upper Paleolithic

come   L'archive ouverte pluridisciplinaire

Introduction

Since 2010, Cussac cave (Dordogne, France) underwent several extensive research campaigns that focused on several non-invasive archaeological and anthropological investigations [START_REF] Ferrier | Cussac Cave (Dordogne, France): The role of the rock support in the parietal art distribution, technical choices, and intentional and unintentional marks on the cave walls[END_REF][START_REF] Henry-Gambier | Grotte de Cussac (Le Buisson-de-Cadouin, Dordogne) : un exemple de comportement original pour le Gravettien[END_REF], Ledoux, et al., 2016[START_REF] Villotte | In situ study of the Gravettian individual from Cussac cave, locus 2 (Dordogne, France)[END_REF]. The cave hosts a unique combination of parietal art and human remains dating to the Gravettian (31,000 -22,000 BP). The striking artistic elements consist of more than 150 engravings, mostly animal and human depictions, which represents one the richest corpora of European parietal art for the Gravettian [START_REF] Aujoulat | La grotte ornée de Cussac (Dordogne). Observations liminaires[END_REF], Jaubert, et al., 2016). Similarly, human skeletal remains from this site significantly enriched the fossil record of the Upper Paleolithic of Europe, a period for which well-preserved human skeletons are relatively rare [START_REF] Henry-Gambier | Les populations européennes du Paléolithique supérieur[END_REF]. Several hundred human remains (fragments or complete bones) are scattered on the surface of the cave in at least three different loci. One of the loci, Locus 2, consists of a well-represented skeleton -partially covered by claycontained in a bear nest. Given the absence of duplicate elements, all the bones are attributed to the same individual, subject L2A. The current arrangement of the bones in the locus follows a logical anatomical pattern (Figure 1), indicating that the individual was deposited in ventral decubitus [START_REF] Henry-Gambier | Grotte de Cussac (Le Buisson-de-Cadouin, Dordogne) : un exemple de comportement original pour le Gravettien[END_REF][START_REF] Villotte | In situ study of the Gravettian individual from Cussac cave, locus 2 (Dordogne, France)[END_REF]. The other loci contain commingled human remains from several individuals. The consistency between the artistic style of the parietal engravings and two 14 C dates suggests that the human activity in the cave spans between 29,000 and 28,000 cal BP [START_REF] Jaubert | The chronology of human and animal presence in the decorated and sepulchral cave of Cussac (France)[END_REF]. Due to its unique characteristics, the Cussac cave is now protected under the national heritage status. The site is closed to the public, and even the scientific team (Programme Commun de Recherche, PCR Cussac) has a relatively short window of time for intervention each year due to the high concentration of CO 2 in the cave between spring and fall. Moreover, to ensure the protection of the preserved floor, no excavations nor samplings are planned in the near future. Given these limitations, the bioanthropological studies focus on non-invasive analyses, consisting, to date, of in situ observations and measurements of the visible elements. Using these methods, a preliminary study of L2A had been performed from a metallic walkway specially extended above the Locus 2 [START_REF] Villotte | In situ study of the Gravettian individual from Cussac cave, locus 2 (Dordogne, France)[END_REF]. Before and after this study, the locus was virtually documented via 3D photogrammetry, allowing for additional research on these human remains. Based on the characteristics of the auricular surface [START_REF] Schmitt | Une nouvelle méthode pour estimer l'âge au décès des adultes à partir de la surface sacro-pelvienne iliaque[END_REF], the first evaluation of L2A estimated an age-at-death between 20 to 49 years [START_REF] Villotte | In situ study of the Gravettian individual from Cussac cave, locus 2 (Dordogne, France)[END_REF]. The sex determination was done through visual scoring of the visible morphology of the os coxae [START_REF] Brůžek | A method for visual determination of sex, using the human hip bone[END_REF], and from measurements computed in a probabilistic tool for sexual diagnosis, DSP [START_REF] Murail | DSP: a tool for probalistic sex diagnosis using worldwide variation in hip-bone measurements[END_REF]. These morphometric methods applied to the left os coxae gave a male diagnosis. However, the presence of shared sexual traits in the morphoscopic evaluation, as well as a very low stature for this individual (outside the Gravettian range) called for further investigation into the sex of L2A [START_REF] Villotte | In situ study of the Gravettian individual from Cussac cave, locus 2 (Dordogne, France)[END_REF]. Moreover, the body proportions of this individual raised the question of a possible developmental anomaly [START_REF] Villotte | In situ study of the Gravettian individual from Cussac cave, locus 2 (Dordogne, France)[END_REF], an occurrence which appears to be abnormally frequent in the late Pleistocene sample [START_REF] Wu | An Enlarged Parietal Foramen in the Late Archaic Xujiayao 11 Neurocranium from Northern China, and Rare Anomalies among Pleistocene Homo[END_REF]. In addition, L2A was laid in a unique context for an Upper Paleolithic burial (a bear nest), and in ventral decubitus, a very uncommon position for the Gravettian period [START_REF] Henry-Gambier | Comportement des populations d'Europe au Gravettien : pratiques funéraires et interprétations[END_REF]. A new analysis of this individual will therefore contribute to the debate on possible specific burial treatments of pathological individuals during the Upper Paleolithic [START_REF] Formicola | From the Sunghir Children to the Romito Dwarf: Aspects of the Upper Paleolithic Funerary Landscape[END_REF][START_REF] Formicola | Bilateral absence of the lesser trochanter in a Late Epigravettian skeleton from Arene Candide (Italy)[END_REF][START_REF] Formicola | The Upper Paleolithic triple burial of Dolní Vĕstonice: pathology and funerary behavior[END_REF][START_REF] Mallegni | The human skeletal remains from the upper palaeolithic burials found in Romito cave (Papasidero, Cosenza, Italy)[END_REF][START_REF] Pettitt | The Palaeolithic origins of human burial[END_REF], Sparacello, et al., n.d., Villotte, et al., 2017). To this end, we employed a high-resolution photogrammetric tridimensional (3D) reconstruction to obtain virtual models of L2A's bones, and applied the now commonly used methods of virtual anthropology [START_REF] Weber | Virtual anthropology: a guide to a new interdisciplinary field[END_REF]. We present in this report a new analysis of the L2A pelvis with the DSP method, as well as an analysis of its cranial morphology in comparison with a sample from the Upper Paleolithic.

Material and Methods

The 3D photogrammetric model

Several photogrammetric acquisitions of the Cussac cave areas that include human remains were performed between 2012 and 2015, with global views and close-ups of the visible osseous elements. This allows for a 3D modelling and isolation of the different bones; the resulting 3D model of Locus 2 is displayed in Figure 1. More than 700 photographs were processed with a Nikon d700 (zoom 24-70 2.8 used at 24 mm) during three different campaigns using Photoscan (Agisoft©) and PMVS [START_REF] Furukawa | Accurate, Dense, and Robust Multi-View Stereopsis[END_REF] software. This resulted in several 3D point clouds with a spatial resolution ranging from 1 to 0.5 mm. Points clouds were meshed and scaled using local physical scales and topographic points (with x, y, z coordinates, acquired by a Leica© total station).

Methodology for the ossa coxarum analysis

The 3D photogrammetric models of the ossa coxarum of L2A were extracted from the global model to allow more accurate measurements (Figure 2). In order to reveal key anatomical areas for age and sex estimation, a layer of clay was removed from the left os coxae during one of the campaigns [START_REF] Villotte | In situ study of the Gravettian individual from Cussac cave, locus 2 (Dordogne, France)[END_REF]. Figure 2 displays the 3D models of the bones after cleaning (i.e. posterior to the record in Figure 1). The ossa coxarum were isolated with the Geomagic® Wrap 2014 software. One of the most objective and reliable morphometric methods to assess the sex from the pelvic region is DSP [START_REF] Murail | DSP: a tool for probalistic sex diagnosis using worldwide variation in hip-bone measurements[END_REF]. DSP has been used in the preliminary in situ study of Cussac L2A, and resulted in a male assessment from direct measurements of the left os coxae [START_REF] Villotte | In situ study of the Gravettian individual from Cussac cave, locus 2 (Dordogne, France)[END_REF]. As this method was already validated for use in a virtual approach [START_REF] Chapman | Sex determination using the Probabilistic Sex Diagnosis (DSP: Diagnose Sexuelle Probabiliste) tool in a virtual environment[END_REF], measurements were recorded for this study on the virtual model of the ossa coxarum in TIVMI software [START_REF] Dutailly | TIVMI, Treatment and Increased Vision for Medical Imaging[END_REF] to confirm the results previously obtained with traditional measures [START_REF] Villotte | In situ study of the Gravettian individual from Cussac cave, locus 2 (Dordogne, France)[END_REF].

Methodology of analysis for the cranial morphology

The refined photogrammetric acquisition of the cranium (Figure 3) was also cropped in Geomagic® Wrap 2014 for virtual processing. A layer of clay deposit covers most of the human remains in Locus 2, with a variable thickness; the teeth and spheno-occipital suture were cleaned during the in situ analysis for age estimation purposes [START_REF] Villotte | In situ study of the Gravettian individual from Cussac cave, locus 2 (Dordogne, France)[END_REF]. The mandible, although present, was not included in the present study, as it would require an extensive reconstruction of the non-visible parts (this element lies on its inferior surface and is partly covered by clay). The cranium of L2A is lying on its right side, thus several anatomical parts are non-visible on the model (Figure 3a). Consequently, the available landmarks on the left side were mirrored in TIVMI, using the sagittal plane as a reference (passing through the landmarks nasion, prosthion and bregma), to obtain the corresponding right points and ensure the collection of sufficient craniometrics data (see Table 2). The linear measurements were extracted in TIVMI through the Segment3D plugin, or using projections on planes when necessary [START_REF] Guyomarc'h | Anthropological Facial Approximation in Three Dimensions (AFA3D): Computer-Assisted Estimation of the Facial Morphology Using Geometric Morphometrics[END_REF]. The landmarks that were covered by a layer of clay were manually corrected to account for this layer, or were not collected when the thickness of the sediment could not be evaluated. At this stage of the analysis, this correction can only be an estimation based on comparisons with the areas not showing clay deposits (i.e. 1 to 2 mm thickness on the left lateral vault).

From the list of variables defined by Martin [START_REF] Bräuer | Anthropologie, Band I: Wesen und Methoden de Anthropologie[END_REF], a subsample of measurements were selected following their availability on L2A and their representation in the comparative sample. Table 1 lists the codes, measurements, and corresponding landmarks of the 30 retained variables, along with 18 cranial indices.

Cranial comparative sample

Data on the comparative sample were mainly collected from the literature from specimens with relatively well-preserved crania, with a known context of discovery, and available direct dating if possible. All the selected individuals date to a period covering the Middle-Upper Paleolithic (MUP) and the Late Upper Paleolithic (LUP). Some of these remains have recently been tested genetically (see [START_REF] Fu | The genetic history of Ice Age Europe[END_REF] for V1, O1, DV13, DV14, B1; [START_REF] Tarsi | Ricostruzione della storia genetica per via materna delle comunità paleolitiche dei Balzi Rossi, delle Arene Candide e del Romito, e di quelle neolitiche ed eneolitiche di Samari e di Fontenoce di Recanati, La cultura del morire nelle società preistoriche e protostoriche italiane dal paleolitico all'età del rame, origines, progetti[END_REF] for AC5), and the sex of most subjects was identified from their ossa coxarum morphology (CM1, AP1, BT2, BC1, BG1, BG5, GE5, PG25, DV3, DV16, PV1, S1, R6, AC2, AC4, AC3, AC12, IA, L1, LB4, LP5, SG1, ST1, OK1, OK2, S1). Additional specimen are used in the general comparison analyses, but are not included in the sex-related test (CM2, CM3, DV15, P3, P4, P9, P10, S5, B2, B3, R4, R6, AC12, C, RR1, CB). [START_REF] Villotte | The paleobiology of two adults skeletons from Baousso da Torre (Bausu da Ture) (Liguria, Italy): Implications for Gravetiian lifestyle[END_REF]) [START_REF] Villotte | Révision de l'âge au décès et du sexe des sujets adultes gravettiens[END_REF] Romito R6 (F), R4*, R6* Italy (S) LUP [START_REF] Mallegni | The human skeletal remains from the upper palaeolithic burials found in Romito cave (Papasidero, Cosenza, Italy)[END_REF]) [START_REF] Villotte | Enthésopathies et activités des hommes préhistoriques -Recherche méthodologique et application aux fossiles européens du Paléolithique supérieur et du Mésolithique[END_REF] Arene Candide AC2 (M), AC3 (M), AC4 (M), AC5 (M), AC12 (M) Italy (N) LUP [START_REF] Paoli | Gli scheletri mesolitici della caverna delle Arene Candide (Liguria)[END_REF][START_REF] Sergi | Il giovane paleolitico della caverna delle Arene Candide[END_REF]) [START_REF] Tarsi | Ricostruzione della storia genetica per via materna delle comunità paleolitiche dei Balzi Rossi, delle Arene Candide e del Romito, e di quelle neolitiche ed eneolitiche di Samari e di Fontenoce di Recanati, La cultura del morire nelle società preistoriche e protostoriche italiane dal paleolitico all'età del rame, origines, progetti[END_REF][START_REF] Villotte | Enthésopathies et activités des hommes préhistoriques -Recherche méthodologique et application aux fossiles européens du Paléolithique supérieur et du Mésolithique[END_REF] Le [START_REF] Vercellotti | The Late Upper Paleolithic skeleton Villabruna 1 (Italy): a source of data on biology and behavior of a 14.000 year-old hunter[END_REF]) [START_REF] Fu | The genetic history of Ice Age Europe[END_REF][START_REF] Villotte | Enthésopathies et activités des hommes préhistoriques -Recherche méthodologique et application aux fossiles européens du Paléolithique supérieur et du Mésolithique[END_REF]) Oberkassel OK1 (M), OK2 (F) Allemagne LUP Personal observation [START_REF] Trinkaus | The appendicular skeletal remains of Oberkassel 1 and 2[END_REF] Table 2. Reference material: code, sex, localization, period, data source and sex source of the specimens. * indicates the specimen is considered as of undetermined sex.

To summarize, the MUP subsample includes 10 males, 6 females and 9 undetermined sex individuals; the LUP group is composed of 10 males, 5 females, and 6 undetermined sex individuals. The 15 subjects without assessed sex were not used in the sex-related tests.

Statistical analysis of the cranium

The morphometric variations in terms of size are related to sexual dimorphism, but the morphology also contains variations in shape that may be related to sex: applying a log correction to linear variables is a method commonly used to extract shape information (e.g. [START_REF] Churchill | The Upper Paleolithic population of Europe in an evolutionary perspective[END_REF]). We use size-adjusted measurements following the ratio for shape analysis [START_REF] Jungers | Shape, relative size, and size-adjustments in morphometrics[END_REF], as proposed by [START_REF] Darroch | Canonical and principal components of shape[END_REF] to obtain the log shape ratios of the linear measurements. Including the linear variables, the indices, and the log shape ratios, a total of 78 variables are used to assess the position of L2A among the Upper Paleolithic (UP). As some variables do not follow a normal distribution, the nonparametric test was used (Wilcoxon-Mann-Whitney) to assess the cranial differences between groups on the measurements, indices, and log shape ratios. Principal Component Analyses (PCA) allow for a synthetization of the variables and help with the interpretation of the positioning of the L2A specimen within the comparative samples. Instead of replacing the missing values with the mean of variables (the default option in most software), we use the missMDA R package to perform more specific estimations of these unknown values following the subject's pattern [START_REF] Josse | Handling missing values in exploratory multivariate data analysis methods[END_REF]. Additionally, probabilistic distances, or adjusted Z-scores [START_REF] Maureille | Les dents inférieures du néandertalien Regourdou 1 (site de Regourdou, commune de Montignac, Dordogne)[END_REF], calculated between the means of the UP, MUP, and LUP samples and the variables from L2A will provide a visual assessment of its position relative to the comparative sample. Non-parametric tests were performed in Statistica© (v. 7.1, StatSoft©, Tulsa, OK), and PCA with R (R Core Team, 2016).

Results

DSP results

The virtual analysis of the ossa coxarum allowed for the collection of additional data compared to the in situ preliminary study (one variable on the left os coxae and three new measurements on the right os coxae). The DSP results are displayed in Table 3. [START_REF] Villotte | In situ study of the Gravettian individual from Cussac cave, locus 2 (Dordogne, France)[END_REF], and the measurements collected on the 3D photogrammetric records in TIVMI; n/o = not observable, n/a = not applicable; PF = probability of female sex; PM = probability of male sex. See [START_REF] Murail | DSP: a tool for probalistic sex diagnosis using worldwide variation in hip-bone measurements[END_REF] for measurements abbreviations.

The negligible differences between the in situ and photogrammetric measurements confirm the validity of the virtual application of the DSP method, with an average difference in the measurements of 2.4 % (sd: 2.2 %, min: 0.01 %, max: 7.9 %). The result on the left os coxae is similar to the conclusions of [START_REF] Villotte | In situ study of the Gravettian individual from Cussac cave, locus 2 (Dordogne, France)[END_REF], but the DSP failed to attribute a sex using the dimensions of the right os coxae. Any combination of four available measurements provides a higher probability for the male sex, but never above the acceptable threshold of 0.95.

Craniometric results

The visible cranial morphology of L2A is not particularly indicative of its sex, with some masculine (no frontal or parietal eminences, square orbits) and some feminine (V-shaped palate, smooth nuchal crests, small occipital condyles) traits displayed. The other commonly used features are either moderately expressed (mastoid process, supraorbital ridge) or not observable. A morphometric analysis was thus performed.

The 78 cranial variables were tested for significant differences between males and females with the Wilcoxon-Mann-Whitney test. Eleven variables show a significant sexual dimorphism in the UP sample (Supplementary Table 1), with males consistently displaying larger values than females. Overall, the values for L2A tend to indicate a relatively small size. This is illustrated by the PCA on the sexually dimorphic size variables (Figure 4). The first principal component (PC) of the size PCA (50.6%) allows for a distinctive separation between sexes, and L2A is within the male variability (positive values), but also close to the female group (negative values). Regarding the PCA on sexually dimorphic shape variables, the first PC (65.2%) shows a stronger overlap of the groups, and L2A is closer to the male variability (positive values) than the female one (negative values). A linear discriminant analysis computed on the sexually dimorphic variables (size and shape) classified L2A as male with a posterior probability of 0.97 (see Supplementary Data for details on the analysis). When considering the main cranial variables for which no sexual dimorphism was detected in the UP sample, L2A shows some very low (M17, M29, M52, M55) and very high (M44b, M54, M61) values, which are outside the range of variation known for the UP (Figure 5). L2A thus presents a very peculiar morphology, with a relatively flattened cranium, a wide face, a wide and short nasal aperture, a wide palate, and low orbital height (e.g. Figure 6). To illustrate this peculiar cranial shape, a selection of indices was processed in a PCA (I1, I4, I5, I12, I13, I18, I39, I42, I48, I54); because the results displayed a wide variability, the three first PCs are presented in Figure 7. On PC1 and PC3, L2A appears relatively separated from the UP sample, although cannot be considered a clear outlier.

Discussion and conclusions

The extensive photogrammetric record of Cussac cave allowed for the reconstruction of 3D models that could be used for new bioanthropological study, while minimizing the manipulation of those protected human remains. New metric data were generated for the Gravettian individual L2A, allowing for further evaluation of the cranial and pelvic remains through a virtual anthropology approach [START_REF] Weber | Virtual anthropology: a guide to a new interdisciplinary field[END_REF]. Although this kind of virtual analysis cannot fully replace the traditional study of the remains (i.e. we have no access to the whole surface of the bones or their internal structures), our research confirmed a male sex diagnosis for L2A. Despite his relatively small-sized cranium, the presence of some female non-metric cranial features, and a short stature relative to other Upper Palaeolithic individuals, the pelvis clearly indicates a male sex. Several Upper Palaeolithic fossils were traditionally diagnosed solely from the cranial morphology [START_REF] Gambier | Révision du sexe et de l'âge au décès des fossiles de Cro-Magnon (Dordogne, France) à partir de l'os coxal[END_REF][START_REF] Villotte | Révision de l'âge au décès et du sexe des sujets adultes gravettiens[END_REF], and following this method L2A may have been wrongly sexed as a female. The example of L2A highlights the risk of obtaining the wrong sex diagnosis using non-pelvic traits.

L2A's stature is well below the range of variation seen for Gravettian males, and its body proportions are uncommon for the Middle Upper Palaolithic, with the femur being rather long compared to the humerus and tibia [START_REF] Villotte | In situ study of the Gravettian individual from Cussac cave, locus 2 (Dordogne, France)[END_REF]. The peculiarity of this individuals is confirmed by the results of the present analysis, which show an unusual cranial morphology. L2A body proportions and cranial morphological peculiarities may be related to a growth delay, or a pathology. Congenital disorders are common in the Late Pleistocene human fossil record, and this has been interpreted as a possible indicator of high levels of inbreeding (Sparacello, et al., n.d., Wu, et al., 2013). L2A may be another case in support of this hypothesis. No definitive conclusion can be drawn without excavating the remains; however current evidence on cranial morphology and body proportions suggests that L2A can be considered as an outlier among the Gravettian variability.

The joint analysis of the biological profile of buried individuals and of their depositional context is essential for a proper interpretation of Paleolithic funerary behaviors (Mittnik, et al., 2016, Trinkaus and[START_REF] Trinkaus | Early modern human evolution in Central Europe -The people of Dolní Vĕstonice and Pavlov[END_REF]). L2A's biological characteristics, in the specific context of Cussac, raise the question of the social status of this male individual. L2A is one of the few Upper Paleolithic subjects from the southwest of France discovered in relative anatomical connection [START_REF] Henry-Gambier | Les populations européennes du Paléolithique supérieur[END_REF], 2008), and although several flooding episodes are attested to the layer of clay covering most of its bones, their position strongly advocates for a primary deposit in ventral decubitus [START_REF] Jaubert | The chronology of human and animal presence in the decorated and sepulchral cave of Cussac (France)[END_REF], a very rare position for the Gravettian. Indeed, out of ca. 33 attested Gravettian burials (with ca. 100 individuals), only three other subjects are in a ventral position (BT3, GE5 in the Grotte des Enfants double burial, and DV14 in the triple burial of Dolní Vĕstonice) (Henry- Gambier, 2008, Riel-Salvatore and[START_REF] Riel-Salvatore | Upper Palaeolithic mortuary practices in Eurasia: a critical look at the burial record[END_REF]. Moreover, very few Paleolithic human remains were found in a decorated cave or shelter [START_REF] Bartolomei | Una sepoltura epigravettiana nel deposito pleistocenico del Riparo Tagliente in Valpantena (Verona)[END_REF][START_REF] Foucher | Découverte de la mandibule d'un jeune enfant dans un niveau gravettien de la grotte de Gargas[END_REF][START_REF] Henry-Gambier | New hominid remains associated with Gravettian parietal art (Les Garennes, Vilhonneur, France)[END_REF][START_REF] Mallegni | The human skeletal remains from the upper palaeolithic burials found in Romito cave (Papasidero, Cosenza, Italy)[END_REF][START_REF] Mussi | Parietal art discovered at Arene Candide Cave (Liguria, Italy)[END_REF]. The body position of L2A, the fact that the body was deposited in a decorated cave, and its peculiar morphology echo with other UP burials displaying an association between a specific funerary treatment and marked osseous modifications, i.e. the burials of Dolní Vĕstonice, Sunghir, Arene Candide, and Romito sites [START_REF] Formicola | From the Sunghir Children to the Romito Dwarf: Aspects of the Upper Paleolithic Funerary Landscape[END_REF][START_REF] Formicola | The Upper Paleolithic triple burial of Dolní Vĕstonice: pathology and funerary behavior[END_REF][START_REF] Mallegni | The human skeletal remains from the upper palaeolithic burials found in Romito cave (Papasidero, Cosenza, Italy)[END_REF], Sparacello, et al., n.d., Trinkaus, et al., 2014).

The continuation of the biological analyses of the Cussac human remains -in particular the infracranial skeleton of L2A in the event of new excavation of the Locus 2 in the near futureas well as their integration with data from other scientific disciplines, will deepen our understanding of the Gravettian funerary practices in the southwest of France. The analysis (and reanalysis) of the other skeletal remains from this region will also be necessary. A laboratory study of these remains is essential to bridge the gaps in our knowledge about the biological variability during the Gravettian, and its possible relations with funerary behavior, which is still poorly documented.

Figure 1 .

 1 Figure 1. Visualization of the global 3D photogrammetric record of the Locus 2 with textures(in MeshLab©, ISTI, v.1.3.4)

Figure 2 .

 2 Figure 2. Visualization of the isolated 3D models of the ossa coxarum of L2A with textures (in MeshLab©, ISTI, v.1.3.4): (a) right os coxae in lateral view; (b) left os coxae in medial view.

Figure 3 .

 3 Figure 3. Visualization of the 3D model of L2A cranium with textures (in MeshLab©, ISTI, v.1.3.4): (a) isolated cranium in frontal view; (b) cranium in left lateral view; (c) cranium in inferior view.

Figure 4 .

 4 Figure 4. PCAs on sexually dimorphic size (a) and shape (b) variables (blue = males, red = females, undetermined sex = green). L2A is indicated with a star. Dimorphic variables used are M1, M8, M23, M26, M45, M48 and M51 for size morphology, and I2, I3, logM44b, logM46b and logM48 for shape morphology.

Figure 5 .

 5 Figure 5. Adjusted Z-scores between L2A, and the UP (green line), MUP (gray dashed line), and LUP (orange double-dashed line) comparative samples, computed on main cranial variables. The values -1 and 1 on the y axis represent 95% of the group variabilities; the position of the groups is relative to the L2A measurements.

Figure 6 .

 6 Figure 6. Bivariate plots for orbital (a) and nasal (b) height and breadth (in mm) in the UP sample (MUP = gray, LUP = black, male = square, female = circle, indeterminate = triangle, L2A = star).

Figure 7 .

 7 Figure 7. 3D scatterplot of PC1, PC2 and PC3 of the PCA performed on 10 cranial indices (MUP = blue, LUP = green, L2A = star).

Table 1 .

 1 Variables used in the study and raw data for L2A (in mm). * indicates the missing landmarks estimated for L2A.

	Code	Variable	Landmarks/measurements L2A	Code	Variable	Landmarks/measurements L2A
	(Bräuer,				(Bräuer,
	1988)				1988)
	M1	Maximum cranial length	Glabella -Opisthocranion	191.1 M54	Nasal breadth Apertion L -Apertion R	27.7
	M2	Glabella-inion length	Glabella -Inion	179.7 M55	Nasal height	Nasospinale -Nasion	42.3
	M3	Glabella-lambda length	Glabella -Lambda	179.9 M60	Maxillo-alveolar length	Prosthion -Alveolare	52.7
		Length of the skull	Basion -Nasion			Maxillo-	Ectomolare L -Ectomolare
	M5	base		98.6	M61	alveolar	R*	70.2
						breadth

  Table 2 lists the 46 subjects of the comparative sample with their source publications and the acronyms used in this study.

	Site	Specimen (Sex)	Localization	Period Data source	Sex diagnosis
	Cro-Magnon	CM1 (M), CM2*, CM3*	France (SW)	MUP	(Vallois and Billy, 1965)	(Gambier, et al., 2006)
	Abri Pataud	AP1 (F)	France (SW)	MUP	(Billy, 1975)	(Villotte, 2009)
	Dolní Vĕstonice and Pavlov	DV3 (F), DV13 (M), DV14 (M), DV15*, DV16 (M), PV1 (M)	Czech Republic	MUP	(Sládek, et al., 2000)	(Fu, et al., 2016, Trinkaus and Svoboda, 2006)
	Předmostí	P3*, P4*, P9*, P10*	Czech Republic	MUP	Brůžek, 2008) (Velemínská and	-
	Sunghir	S1 (M), S5*	Russia	MUP	(Trinkaus, et al., 2014) (Villotte, 2009)
	Ostuni	O1 (F)	Italy (S)	MUP	(Coppola, 2013)	(Villotte, 2009)
	Paglicci	PG25 (F)	Italy (S)	MUP	(Coppola, 2013)	(Villotte, 2009)
	Brno	B2 (I), B3 (I)	Czech Republic	MUP	(Coppola, 2013)	-
	Grotte des Enfants	GE5 (F)	Italy (N)	MUP	(Coppola, 2013)	(Villotte, et al., 2011)
	Barma del Caviglione	BC1 (F)	Italy (N)	MUP	(De Lumley, 2016, Guipert, et al., 2014)	(De Lumley, 2016)
	Barma Grande	BG1 (M), BG5 (M)	Italy (N)	MUP	(De Lumley, 2016, Graziosi, 1942)	(De Lumley, 2016, Villotte, et al., 2011)
	Baousso da Torre	BT2 (M)	Italy (N)	MUP

Table 3 .

 3 DSP input and results for the left and right ossa coxarum, including the measurements collected in situ by two observers

	Os Coxae	Observation	PUM SPU DCOX IIMT ISMM SCOX SS	SA SIS VEAC PF	PM	Sex
		Villotte et al. 2015, obs 1 n/o 28.7 214.0 37.0 n/o	n/o 75.6 76.4 n/o n/o 0.018 0.982 Male
	Left	Villotte et al. 2015, obs 2 n/o 30.3 216.0 38.2 n/o	n/o 75.1 78.2 n/o n/o 0.012 0.988 Male
		3D photogrammetry	n/o 29.4 216.9 38.4 n/o 152.1 75.1 76.0 n/o n/o 0.001 0.999 Male
		Villotte et al. 2015, obs 1 n/o n/o n/o n/o 100.0 n/o n/o n/o 38.8 47.6 n/a	n/a	n/a
	Right	Villotte et al. 2015, obs 2 n/o n/o n/o n/o 103.5 n/o n/o n/o 37.2 52.5 n/a	n/a	n/a
		3D photogrammetry	n/o n/o n/o 39.4 106.7 155.9 74.8 n/o 37.2 51.3 0.386 0.614 ND

Acknowledgements

Authors would like to acknowledge Vitale Stefano Sparacello for sharing additional data, providing insights on this study, and for a thorough editing of the paper, as well as Frédéric Santos for statistical support. We are grateful to the following institutions, curators and researchers that provided MS access to their facilities and granted her permission to examine the palaeoanthropological material in their care: Institut de Paléontologie Humaine (H. de Lumley), Museo Geologico G. G. Gemmellaro di Palermo (C. di Patti), Musée d'Archéologie Nationale (C. Schwab), Musée National de Préhistoire (J.-J.-Cleyet-Merle), Musée d'Archéologie Tricastine (M. Lert) and Service Régional d'Archéologie Auvergne-Rhône-Alpes (P. Chapuis), Musée d'Art et d'Archéologie du Périgord, Périgueux (V. Merlin-Anglade), Musée Labenche d'Art et d'Histoire, Brive (L. Michelin), Museum d'Histoire Naturelle Victor Brun, Montauban (A. Bergeret), R. W. Schmitz (LVR-LandesMuseum Bonn) and J. J. Hublin (Department of Human Evolution Max Planck Institute for Evolutionary Anthropology, Leipzig). Authors are also thankful for the Cussac Programme Commun de Recherche for their productive teamwork on this exceptional site.

This work was supported by the Agence Nationale de la Recherche (ANR) GRAVETT'OS (grant number ANR-15-CE-0004)

Supplementary Data

The following Supplementary Table 1 contains the descriptive statistics of the comparative sample (n=46) and the raw values for the subject L2A. 1. Descriptive statistics for the MUP and LUP groups and results of the Wilcoxon-Mann-Whitney test on linear distances, indices, and log shape ratios (log) with mean, standard deviations (sd), and sample size (n). The adjusted p level is reported with a significant level inferior to 0.05 (in bold). The raw values (mm for the linear variables) of Cussac L2A are reported for comparison.

Supplementary Table

The Supplementary Table 2 presents the specificities of a Linear Discriminant Analysis (LDA) based on dimorphic size and shape variables. Using this model, individuals with an undetermined sex where classified as male for P10, P3, C, and CM3; were classified as female RR1, R5, and B3; six subjects remained undetermined under the posterior probability of 0.90: CM2, P3, P4, B2, S5, and CB. These results indicate a trend in the morphology of this sub-sample of MUP individual, and are globally consistent with the sexual diagnosis based on the cranial remains of these individuals. However, due to the small sample size and the varying sexual dimorphism of the skull, such model cannot be considered per se as a reliable sex determination method.