Thierry Tribout

Pascal Croiseau

Rachel Lefebvre

Anne Barbat

Mekki Boussaha

Sébastien Fritz

Didier Boichard

Chris Hoze

Marie-Pierre Sanchez
email: marie-pierre.sanchez@inrae.fr

Confirmed effects of candidate variants for milk production, udder health, and udder morphology in dairy cattle

Background: Over the last years, genome-wide association studies (GWAS) based on imputed whole-genome sequences (WGS) have been used to detect quantitative trait loci (QTL) and highlight candidate genes for important traits. However, in general this approach does not allow to validate the effects of candidate mutations or determine if they are truly causative for the trait(s) in question. To address these questions, we applied a two-step, within-breed GWAS approach on 15 traits (5 linked with milk production, 2 with udder health, and 8 with udder morphology) in Montbéliarde (MON), Normande (NOR), and Holstein (HOL) cattle. We detected the most-promising candidate variants (CV) using imputed WGS of 2515 MON, 2203 NOR, and 6321 HOL bulls, and validated their effects in three younger populations of 23,926 MON, 9400 NOR, and 51,977 HOL cows.

Results: Bull sequence-based GWAS detected 84 QTL: 13, 10, and 30 for milk production traits; 3, 0, and 2 for somatic cell score (SCS); and 8, 2 and 16 for udder morphology traits, in MON, NOR, and HOL respectively. Five genomic regions with effects on milk production traits were shared among the three breeds whereas six (2 for production and 4 for udder morphology and health traits) had effects in two breeds. In 80 of these QTL, 855 CV were highlighted based on the significance of their effects and functional annotation. The subsequent GWAS on MON, NOR, and HOL cows validated 8, 9, and 23 QTL for production traits; 0, 0, and 1 for SCS; and 4, 1, and 8 for udder morphology traits, respectively. In 47 of the 54 confirmed QTL, the CV identified in bulls had more significant effects than single nucleotide polymorphisms (SNPs) from the standard 50K chip. The best CV for each validated QTL was located in a gene that was functionally related to production (36 QTL) or udder (9 QTL) traits.

Conclusions:

Using this two-step GWAS approach, we identified and validated 54 QTL that included CV mostly located within functional candidate genes and explained up to 6.3% (udder traits) and 37% (production traits) of the genetic variance of economically important dairy traits. These CV are now included in the chip used to evaluate French dairy cattle and can be integrated into routine genomic evaluation.

applications and outcomes of the "1000 Bull Genomes" project, Hayes and Daetwyler [START_REF] Hayes | 1000 bull genomes project to map simple and complex genetic traits in cattle: applications and outcomes[END_REF] noted that, even if the majority of polymorphisms within a cattle population can be tested using readily available whole-genome sequence data, the unambiguous identification of an individual mutation as causative for a complex trait remains the exception rather than the norm. GWAS on imputed WGS enables the targeting of small genomic regions such as genes, but the identification of causal polymorphisms is much less straightforward. Difficulties in pinpointing causal mutations arise from (i) the long-range linkage disequilibrium (LD) that exists in cattle breeds, which usually results in the detection of a set of variants in high LD rather than a single causal variant, (ii) variability in imputation accuracy, which may favor a variant in LD with the causal mutation rather than the mutation itself, and (iii) poor annotation of the bovine genome, in particular in regulatory regions, which makes it difficult to distinguish the best functional candidate in a set of variants.

Beyond providing a better understanding of the underlying biology of complex traits, the identification of causal mutations could be beneficial for genomic evaluation, especially across populations. The integration of causal mutations into genomic evaluation models could increase the accuracy of predictions and ensure the persistence of these models across generations or for distantly related individuals [START_REF] Liu | Improvement of genomic prediction by integrating additional single nucleotide polymorphisms selected from imputed whole genome sequencing data[END_REF]. Models that have been developed in major breeds might then be more easily transposed to smaller breeds, for which accurate genomic evaluation is difficult to implement. In addition, models with causal variants can account for interactions between genes more easily. However, to avoid the integration of false-positive candidate variants into models, their effects must first be validated in other populations that are as independent as possible. The Eurogenomics custom single nucleotide polymorphism (SNP) chip, which has been developed for bovine genomic selection, appears to be an ideal tool for this purpose. It contains an add-on feature that can be updated once or twice a year, and it is widely used in multiple breeds [START_REF] Boichard | Experience from large scale use of the EuroGenomics custom SNP chip in cattle[END_REF], which makes it possible to validate the effects of candidate variants detected by GWAS in different large populations.

In dairy cattle in particular, production traits are of major importance. First and foremost, high milk production is conditioned by a good and healthy udder. Mastitis is the most important health problem in dairy cattle and has an unfavorable genetic correlation with milk yield [START_REF] Oltenacu | The impact of genetic selection for increased milk yield on the welfare of dairy cows[END_REF][START_REF] Rupp | Genetic parameters for clinical mastitis, somatic cell score, production, udder type traits, and milking ease in first lactation Holsteins[END_REF]. Udder morphology is closely linked to sustainable milk production and is also associated with mastitis resistance [START_REF] Rupp | Genetic parameters for clinical mastitis, somatic cell score, production, udder type traits, and milking ease in first lactation Holsteins[END_REF] and longevity [START_REF] Larroque | Relationships between type and longevity in the Holstein breed[END_REF]. Thus, there are great benefits to considering all of these traits in the same study.

In order to disentangle the biological relationship between these complex traits and propose candidate causative variants, the objectives of this study were to identify genes, and the polymorphisms within them, that are responsible for the genetic variation in traits related to milk production, udder health, and udder morphology in the three main French dairy cattle breeds: Holstein (HOL), Montbéliarde (MON), and Normande (NOR). First, we conducted within-breed GWAS using imputed WGS of bulls with performances (Part I); then, we validated the effects of the candidate causal variants highlighted in the initial detection by performing withinbreed GWAS in statistically independent populations of cows (Part II).

Methods

This study comprised two parts. Part I consisted of identifying QTL and candidate variants from sequence-based GWAS of three bull populations. Part II aimed at confirming their effects by conducting a GWAS using the candidate variants from Part I and SNPs from the 50K SNP chip in three cow populations. For this study, we did not perform any experiments on animals; thus, no ethical approval was required. higher than 0.01 in at least one breed, and genotype frequencies had to be in Hardy-Weinberg equilibrium with P > 10 -4 .

In total, we analyzed 16 (HOL and MON) or 15 (NOR) routinely collected traits:

• Five milk production traits: milk yield (MY), protein yield (PY), fat yield (FY), protein content (PC), and fat content (FC); • Two udder health traits: somatic cell score (SCS) and clinical mastitis (CM). SCS was defined as SCS = 3 + log 2 (SCC/100,000) and averaged over monthly measures within lactation, with SCC being the number of somatic cells per ml of milk. CM was defined within lactation as a 0/1 trait with 1 corresponding to the occurrence of at least one clinical case before 150 days in milk; • Eight udder morphology traits, recorded by a type classifier during a classification visit: udder support (US), udder depth (UD), fore udder attachment (FUA), rear udder height (RUH), fore teat distance (FTD), udder balance (UB), and teat orientation (TO) in all breeds, teat length (TL) in MON and HOL. Scores, ranging from 1 to 9, were recorded only once per cow in first lactation; • Milking speed score (MSS), a subjective appraisal ranging from 1 to 5, given by the farmer and recorded with morphology traits.

In this paper, for convenience, health traits, type traits and milking speed are referred to as udder traits.

For all traits, the phenotypes used in the analyses were the daughter yield deviations (DYD) of each bull, defined as the average value of daughters' performance, adjusted for fixed and non-genetic random effects and for the breeding value of their dams [START_REF] Vanraden | Derivation, calculation, and use of national animal-model information[END_REF]. DYD are produced by the French national genetic evaluation systems of HOL, MON, and NOR populations with the models described at https ://inter bull.org/ib/gefor ms [11]. Mean reliabilities for all traits, excluding CM, ranged from 0.74 to 0.94, depending on the breed and on the trait (Table 2). Mean reliabilities for CM were lower (0.40 for NOR, 0.43 for MON and HOL), which was a result of both the lower heritability (about 0.02) of this trait and the fact that it began being recorded on farms more recently than other traits.

Imputation to whole-genome sequences

Using the UMD3.1 assembly, genotypes of all bulls were imputed to WGS with the FImpute software, which accurately and quickly processes large datasets [START_REF] Sargolzaei | A new approach for efficient genotype imputation using information from relatives[END_REF]. A two-step process was performed in order to improve imputation accuracy: from 50K to HD, and then from HD to WGS [START_REF] Binsbergen | Accuracy of imputation to whole-genome sequence data in Holstein Friesian cattle[END_REF]. All imputations were performed separately for each breed using either a breed-specific (from 50K to HD SNPs) or a multi-breed (from HD SNPs to WGS) reference panel depending on the targeted density [START_REF] Bouwman | Consequences of splitting whole-genome sequencing effort over multiple breeds on imputation accuracy[END_REF]. In each breed, imputations to the HD SNP level were performed using a within-breed reference population that included, respectively, 522 MON, 546 NOR, and 776 HOL bulls that had been genotyped with the Illumina BovineHD BeadChip (Illumina Inc., San Diego, CA). WGS variants were imputed from HD SNP genotypes using WGS variants of the 1147 Bos taurus bulls from Run4 of the 1000 Bull Genomes Project [START_REF] Bouwman | Meta-analysis of genome-wide association studies for cattle stature identifies common genes that regulate body size in mammals[END_REF]; these bulls represented 27 cattle breeds, and included 288 HOL, 28 MON, and 24 NOR individuals. WGS variants were selected by applying the protocol defined by the 1000 Bull Genomes consortium [START_REF] Bouwman | Meta-analysis of genome-wide association studies for cattle stature identifies common genes that regulate body size in mammals[END_REF][START_REF] Daetwyler | Whole-genome sequencing of 234 bulls facilitates mapping of monogenic and complex traits in cattle[END_REF]. First, short reads were filtered for quality and aligned to the UMD3.1 reference sequence [START_REF] Boussaha | Construction of a large collection of small genome variations in French dairy and beef breeds using whole-genome sequences[END_REF], and small genomic variations (SNPs and InDels) were detected using SAMtools 0.0.18 [START_REF] Li | The sequence alignment/map format and SAMtools[END_REF]. Raw variants were then filtered as described in Boussaha et al. [START_REF] Boussaha | Construction of a large collection of small genome variations in French dairy and beef breeds using whole-genome sequences[END_REF] to produce a dataset of 26,738,438 autosomal variants. Finally, filtered variants were annotated using the Ensembl variant effect predictor pipeline v81 [START_REF] Mclaren | The Ensembl Variant Effect Predictor[END_REF], and the effects of amino-acid changes were predicted using the SIFT tool [START_REF] Kumar | Predicting the effects of coding non-synonymous variants on protein function using the SIFT algorithm[END_REF]. Imputation accuracies were estimated in the MON and HOL datasets by calculating genotypic concordance rates; these values reached 0.90 and 0.94, respectively [START_REF] Sanchez | Within-breed and multi-breed GWAS on imputed whole-genome sequence variants reveal candidate mutations affecting milk protein composition in dairy cattle[END_REF]. Although the number of sequenced bulls was slightly lower in NOR than in MON, they contributed a higher proportion of the genes of the population and we assumed that imputation accuracy was similar in both breeds. Only variants with a MAF higher than 0.1% were retained for within-breed association analyses, i.e. around 12 million variants in each breed.

Table 2 Number of bulls with genotypes and phenotypes (DYD) and average reliability of their phenotypes for each trait, in Montbéliarde (MON), Normande (NOR), and Holstein (HOL) cattle

Whole-genome sequence association analyses

We performed within-breed and single-trait association analyses between all 12 million polymorphic variants (MAF ≥ 0.001) and the traits described in Table 2. All association analyses were performed using the mlma option of GCTA software (version 1.24), which applies a mixed linear model that includes the variant to be tested [START_REF] Yang | GCTA: a tool for genome-wide complex trait analysis[END_REF]:

where y is the vector of DYD standardized by the genetic standard deviation of the trait in the considered breed ( σ u_pop ); µ is the overall mean; b is the additive fixed effect of the variant to be tested for association; x is the vector of imputed genotypes, coded 0, 1, or 2 (number of copies of the second allele); u ∼ N 0, Gσ 2 u is the vector of random polygenic effects, with G the genomic rela- tionship matrix (GRM) calculated using the HD SNP genotypes, and σ 2 u the polygenic variance, estimated based on the null model y = 1µ + u + e and then fixed while testing for the association between each variant and the trait of interest; and e ∼ N 0, Iσ 2 e is the vector of random residual effects, with I the identity matrix and

(1)

y = 1µ + xb + u + e σ 2
e the residual variance. Because the variability of DYD reliability was limited, residuals were assumed to have a homogeneous variance.

In order to account for multiple testing, the Bonferroni correction was applied to the thresholds by considering 8 million independent tests, after pruning for complete linkage disequilibrium. Therefore, the 5% genome-wide threshold of significance corresponded to a nominal P-value of 6.3 × 10 -9 (-log 10 (P) = 8.2). When a given trait was significantly affected by multiple variants, the variants that were located less than 1 million base-pairs (Mbp) apart were grouped together. The bounds of the confidence intervals (CI) of each region were then determined by considering the positions of variants that were included in the upper third of the peak (individual CI). For a given trait in a given breed, CI that overlapped or were less than 1 Mbp away from each other were grouped in a QTL region. For each QTL, we then defined two CI: (1) a TOP-CI determined by the bounds of the individual CI in which we found the most significant results in the region and (2) an EXT-CI with bounds determined by the outermost positions after all overlapping individual CI were grouped. When only a single individual CI was present in a given region, TOP-CI and EXT-CI were identical. For each trait, the proportion of genetic variance explained by each QTL was estimated by σ 2 g_QTL = 2p ms (1p ms ) b2 ms , with p ms and bms the frequency and the estimated allelic substitution effect in genetic standard deviation units, respectively, of the variant with the most significant effect ( ms ) in the QTL region.

Selection of candidate variants from sequence-based GWAS results

Within each of the QTL regions detected in the sequence-based GWAS, we selected the most plausible variants (SNPs or small InDels) to explain the effects we observed. About 900 variants could be added on the custom part of the chip. Variant selection was performed within breed, trait and individual QTL. A similar number of variants was a priori allocated to each individual QTL. Consequently, due to the number of QTL finally detected, about 10 variants were selected for each individual QTL. Candidate variants with a MAF higher than 0.02 were chosen based first on the level of significance of their effect. For top variants with similar significance levels, the best candidates were discriminated based on their functional annotation with a priority for genic variants in coding (missense and loss of function) and regulatory regions. The selected variants, 855 in total, were then included on the custom part of version 6 of the Illumina EuroG10K BeadChip [START_REF] Boichard | Experience from large scale use of the EuroGenomics custom SNP chip in cattle[END_REF]. When these variants were InDels, their breakpoints were tested as done for SNPs, as described in Fig. 1 [START_REF] Boichard | Experience from large scale use of the EuroGenomics custom SNP chip in cattle[END_REF].

Part II: Validation of the effects of candidate causative variant in cows

The second part of this study was dedicated to validating the effects of these QTL regions and the candidate variants identified within them. To this end, we tested the effects of the candidate variants, as well as those of the 50K SNPs, on the performance of three statistically independent datasets of HOL, MON, and NOR cows.

Genotyping and imputations

Of all the cows genotyped for the purpose of genomic selection in France, we found 51,977 HOL, 23,926 MON, and 9400 NOR cows, born from 2014, whose production and udder phenotypes were not included in the DYD calculations of bulls used in Part I. Thus, phenotypes of bulls used in Part I and cows used in Part II were statistically independent. These cows were genotyped using the BovineSNP50 BeadChip (Illumina Inc.) or the customized low-density EuroG10K BeadChip (versions 1 to 5; Illumina Inc.). Missing genotypes were imputed with the FImpute software [START_REF] Sargolzaei | A new approach for efficient genotype imputation using information from relatives[END_REF] in a two-step procedure. Generic markers from the BovineSNP50 Beadchip were imputed using all 50K genotyped animals as the reference, as per the routine procedure of the French evaluation system. Then, customized markers were imputed using as a reference all males and females (with and without phenotypes) that had been genotyped using the EuroG10K BeadChip (versions 1 to 6), i.e. 52,630 HOL, 32,373 MON, and 12,316 NOR animals. After the imputation process, all cows with phenotypes had genotypes for the variants of both the 50K Beadchip and EuroG10K BeadChip version 6, including the candidate variants detected in Part I. The accuracy of imputation was assessed by calculating mean squared correlations (R 2 ) between imputed and true genotypes in a validation set of variants with MAF ≥ 1%; these values were equivalent in the three breeds and reached on average 97% for the 50K SNPs and 96% for the CV.

GWAS analyses

Single-trait association analyses were performed between all of the polymorphic variants of the 50K and EuroG10K Beadchips with MAF ≥ 1% (46,753, 44,832, and 44,659 SNPs in HOL, MON and NOR, respectively) and the 16 (HOL and MON) or 15 (NOR) traits described in Table 2. The phenotypes considered were the yield deviations (YD) of each cow, as estimated in the French national genetic evaluation programs of the HOL, MON, and NOR populations. YD can be interpreted as a cow's performance, adjusted for environmental effects; for traits with repeated measures, it is the weighted mean of the cow's performance, adjusted for non-genetic effects. As for bulls DYD, YD are by-products of the French evaluation system [11]. As in Part I, we used GCTA [START_REF] Yang | GCTA: a tool for genome-wide complex trait analysis[END_REF] and applied model (1) on the vector y of the YD of the cows, considering G , the genomic relationship matrix (GRM), calculated with the 50K SNP genotypes. The SNP effect was considered significant if its -log 10 (P) value was higher than 6 (5% genome-wide threshold after Bonferroni correction, i.e. 10 -6 ). As before, all variant positions were from the UMD3.1 assembly.

Results

Part I: Results from the bull sequence-based GWAS

GWAS of imputed whole-genome sequences of MON, NOR, and HOL bulls revealed 24, 12, and 48 QTL, respectively, with significant effects (-log 10 (P) ≥ 8.2; Table 3) on production (Fig. 2) or udder morphology and udder health traits (Figs. 3 and4). At least one QTL was identified for all traits in the HOL dataset with the exception of CM, but no QTL was found for five traits in MON (PY, CM, TL, FTD, and TO) and nine traits in NOR (PY, CM, SCS, UD, FUA, FTD, UB, TO, and MSS). For the three breeds, we detected a larger number of QTL linked with milk production (13, 10, and 30 in MON, NOR and HOL, respectively), than with udder morphology (8, 2, and 16, respectively) or udder health traits (3, 0, and 2, respectively). Each QTL explained from 1.1 to 11.1% of the genetic variance of its associated trait in MON, 1.7 to 18.4% in NOR, and 0.3 to 26.8% in HOL. In each of the three breeds, the largest number of QTL was found for PC (6, 5, and 11 in MON, NOR, and HOL, respectively; Fig. 2), and their cumulative effects explained 17.2% (MON), 20.0% (NOR), and 27.7% (HOL) of the genetic variance of this trait. In each breed, the QTL that explained the largest percentage of the genetic variance of a trait was associated with FC, and the cumulative effects Fig. 1 Design of the molecular test for a structural variant in a SNP chip (example of an insertion). Using the black arrow as a primer, the G allele reveals the insertion and the A allele the absence of insertion. A confirmation can be obtained with a second test on the other side of the insertion (primers = red arrows). Allele C reveals the insertion and T the absence of insertion of all the QTL detected for this trait accounted for 19.7%, 23.2%, and 37% of the total genetic variance in MON, NOR, and HOL, respectively. In the three breeds, both the number of QTL and their individual estimated effects were lower for udder traits than for production traits; consequently, the QTL that were identified for these traits explained a smaller part of their genetic variance. In addition, in contrast to the results for production traits, the udder morphology or health trait which had the largest percentage of genetic variance explained by QTL was different among breeds: SCS with 6.3% in MON, RUH with 5.2% in NOR, and UD with 6.3% in HOL.

As described in the "Methods" section, for each QTL we defined two confidence intervals (CI) using either the CI of the most significant individual QTL (TOP-CI) or by the inclusion of all individual CI of the QTL within the region (EXT-CI). Genomic annotations of the variants located in the 84 QTL regions (TOP-CI or EXT-CI) are summarized over all traits and breeds in Table 4. Considering all QTL together, 11,696 and 20,798 distinct variants with significant effects were located within the TOP-CI and EXT-CI regions, respectively. These variants were mainly located in intergenic regions (56.8 and 59.7% for TOP-CI and EXT-CI, respectively) or in introns of genes (28.1 and 29.2%, respectively) of the bovine genome. Only 50 (0.43%, TOP-CI) and 66 (0.32%, EXT-CI) of the variants were missense. The remaining variants were located in putative regulatory regions of the bovine genome: mainly, the upstream and downstream regions and, to a lesser extent, the 3′ UTR, 5′ UTR, and splicing regions of genes.

Within a given breed, QTL for multiple production traits or udder traits were sometimes located in the same genomic region. When we grouped QTL based on their location on the genome, the 84 QTL corresponded to 61 distinct regions, referred to as QTL ID in the first column of Tables 5 and6. Of these, 36 regions had effects on production traits and 25 had effects on udder morphology and/or health traits. With respect to production traits, the 36 distinct genomic regions corresponded to 53 QTL, which were located on Bos taurus (BTA) autosomes [START_REF] Qanbari | Classic selective sweeps revealed by massive sequencing in cattle[END_REF][START_REF] Hayes | 1000 bull genomes project to map simple and complex genetic traits in cattle: applications and outcomes[END_REF][START_REF] Liu | Improvement of genomic prediction by integrating additional single nucleotide polymorphisms selected from imputed whole genome sequencing data[END_REF][START_REF] Boichard | Experience from large scale use of the EuroGenomics custom SNP chip in cattle[END_REF]11,[START_REF] Bouwman | Consequences of splitting whole-genome sequencing effort over multiple breeds on imputation accuracy[END_REF][START_REF] Boussaha | Construction of a large collection of small genome variations in French dairy and beef breeds using whole-genome sequences[END_REF][START_REF] Li | The sequence alignment/map format and SAMtools[END_REF][START_REF] Sanchez | Within-breed and multi-breed GWAS on imputed whole-genome sequence variants reveal candidate mutations affecting milk protein composition in dairy cattle[END_REF][START_REF] Yang | GCTA: a tool for genome-wide complex trait analysis[END_REF][START_REF] Govignon-Gion | Multiple trait genetic evaluation of clinical mastitis in three dairy cattle breeds[END_REF][START_REF] Meredith | Genome-wide associations for milk production and somatic cell score in Holstein-Friesian cattle in Ireland[END_REF], and 29 (Table 5). The 25 regions (31 QTL) with effects on udder morphology and health were found on BTA1, 2, 4, 5, 6, 8, 9, 14, 17, 19, 20, 24, 26, 28, and 29 (Table 6). The largest numbers of QTL were found on BTA5, 6, 14, 20, and 19. In addition, 15 of the 61 genomic regions had effects on more than one trait (corresponding to 38 of the 84 QTL); however, no genomic region had pleiotropic effects on both production and udder morphology or health traits (i.e. there was no overlap between the CI of QTL for production traits and any of the other traits). Instead, there were 10 regions that each affected from two to five different production traits and five regions that affected two to three different udder morphology or health traits. Within a breed, even if more than one trait could be linked to a single region, the specific variants with the most significant effects on each trait were largely different. However,

Table 3 Number of QTL and total (TOT), lowest (Min), and largest (Max) percentages of genetic variance of the trait explained by the QTL detected in sequence-based GWAS performed on bulls in each breed

*For the description of the traits see there were some cases in which a single variant, always located within a gene, had significant effects on different traits, i.e. the variant with the most significant effects was the same for multiple traits. We identified 15 QTL ID, all linked with production traits, that were shared among the three breeds; they were located in five genomic regions and affected PC on BTA3 (at ~ 15 Mbp) and BTA6 (at ~ 87 Mbp, Fig. 2) and FC on BTA5 (at ~ 94 Mbp), BTA14 (at ~ 1.8 Mbp), and BTA27 (at ~ 36.2 Mbp). Six other regions (two for production traits and four for udder traits) had effects in two different breeds: in HOL and NOR on BTA5 (at ~ 118 Mbp for PC); in HOL and MON on BTA6 (at ~ 88.8 Mbp for UD, Fig. 4; at ~ 93 Mbp for PC), BTA19 (at ~ 60 Mbp for MSS and UD), and BTA24 (at ~ 34 Mbp for UB and RUH) and in MON and NOR on BTA17 (at ~ 62.7 Mbp for FUA and US). Although regions were shared among breeds, the variants with the most significant effect were different in each breed, with one exception: in one region located on BTA19, the intergenic variant rs109603247 had the most significant effect on MSS in both HOL and MON.

In all the QTL detected, the size of the TOP-CI ranged from 36.7 kb (BTA4 in HOL for FTD) to 1.9 Mb (BTA24 in MON for UB), with mean and median values being equal to 931 and 700 kb, respectively. TOP-CI contained from 0 to 31 genes (mean = 6.1; median = 3). As expected, EXT-CI were often broader, up to 12.8 Mb (mean = 2.1 Mb; median = 1.0 Mb) with a larger number of genes, up to 55 (mean = 8.5; median = 4). When we analyzed the EXT-CI of QTL, we observed that the majority contained at least a gene; only nine QTL (1 for production and 8 for udder morphology and/or health traits) were located entirely within intergenic regions. The variant with the most significant effect was located in an intergenic region for 19 out of 53 QTL identified for production traits and for 20 out of 31 QTL found for udder traits. All other variants presenting the most significant effects were located in intronic (28 for production traits and 9 for udder traits), upstream (4 for production traits and 1 for udder traits), downstream (2 for production traits and 1 for udder traits) or 5′UTR (1 for production trait) regions of genes. The genes in which these variants were located are indicated in Tables 5 and6.

Part II: Confirmation results on cows' performances

Within each of the 84 QTL detected in Part I, we selected the variants that best explained the observed results, hereafter named candidate variants, from sequencebased GWAS results from bulls. For technical reasons, a few of these candidate variants could not be included on the customized EuroG10K chip. In the end, one to 192 candidate variants from each of 80 of the 84 QTL (855 different variants in total) were added to the chip and tested for validation together with the standard 50K SNPs. As a consequence, even for the four QTL for which no candidate variant was added (two for production and two for udder traits), there were SNPs from the standard 50K chip that were located in the EXT-CI and were thus included in this confirmation study. We confirmed-i.e. found significant effects in the corresponding breed x trait analysis-the effects in cows of 54 out of the 84 QTL described in Tables 5 and6 (40 of 53 QTL for production traits, Table 7; 14 out of 31 QTL for udder traits, Table 8). In each of the validated QTL regions, we found significant effects (-log 10 (P) ≥ 6) for up to 99 candidate variants and up to 33 50K SNPs. Of the 80 QTL for which we tested candidate variants, the mean rank of the best candidate variant was 1.8 for all the QTL, for both production and udder traits, and 1.5 for the validated QTL (1.6 for production traits and 1.1 for udder traits). Thus, for the majority of the validated QTL, the variant with the most significant effect was one of the candidate variants selected in Part I for its level of significance and/or annotation; the exceptions were seven QTL that corresponded to four different genomic regions. Of these four regions, we found one in which only one candidate variant was present (at ~ 78 Mb on BTA4 for PC in NOR); another one in which the best candidate variant was ranked 2nd (at ~ 12 Mb on BTA5 for TL in HOL); and two regions linked with production traits in HOL, both located on BTA14 (~ 1.8 Mb and ~ 67.4 Mb). The first region on BTA14 (~ 1.8 Mb) had very significant effects on all five production traits, and for three of them (FC, FY, and PC), one of the candidate variants was ranked first in the peak. In contrast, for the second region (~ 68 Mb), the candidate variant with the most significant effects on FC, MY, and PC was ranked 3rd, 4th and 4th, respectively, in the peak, meaning that the top two or three variants were from the set of 50K SNPs. Therefore, for almost all the QTL for which the effects were validated in Part II of this study, candidate variants from Part I had more significant effects than the SNPs from the 50K chip.

Table 4 Genomic annotations of variants included within the confidence intervals (CI) of the 84 QTL, defined using either the CI of the most significant individual QTL (TOP-CI) or all the individual CI (EXT-CI) within each QTL region

In 47 of the validated QTL, a candidate variant from Part I presented the most significant effects; these corresponded to 39 unique variants. Six of these had the most significant effects on two or three different traits and/or different breeds. In particular, we identified three candidate variants having the most significant effects in different breeds: an intronic variant in the MGST1 gene (rs211210569) for FC in all three breeds; an intergenic variant on BTA6 (rs134776019) for PC in MON and NOR; and an intronic variant in the GC gene (rs436532576) for SCS or UD in HOL and MON (Fig. 5). Two additional variants presented the most significant effects on different traits within a single breed: the rs379230475 variant, located in the 5′UTR region of DGAT1 (BTA14), was the top variant for FC and MY in MON and the missense rs385640152 variant in GHR (BTA 20) was the top variant for PC and FC in HOL. In most of the genomic regions for which effects were observed in different breeds or traits, the variants that had the most significant effects were distinct, and multiple variants were located in the same gene in only a few cases (MGST1, DGAT1, and GPAT4). Of the 39 variants with the most significant effects, 10 were in intergenic regions (5 for production traits and 5 for udder traits) while 29 were located in genes listed in Tables 7 and8.

Discussion

The approach used in Part I of this study-GWAS on imputed whole genome sequences in bulls and the selection of candidate variants in QTL regions-led to the identification of 84 QTL for traits related to production [START_REF] Iso-Touru | Genome-wide association analysis of milk yield traits in Nordic Red Cattle using imputed whole genome sequence variants[END_REF], udder morphology [START_REF] Sodeland | Molecular characterization of a long range haplotype affecting protein yield and mastitis susceptibility in Norwegian Red cattle[END_REF], and udder health (5) in the three main French dairy breeds. In Part II, we investigated these QTL in statistically independent populations of cows, and confirmed the effects of 54 of them (40 of 53, 75%, for production traits and 14 of 31, 45%, for udder traits). In addition, by performing a GWAS with sequence-level resolution on thousands of bulls for which accurate phenotypes were available, we were able to propose 855 candidate causative variants in the QTL regions, of which 452 were validated in large populations of cows (9400 to 51,977, depending on the breed).

However, the number of QTL detected and validated differed between breeds. The sequence-based GWAS identified twice as many QTL in HOL [START_REF] Goddard | Genetics of complex traits: prediction of phenotype, identification of causal polymorphisms and genetic architecture[END_REF] as in MON [START_REF] Boichard | Genetic parameters for 1st lactation dairy traits in Friesian, Montbéliarde and Normande breeds[END_REF], and four times more than in NOR [START_REF] Sargolzaei | A new approach for efficient genotype imputation using information from relatives[END_REF]. Likewise, the proportion of validated QTL was also higher in HOL than in MON or NOR [START_REF] Pausch | A multi-trait metaanalysis with imputed sequence variants reveals twelve QTL for mammary gland morphology in Fleckvieh cattle[END_REF][START_REF] Sargolzaei | A new approach for efficient genotype imputation using information from relatives[END_REF], and 10 QTL, respectively). Furthermore, regardless of the breed, both the number of QTL and their level of significance varied among traits. In the sequence-based GWAS, the 84 QTL corresponded to 36 different genomic regions linked with production traits (meanlog 10 (P) value of 27.1), 23 regions associated with udder morphology (meanlog 10 (P) value of 11.2), and five regions for udder health traits (mean log 10 (P) value of 8.9).

Factors that can affect GWAS results

The differences that we observed in GWAS results may, at least in part, have been due to factors that were unique to the breeds, traits, and populations (bulls and cows) analyzed here.

Number of animals with phenotypes and genotypes

HOL, MON, and NOR cows represent 64, 19, and 9% of French dairy herds, respectively. For this reason, the number of animals with phenotypes was much larger in HOL than in MON or NOR (6262 HOL bulls vs 2434 MON and 2175 NOR for the primary detection; 51,977 HOL cows vs. 23,926 MON and 9400 NOR for the validation). This discrepancy clearly affected the power of detection in both sets of analyses: we were able to detect and validate QTL with smaller effects in the HOL population, and consequently, identified more QTL in total in HOL than in the other two breeds.

Imputation accuracies

The number of sequenced bulls included in RUN4 of the 1000 Bull Genomes Project, and therefore in the reference population for sequence-level imputation, were 288, 28, and 24 in HOL, MON, and NOR, respectively. Unsurprisingly, the estimated imputation accuracies were then higher in HOL than in MON [START_REF] Sanchez | Within-breed and multi-breed GWAS on imputed whole-genome sequence variants reveal candidate mutations affecting milk protein composition in dairy cattle[END_REF] and NOR. In addition, MON is related to the Simmental breed that is well represented in the 1000 bull genome population, whereas NOR is quite specific and likely benefits less from the sequences of the other breeds. In addition, we estimated that the 28 MON and 24 NOR bulls whose sequences were included in the reference population had cumulative contributions to the French populations of 64 and 59%, respectively. These differences may have also promoted a higher imputation accuracy in MON than in NOR and therefore explain the smaller number of QTL found for the NOR bulls.

Table 8 GWAS validation results for udder traits in

Between the bull and cow populations, missing genotypes were imputed based on different reference populations. For imputations of bull genotypes (WGS), we used a multi-breed reference population that consisted of their major ancestor bulls. This reference population was of limited size, especially within breed, which likely affected the accuracy of imputation, especially for breedspecific and/or low-MAF variants. For imputations of cow genotypes (50K SNPs + candidate variants), we used large within-breed reference populations that consisted of all animals genotyped with the EuroG10k chip; thus, imputation accuracy was much higher than that at the sequence level.

Heritability and reliability of traits

Differences among traits in the numbers of QTL detected in the sequence-based GWAS could also be explained by differences in DYD reliabilities. DYD is considered as a bull's own performance for a trait, the heritability of which would be equal to the reliability of the DYD value. The higher the reliability, the smaller the residual variance and the higher the detection power. In addition to the heritability of the trait, the reliability of the DYD also depends on the effective daughter contribution [START_REF] Fikse | Weighting factors of sire daughter information in international genetic evaluations[END_REF], and on average, progeny groups were a little larger in HOL than in MON and NOR. Because udder health traits had lower heritabilities (h 2 = 0.018 to 0.15), the reliability (REL) of their DYD values was lower (REL = 0.40 to 0.88) than for udder morphology traits (h 2 = 0.15 to 0.45; REL = 0.74 to 0.95) and production traits (h 2 = 0.30 to 0.50; REL = 0.89 to 0.95) (Tables 2 and Additional file 1: Table S1). In addition, morphological traits were recorded only once for each daughter, whereas DYD calculations for production and health traits included up to three lactations per cow. Finally, recording of CM started only recently and is not exhaustive [START_REF] Govignon-Gion | Multiple trait genetic evaluation of clinical mastitis in three dairy cattle breeds[END_REF], meaning that DYD information is available for fewer bulls, with smaller informative progeny groups. All these reasons explain why the power of detection decreased from analyses of milk composition to those of milk yield, udder type, SCS, and finally CM.

In the cow confirmation study, the sample size was larger than in the bull populations, but the reliability of the traits, equal to the heritability (for non-repeated records), was always lower than reliability of the DYD, and for CM, considerably so. Depending on the trait and the population in question, the power of detection in the cow populations was either higher (e.g., for HOL and MON and high or medium heritability traits) or lower (e.g., for CM). The resulting lower power of the validation #: number of a ID number associated with a group of QTL linked with udder traits that had overlapping confidence intervals or less than 1 Mbp distance between the bounds of the confidence intervals b Individual QTL; somatic cell score (SCS), clinical mastitis score (CM), udder support (US), udder depth (UD), fore udder attachment (FUA), rear udder height (RUH), fore teat distance (FTD), udder balance (UB), teat orientation (TO), teat length (TL), and milking speed score (MSS); c candidate variants selected from sequence-based GWAS results dataset in some cases could be a possible explanation why certain variant effects were unconfirmed. For these reasons, we were able to explain a higher percentage of genetic variance for the most heritable traits and to detect small QTL for the phenotypes with the highest accuracy. Regardless of the breed analyzed, our results varied widely among traits. We detected no QTL for CM, the trait with the lowest values of heritability (≤ 0.023) and reliability (≤ 0.43), while for FC and PC, which had the highest heritability (0.50) and reliability (0.92-0.94), we recovered the largest number of QTL (up to 11 for PC in HOL) which together explained the highest percentage of genetic variance of any trait (up to 37% for FC in HOL). Because significant effects are likely to be overestimated, it is possible that the percentage of variance explained by each QTL may have been artificially high. The number of detected QTL was rather limited. This is explained by the very conservative detection threshold used (P ≤ 6.10 -9 ; log 10 (P) ≥ 8.2) that decreased power of detection and excluded the QTL with smaller effects. For example, by decreasing the detection threshold tolog 10 (P) = 7 (P ≤ 6.10 -7 ), we identified two additional QTL for CM in MON and HOL. These were located in a single genomic region around 88.5 Mb on BTA6 in the region of the GC gene, where we had found significant effects on udder morphology traits and SCS (Fig. 5).

Table 8 (continued)

QTL confirmation rate

In spite of the application of a very strict detection threshold to the bull GWAS results, about one-third of these QTL were not found in the cow populations. Several explanations could explain this situation. First, it is important to note that nearly all the highly significant QTL and all the QTL present in several breeds or affecting several traits were confirmed. A few QTL with log 10 (P) > 10 were not confirmed but this was due to technical problems, the best selected variants being lost during the design of the chip. Most unconfirmed QTL, especially for udder conformation (10 out of 16), were detected in the bull population withlog 10 (P) values between 8.2 and 10 and hadlog 10 (P) < 6 in the cow populations. Two reasons may be advocated. For these QTL, annotations were frequently very poor and we may have selected inappropriate variants. This point is especially critical when a small number of variants was selected. Indeed, due to our selection strategy, the enrichment in candidate variants increased with significance level in the bull populations, QTL size, and QTL sharing across breeds and traits, and the smallest QTL received a small number of candidate variants. In addition, we cannot exclude that some results that were unconfirmed in the cow population represent false positives.

QTL shared among breeds and traits

Although our results may have been shaped by factors specific to the breeds, traits, and populations analyzed here, still we successfully identified and validated QTL which were shared among more than one breed or related trait.

QTL shared among breeds

Five QTL associated with milk production traits were shared among all three breeds, while six other QTL were found in two of the three breeds (2 QTL for production and 4 for udder traits). Most of the QTL found in two breeds were shared between HOL and MON (4 QTL), probably because these were the two breeds FOR which we found the largest number of QTL. As mentioned above, the very strict detection threshold applied for the bull GWAS excluded some potential variants that also mapped at the same location for the same trait in another breed; thus, this reduced the number of significant results shared between breeds. For example, in the CI of QTL ID 42 (BTA5), detected for UD in HOL, we found a variant at 88,862,824 bp that also had, for the same trait, a significant effect in MON cows (-log 10 (P) = 7.4, results not shown) and an effect close to significance in MON bulls (-log 10 (P) = 7.2).

With these results, we were able to validate QTL shared among breeds for certain traits of interest. However, as previously reported from other studies conducted in multiple breeds at the nucleotide-level resolution [START_REF] Sanchez | Within-breed and multi-breed GWAS on imputed whole-genome sequence variants reveal candidate mutations affecting milk protein composition in dairy cattle[END_REF][START_REF] Pausch | Meta-analysis of sequence-based association studies across three cattle breeds reveals 25 QTL for fat and protein percentages in milk at nucleotide resolution[END_REF], the variants with the most significant effects for a given trait differed largely among breeds. The reason for this result remains unclear. This could indicate that the causal mutations differed across breeds, but it may also be the result of differences in the quality of imputation of candidate variants among breeds. Within a QTL region, the effects of variants with the highest imputation accuracy, which are not necessarily the same across breeds, were probably estimated more accurately and were thus more likely to be significant. As shown later in Discussion, this hypothesis seems to be supported by the fact that, for several QTL detected in more than one breed, a shared variant often ranked highly among the best significant variants, even if it was not the very best. Precise identification of causal variants is further complicated by the presence of strong LD over large regions beyond the gene level.

QTL shared among traits

Within a breed, many of the QTL that we detected had effects on more than one production or udder trait (morphology and health). For example, in HOL, the QTL linked with production traits on BTA20 had effects on MY, PC, and FC, whereas the QTL identified on BTA29 for udder traits affected UD, UB, and SCS. These results are consistent with estimates of genetic correlations between milk yield and milk composition traits [START_REF] Boichard | Genetic parameters for 1st lactation dairy traits in Friesian, Montbéliarde and Normande breeds[END_REF] and between udder health and udder morphology traits [START_REF] Oltenacu | The impact of genetic selection for increased milk yield on the welfare of dairy cows[END_REF][START_REF] Rupp | Genetic parameters for clinical mastitis, somatic cell score, production, udder type traits, and milking ease in first lactation Holsteins[END_REF]. However, although significant genetic correlations were reported between milk production and udder morphology or udder health traits [START_REF] Rupp | Genetic parameters for clinical mastitis, somatic cell score, production, udder type traits, and milking ease in first lactation Holsteins[END_REF], we were not able to identify any QTL that had overlapping CI for production and udder (morphology or health) traits, even when we considered those with less-significant effects for CM (P ≤ 6•10 -7 ). As an example, the QTL found on BTA6, which had effects on both udder health and morphology in all three breeds (Fig. 5), was located in the vicinity of another QTL that was detected in all three breeds and had effects on PC or PY. Depending on the breed, the variants with the most significant effects on udder traits were located between 88.5 and 88.9 Mb, while those with the most significant effects on production traits were located between 87 and 87.6 Mb. This region of the bovine genome (86-90 Mb) has been the subject of particular interest over the last ten years for its effects on milk production and udder health (clinical mastitis and somatic cell scores) [START_REF] Olsen | Fine mapping of a QTL on bovine chromosome 6 using imputed full sequence data suggests a key role for the group-specific component (GC) gene in clinical mastitis and milk production[END_REF][START_REF] Sodeland | Molecular characterization of a long range haplotype affecting protein yield and mastitis susceptibility in Norwegian Red cattle[END_REF][START_REF] Meredith | Genome-wide associations for milk production and somatic cell score in Holstein-Friesian cattle in Ireland[END_REF][START_REF] Fang | Integrating sequence-based GWAS and RNA-Seq provides novel insights into the genetic basis of mastitis and milk production in dairy cattle[END_REF]; the two most recent studies, both performed at the nucleotide level, identified a single variant or distinct but very close variants with the most significant effects on both milk yields and mastitis resistance [START_REF] Sodeland | Molecular characterization of a long range haplotype affecting protein yield and mastitis susceptibility in Norwegian Red cattle[END_REF][START_REF] Fang | Integrating sequence-based GWAS and RNA-Seq provides novel insights into the genetic basis of mastitis and milk production in dairy cattle[END_REF]. Our study did not confirm the existence of a QTL with pleiotropic effects in this region; instead, our data suggest the presence of two neighboring QTL.

Further investigations of QTL regions reveal the best candidate genes and variants

For QTL that were shared between breeds, and that had effects on multiple traits or were identified in both bulls and cows, the results obtained at the nucleotide level appeared to be very sensitive to the accuracies of phenotypes and genotypes. In most cases, the variant with the most significant effects differed among traits, among breeds or among populations within a breed (bulls vs. cows). However, in most of these QTL regions, a detailed investigation of the GWAS results revealed the genes and the variants that are most likely to be causative.

Candidate genes and variants for udder traits

The QTL that were confirmed to have the most significant effects on udder traits were located on BTA5, 6, and 14. In these three regions, GWAS results pinpointed ABCC9, GC, and PLAG1 as the candidate genes, respectively.

The ABCC9 (ATP binding cassette subfamily C member 9) gene, located on BTA5, was associated here with UD and FUA in HOL bulls and cows. In both analyses, the variant(s) with the most significant effects for both traits were located in intronic regions of the ABCC9 gene, but approximately 12 to 24 kb apart: at 88,800,994 bp (rs110461240) in bulls and at 88,812,245 bp and 88,824,857 bp (rs209585944 and rs209893772, respectively, with the same significance level) in cows. This gene has previously been linked with milk production and fertility [START_REF] Pausch | Meta-analysis of sequence-based association studies across three cattle breeds reveals 25 QTL for fat and protein percentages in milk at nucleotide resolution[END_REF] and more recently with udder morphology (UD and FUA), milk production (MY and PY), and daughter pregnancy rate [START_REF] Jiang | Functional annotation and Bayesian fine-mapping reveals candidate genes for important agronomic traits in Holstein bulls[END_REF] in Holstein cattle. However, Jiang et al. [START_REF] Jiang | Functional annotation and Bayesian fine-mapping reveals candidate genes for important agronomic traits in Holstein bulls[END_REF], who performed a multi-trait analysis at the sequence level, failed to detect shared variants associated with different trait groups, suggesting the existence of several causal mutations for the different traits. In their study, the variants with the most significant effects were located at 88,818,703 bp (intron) for FUA and at 88,823,164 bp (splice region) for UD, i.e. between the most plausible candidate variants that we identified here. In our study, the best candidate variants identified by Jiang et al. [START_REF] Jiang | Functional annotation and Bayesian fine-mapping reveals candidate genes for important agronomic traits in Holstein bulls[END_REF] were confirmed to have very significant effects on UD (P = 4.3•10 -13 and 4.4•10 -13 , respectively) and FUA (P = 1.3•10 -11 and 1.4•10 -11 , respectively). In a nearby region, we also found significant effects on PY in HOL bulls but we could not confirm this in cows; the most significant variant in that case was the same as the one we detected for udder traits (rs136903701; 88,830,128 bp) but did not reach the level of significance (-log 10 (P) = 5.5). The ABCC9 gene encodes a protein involved in the formation of the ATP-sensitive potassium channels in different muscles. These channels are expressed in many tissues and regulate different cellular functions; thus, mutations in the ABCC9 gene could have potential effects on many traits.

As mentioned earlier, the region around 88.7 Mb on BTA6 has previously been linked with mastitis resistance [START_REF] Fang | Integrating sequence-based GWAS and RNA-Seq provides novel insights into the genetic basis of mastitis and milk production in dairy cattle[END_REF][START_REF] Cai | Prioritizing candidate genes post-GWAS using multiple sources of data for mastitis resistance in dairy cattle[END_REF][START_REF] Sahana | Genome-wide association study using high-density single nucleotide polymorphism arrays and whole-genome sequences for clinical mastitis traits in dairy cattle[END_REF] or udder morphology [START_REF] Pausch | A multi-trait metaanalysis with imputed sequence variants reveals twelve QTL for mammary gland morphology in Fleckvieh cattle[END_REF]. Here, we detected effects of this region on both type of traits-SCS in HOL and UD in HOL and MON-and, furthermore, it was the only region associated with mastitis resistance that we successfully validated in cows. Interestingly, the candidate variant that had the most significant effects (rs436532576; 88,723,742 bp) on these two traits in these two breeds in the validation GWAS was the most plausible causative variant previously identified in Red Danish [START_REF] Cai | Prioritizing candidate genes post-GWAS using multiple sources of data for mastitis resistance in dairy cattle[END_REF][START_REF] Sahana | Genome-wide association study using high-density single nucleotide polymorphism arrays and whole-genome sequences for clinical mastitis traits in dairy cattle[END_REF] and German Fleckvieh cattle [START_REF] Pausch | A multi-trait metaanalysis with imputed sequence variants reveals twelve QTL for mammary gland morphology in Fleckvieh cattle[END_REF]. The effect of this candidate variant did not reach the level of significance in NOR (P = 4.10 -3 ), but in NOR, the MAF of this variant was lower (0.21) than in MON (0.37) and HOL (0.40). The rs436532576 variant is located in an intronic region of GC (vitamin D binding protein), which was previously proposed as a candidate gene for resistance to mastitis in cattle because it encodes a Gc-globulin that is involved in both the transport of vitamin D to monocytes and phagocytic activity in macrophages [START_REF] Sahana | Genome-wide association study using high-density single nucleotide polymorphism arrays and whole-genome sequences for clinical mastitis traits in dairy cattle[END_REF].

On BTA14, the most plausible causative variants were identified in the PLAG1 (PLAG1 zinc finger) gene in both MON bulls and cows. In the GWAS performed on imputed WGS of bulls, the variant with the most significant effects on RUH was located in the 5′-UTR region of PLAG1 (rs210030313). Unfortunately, for technical reasons, it was not possible to add this variant to the customized chip. Instead, the variant with the most significant effects on this trait in MON cows was an intronic variant in PLAG1, located at 25,015,640 bp on BTA14 (rs109815800). The PLAG1 gene has been associated with stature in cattle [START_REF] Bouwman | Meta-analysis of genome-wide association studies for cattle stature identifies common genes that regulate body size in mammals[END_REF][START_REF] Karim | Variants modulating the expression of a chromosome domain encompassing PLAG1 influence bovine stature[END_REF] and humans [START_REF] Pryce | Polymorphic regions affecting human height also control stature in cattle[END_REF] but also with udder morphology [START_REF] Pausch | A multi-trait metaanalysis with imputed sequence variants reveals twelve QTL for mammary gland morphology in Fleckvieh cattle[END_REF]; variant rs109815800, which is a SNP on the Illumina Bovine HD BeadChip, was the most strongly associated of the whole-genome sequence variants with stature in the bovine meta-analysis of Bouwman et al. [START_REF] Bouwman | Meta-analysis of genome-wide association studies for cattle stature identifies common genes that regulate body size in mammals[END_REF] and with udder depth in the study of Pausch et al. [START_REF] Pausch | A multi-trait metaanalysis with imputed sequence variants reveals twelve QTL for mammary gland morphology in Fleckvieh cattle[END_REF]. However, in our study, this variant was ranked 3rd by the sequence-based GWAS, after the 5′-UTR variants located at 25,052,440 bp (1st) and 25,052,394 bp (2nd). These two variants are also plausible causal variants as they present a higher probability of being located within a transcription binding site. Moreover, Pausch et al. [START_REF] Pausch | A multi-trait metaanalysis with imputed sequence variants reveals twelve QTL for mammary gland morphology in Fleckvieh cattle[END_REF], who found no association when UD was conditioned on body height, suggested that the association between PLAG1 and udder morphology traits could be the result of phenotypic variation in body size rather than a true effect on mammary gland morphology. In our study, the lack of a significant effect of PLAG1 on other udder morphology traits than UD (less dependent on stature than UD) tends to support this hypothesis.

We identified and confirmed the effects on mammary gland morphology of other candidate variants on BTA5 in an intron of TMTC2 (transmembrane O-mannosyltransferase targeting cadherins 2), on BTA20 upstream of ISL1 (ISL LIM Homeobox 1), and on BTA26 in the 3′-UTR of RAB11FIP2 (RAB11 family interacting protein 2). TMTC2 was previously found to be associated with six udder type traits by Jiang et al. [START_REF] Jiang | Functional annotation and Bayesian fine-mapping reveals candidate genes for important agronomic traits in Holstein bulls[END_REF]. Instead, no such relationship has been reported for either ISL1, which encodes a member of the LIM/homeodomain family of transcription factors, or RAB11FIP2.

Candidate genes and variants for production traits

Among the genes that we identified here as being associated with milk production and composition traits, there are a number of well-characterized functional candidate genes: GHR, which encodes a growth hormone receptor, PAEP and CSN2, which encode milk proteins, and DGAT1, GPAT4 and FASN, all of which encode enzymes involved in the metabolism of fatty acids in milk. We also identified several other candidate genes with less well known functions or for which functional links with dairy traits have not yet been established: MGST1, CDDC57, TBC1D22A, VPS13B, PICALM, and GRAMD4.

The F279Y missense mutation in the GHR (growth hormone receptor) gene, which has previously been implicated in the genetic variation of PC and FC [START_REF] Viitala | The role of the bovine growth hormone receptor and prolactin receptor genes in milk, fat and protein production in Finnish Ayrshire dairy cattle[END_REF][START_REF] Blott | Molecular dissection of a quantitative trait locus: a phenylalanine-totyrosine substitution in the transmembrane domain of the bovine growth hormone receptor is associated with a major effect on milk yield and composition[END_REF], had the most significant effects on PC (P = 1.2•10 -102 ) and FC (P = 4.4•10 -55 ) in HOL cows, and was ranked 2nd for MY (P = 1.4•10 -10 ), confirming the QTL region identified in HOL bulls. The allele responsible for a decrease in the protein and fat contents of milk had a frequency of 0.12. In the GWAS performed on bulls, the F279Y variant had very significant effects on PC (P = 7.9•10 -23 ) but the variant with the most significant effects in this region was an intergenic variant located at 32,254,539 bp (P = 2.5•10 -30 ), i.e. relatively distant (~ 250 kb) from the causal mutation; this suggested poor imputation accuracy in the region surrounding GHR. No effects of this region were detected in NOR and MON cows, but the MAF of the F279Y variant was much lower in these two breeds (0.07 and 0.006, respectively). In contrast to Viitala et al. [START_REF] Viitala | The role of the bovine growth hormone receptor and prolactin receptor genes in milk, fat and protein production in Finnish Ayrshire dairy cattle[END_REF], we found no significant effects of the S18N variant in the PRLR (prolactin receptor) gene, located approximately 7 Mb downstream of GHR, on any of the milk production traits, although this missense mutation was polymorphic in MON, NOR, and HOL (MAF = 0.23, 0.42 and 0.16, respectively). Our result corroborates the hypothesis that the S18N mutation in PRLR may not be causative but is instead, at least in populations in which its effects have been demonstrated, in LD with the causal mutation [START_REF] Pausch | Short communication: validation of 4 candidate causative trait variants in 2 cattle breeds using targeted sequence imputation[END_REF].

In HOL, we identified and validated two QTL located near the PAEP (progestagen associated endometrial protein) gene, which encodes β-lactoglobulin (BTA11 at ~ 103.3 Mbp), and likewise confirmed the effects of the cluster of casein genes encoding the αs1 (CSNS1), αs2 (CSNS2), β (CSN2), and κ (CSN3) caseins (BTA6 at ~ 87.2 Mbp) in MON, NOR, and HOL. Although our results differed depending on the breed and population (bulls or cows) analyzed, PAEP and CSN2 were found to be the best candidate genes in HOL cows for the QTL acting on FC and PY, respectively. The best candidate variant in CSN2 was the missense variant responsible for the A1/B and A2 protein variants (at 87,181,619 bp; rs43703011), which has previously been implicated in milk composition and cheese-making quality [START_REF] Caroli | Invited review: milk protein polymorphisms in cattle: Effect on animal breeding and human nutrition[END_REF]. This variant also had very significant effects on PC in all three breeds (MON P = 8.8•10 -28 , MAF = 0.38; NOR P = 7.2•10 -13 , MAF = 0.28; and HOL P = 9.8•10 -11 , MAF = 0.33) but it was not ranked among the top 10 variants of the peak for this trait. We also detected and confirmed another QTL on BTA6 in HOL in the region of the ABCG2 gene previously identified for milk composition [START_REF] Cohen-Zinder | Identification of a missense mutation in the bovine ABCG2 gene with a major effect on the QTL on chromosome 6 affecting milk yield and composition in Holstein cattle[END_REF]. Only two of the 138 variants with significant effects on FC and/or PC in Holstein bulls, located in the EXT-CI of the QTL (37.5-38.5 Mb), were in the ABCG2 gene (rs136230937 at 38,015,146 bp and rs110063427 at 38,020,110 bp). They are intronic and therefore distinct from the rs43702337 missense variant (at 38,027,010 bp) described by Cohen-Zinder et al. [START_REF] Cohen-Zinder | Identification of a missense mutation in the bovine ABCG2 gene with a major effect on the QTL on chromosome 6 affecting milk yield and composition in Holstein cattle[END_REF]. Moreover, both variants were much less significant (-log 10 (P) = 7.4 for FC and 16.8 for PC) than the variant with the most significant effect on both traits, located in the HERC6 gene (intron) (-log 10 (P) = 9.7 for FC and 24.3 for PC). Thus, in our study ABCG2 is not the best candidate gene. However, we cannot completely exclude it because of its low MAF (0.02) and therefore its limited imputation accuracy, which may tend to underestimate its effect. For the QTL on BTA11 that affected FC in HOL cows, the 10 most significant variants were all located in the PAEP gene. Six of them were identified in a 1.5-kb stretch of the upstream region of the gene (103,299,655-103,301,229 bp), and were ranked from 1st (103,300,548 bp; rs109982707) to 8th in the peak; the 4th-ranked variant was in the 5′-UTR region (103,301,694 bp; rs41255686); the 6th-ranked variant was located in the downstream region (103,308,330 bp; rs109087963); the 9th-ranked variant was in a splicing region (103,304,656 bp; rs109990218); and finally, the 10th-ranked variant in the peak was a missense variant (103,303,475 bp; rs110066229). Together with another missense variant located at 103,304,757 bp (rs109625649), variant rs110066229 was previously identified as the functional mutation for protein variants A and B, which are associated with different levels of β-lactoglobulin in milk [START_REF] Ganai | Novel polymorphisms in the bovine beta-lactoglobulin gene and their effects on betalactoglobulin protein concentration in milk[END_REF]. Several nucleotide-level GWAS have found effects of this region on FC [START_REF] Pausch | Meta-analysis of sequence-based association studies across three cattle breeds reveals 25 QTL for fat and protein percentages in milk at nucleotide resolution[END_REF][START_REF] Jiang | Functional annotation and Bayesian fine-mapping reveals candidate genes for important agronomic traits in Holstein bulls[END_REF][START_REF] Raven | Targeted imputation of sequence variants and gene expression profiling identifies twelve candidate genes associated with lactation volume, composition and calving interval in dairy cattle[END_REF] or milk whey proteins [START_REF] Sanchez | Within-breed and multi-breed GWAS on imputed whole-genome sequence variants reveal candidate mutations affecting milk protein composition in dairy cattle[END_REF][START_REF] Sanchez | Sequence-based GWAS, network and pathway analyses reveal genes co-associated with milk cheese-making properties and milk composition in Montbéliarde cows[END_REF], and all have pointed to candidate variants in the PAEP gene. However, each of these studies highlighted a different best candidate variant, and these variants were always distinct from the two missense variants that cause the A and B protein polymorphisms. Moreover, Sanchez et al. [START_REF] Sanchez | Within-breed and multi-breed GWAS on imputed whole-genome sequence variants reveal candidate mutations affecting milk protein composition in dairy cattle[END_REF] found that a peak remained when one of the missense variants was fixed in the GWAS, which suggested that the missense variants described by Ganai et al. [START_REF] Ganai | Novel polymorphisms in the bovine beta-lactoglobulin gene and their effects on betalactoglobulin protein concentration in milk[END_REF] do not explain all the effects of this region on milk composition.

We also identified several genes involved in the metabolism of milk fatty acids (FASN, DGAT1, GPAT4, and MGST1) as good functional candidates to explain the changes observed in milk composition, and in each of these genes, we highlighted the most plausible candidate variants. FASN (fatty acid synthase) encodes a key enzyme in de novo fatty acid synthesis, whereas GPAT4 (glycerol-3-phosphate acyltransferase 4) is paralogous to DGAT1 (diacylglycerol O-acyltransferase 1), with the two genes occupying adjacent nodes of the mammary triglyceride synthesis chain [START_REF] Coleman | Enzymes of triacylglycerol synthesis and their regulation[END_REF]. The MGST1 (microsomal glutathione S-transferase 1) gene plays a role in oxidative stress reaction and although it has typically been associated with milk composition, and in particular with milk fat, its role in lipid metabolism is less clear. It has been shown to reduce lipid peroxidation products in human mammary cell culture [START_REF] Johansson | Multiple roles of microsomal glutathione transferase 1 in cellular protection: a mechanistic study[END_REF], but its functional impact on bovine milk production or composition traits has not been yet demonstrated.

The QTL region that was detected and validated at the centromeric end of BTA14 presented effects on different milk production traits, with the strongest effect on FC in the three breeds. With a frequency of 0.22, the A allele of the K232A mutation in DGAT1, which decreases FC, PC, and FY, and increases MY and PY [START_REF] Grisart | Positional candidate cloning of a QTL in dairy cattle: identification of a missense mutation in the bovine DGAT1 gene with major effect on milk yield and composition[END_REF], was the most significant variant for FC in HOL cows. It ranked 3rd and 13th in the peak for this trait and was much less polymorphic in NOR (MAF = 0.08) and MON (MAF = 0.007), respectively. In the vicinity of DGAT1, many genes have been annotated in the 0.5-Mb region between 1.5 and 2 Mb on BTA14. Our analyses indicated that the best candidate variants for many other traits in different breeds were located in other genes of this region (MROH1, bta-mir-1839, HSF1, RECQL4, MFSD3, GPT, CPSF1, ADCK5, and SLC39A4); further investigations could reveal, as has been suggested in many dairy cattle breeds and in particular in HOL, MON, and NOR [START_REF] Gautier | Characterization of the DGAT1 K232A and variable number of tandem repeat polymorphisms in French dairy cattle[END_REF], the existence of other causal mutations in this region.

We also identified two other candidate genes acting on FC in the QTL detected on BTA19 in HOL bulls and cows. In HOL cows, the variants with the most significant effects were both missense and located in the CCDC57 gene (rs41921161 at 51,319,797 bp, ranked 1st, and rs41921160 at 51,319,759 bp, ranked 2nd). However, five variants in the FASN gene ranked 5th to 9th in the peak with three located in the upstream region and two intronic. Among these, the upstream rs136067046 variant (at 51,383,847 bp, ranked 6th) was also the best candidate variant identified in a previous study for a QTL acting on milk fatty acid composition [START_REF] Sanchez | Sequence-based GWAS, network and pathway analyses reveal genes co-associated with milk cheese-making properties and milk composition in Montbéliarde cows[END_REF]. This region has been extensively studied for its effects on milk fat content and milk fatty acid composition. Although the role of FASN in the regulation of milk fat is more obvious than that of CCDC57, both genes are generally cited to explain the effects of this region [START_REF] Bouwman | Fine mapping of a quantitative trait locus for bovine milk fat composition on Bos taurus autosome 19[END_REF][START_REF] Goddard | Genetics of complex traits: prediction of phenotype, identification of causal polymorphisms and genetic architecture[END_REF][START_REF] Knutsen | Unravelling genetic variation underlying de novo-synthesis of bovine milk fatty acids[END_REF][START_REF] Palombo | Genome-wide association study of milk fatty acid composition in Italian Simmental and Italian Holstein cows using single nucleotide polymorphism arrays[END_REF].

In the three breeds studied here, we found a QTL on BTA27 that was also strongly associated with FC. The results of the cow GWAS directly pointed to five candidate variants, all located in the GPAT4 gene, which ranked in the top 5 in all three breeds. These variants, which were in complete LD in each of the three breeds, had a MAF ranging from 0.47 to 0.49 depending on the breed and are located in the upstream (rs211250281, rs378026790, rs209479876, and rs209855549) or the 5′-UTR (rs208675276) regions of GPAT4. GPAT4, also named AGPAT6, was previously described as a functional gene for milk fat content as well as protein and lactose contents by Littlejohn et al. [START_REF] Littlejohn | Expression variants of the lipogenic AGPAT6 gene affect diverse milk composition phenotypes in Bos taurus[END_REF]. These authors identified 10 linked variants associated with milk composition, which included the rs211250281, rs209855549, and rs208675276 variants found in the top 5 for each of the three breeds analyzed in our study. Moreover, the four variants located in the upstream region have been identified as candidate causal variants for FC in Holstein and Fleckvieh cows [START_REF] Daetwyler | Whole-genome sequencing of 234 bulls facilitates mapping of monogenic and complex traits in cattle[END_REF] whereas the top five variants were found to be the best candidates to explain variations in milk protein composition in a multi-breed analysis [START_REF] Sanchez | Within-breed and multi-breed GWAS on imputed whole-genome sequence variants reveal candidate mutations affecting milk protein composition in dairy cattle[END_REF] and variations in fat content in a meta-analysis [START_REF] Pausch | Meta-analysis of sequence-based association studies across three cattle breeds reveals 25 QTL for fat and protein percentages in milk at nucleotide resolution[END_REF]. All these results are consistent with the existence of a causative mutation located in the promoter region of GPAT4 which could regulate the expression level of this gene. Daetwyler et al. [START_REF] Daetwyler | Whole-genome sequencing of 234 bulls facilitates mapping of monogenic and complex traits in cattle[END_REF] suggested that the InDel rs378026790 was the most likely causal variant because of its high probability to overlap a transcription factor binding site but we cannot exclude rs208675276 which is in the 5′-UTR region and therefore closer to the transcription initiation site.

MGST1 has also been frequently described as a functional candidate gene for the QTL detected at ~ 94 Mb on BTA5 with effects on milk composition traits [START_REF] Sanchez | Within-breed and multi-breed GWAS on imputed whole-genome sequence variants reveal candidate mutations affecting milk protein composition in dairy cattle[END_REF][START_REF] Pausch | Meta-analysis of sequence-based association studies across three cattle breeds reveals 25 QTL for fat and protein percentages in milk at nucleotide resolution[END_REF][START_REF] Raven | Targeted imputation of sequence variants and gene expression profiling identifies twelve candidate genes associated with lactation volume, composition and calving interval in dairy cattle[END_REF][START_REF] Sanchez | Sequence-based GWAS, network and pathway analyses reveal genes co-associated with milk cheese-making properties and milk composition in Montbéliarde cows[END_REF][START_REF] Van Den Berg | Comparing power and precision of within-breed and multibreed genome-wide association studies of production traits using whole-genome sequence data for 5 French and Danish dairy cattle breeds[END_REF][START_REF] Iso-Touru | Genome-wide association analysis of milk yield traits in Nordic Red Cattle using imputed whole genome sequence variants[END_REF][START_REF] Littlejohn | Sequence-based association analysis reveals an MGST1 eQTL with pleiotropic effects on bovine milk composition[END_REF]. In the present study, the variant, which was shared between MON, NOR, and HOL cows and was most strongly associated with fat content or yield, was located in an intronic region of this gene (rs211210569 at 93,945,738 bp). This variant was also found to be responsible for effects on fat yield in the study of van den Berg et al. [START_REF] Van Den Berg | Comparing power and precision of within-breed and multibreed genome-wide association studies of production traits using whole-genome sequence data for 5 French and Danish dairy cattle breeds[END_REF] in both Danish and French Holstein bulls.

In addition to these good functional genes, we also identified and validated other promising genes for which the relationship with milk production or composition traits is less thoroughly understood. PICALM (phosphatidylinositol binding clathrin assembly protein), which was linked with PC in HOL (on BTA29), was previously associated with milk protein composition and lactose content [START_REF] Sanchez | Within-breed and multi-breed GWAS on imputed whole-genome sequence variants reveal candidate mutations affecting milk protein composition in dairy cattle[END_REF][START_REF] Sanchez | Sequence-based GWAS, network and pathway analyses reveal genes co-associated with milk cheese-making properties and milk composition in Montbéliarde cows[END_REF][START_REF] Lopdell | DNA and RNA-sequence based GWAS highlights membrane-transport genes as key modulators of milk lactose content[END_REF]. TBC1D22A (TBC1 domain family member 22A) was associated with PC in HOL and has been previously implicated in milk protein content [START_REF] Pausch | Meta-analysis of sequence-based association studies across three cattle breeds reveals 25 QTL for fat and protein percentages in milk at nucleotide resolution[END_REF][START_REF] Jiang | Functional annotation and Bayesian fine-mapping reveals candidate genes for important agronomic traits in Holstein bulls[END_REF]. VPS13B (vacuolar protein sorting 13 homolog B) had effects on FC, PC, and MY in our study and has been previously associated with milk fat and protein contents [START_REF] Fang | A Multiple-SNP approach for genome-wide association study of milk production traits in Chinese Holstein cattle[END_REF]. Finally, GRAMD4 (GRAM domain containing 4) had effects on PC in MON and was previously identified as a candidate gene for milk protein and mineral composition in the same breed [START_REF] Sanchez | Sequence-based GWAS, network and pathway analyses reveal genes co-associated with milk cheese-making properties and milk composition in Montbéliarde cows[END_REF].

The candidate variants that we identified in this study for both production and udder traits, which were mostly located in the non-coding regions of the genome, are either causative themselves or in LD with causative variants. The discovery of causal variants for complex traits remains challenging but should be facilitated in the next few years by two factors: (i) the most recent run of the 1000 Bull Genomes Project (run8 released in 2020), which contains, in total and within each breed, a larger number of bovine animals with whole-genome sequences that are aligned on the most recent ARS-UCD1.2 bovine genome assembly [START_REF] Rosen | De novo assembly of the cattle reference genome with single-molecule sequencing[END_REF] to enable more accurate imputation, and (ii) improved annotations of regulatory regions of the bovine genome, provided by the FAANG consortium [START_REF] Giuffra | Functional Annotation of Animal Genomes (FAANG): current achievements and roadmap[END_REF].

Conclusions

In the current study, GWAS analyses conducted on 10,871 bulls and 85,303 cows of the three main French dairy cattle breeds, Holstein, Montbéliarde, and Normande, enabled the identification and validation of 54 QTL for economically important traits related to milk production, udder morphology, and udder health. The first set of GWAS was carried out using whole-genome sequence data from bulls for the purpose of primary detection, and these enabled us to directly target candidate genes and candidate variants that were then added to the customized chip used for routine genomic evaluation of French dairy cattle. Analyses conducted in younger populations of cows then enabled us to validate a large number of these genes and variants, and yielded a more comprehensive understanding of the genetic determinism underlying these traits. Because they are now included on the genotyping chip, these candidate causative variants can be used for genomic predictions of production and udder traits in these three dairy cattle breeds.

Fig. 2

 2 Fig. 2 Sequence-based GWAS: -log 10 (P) plotted against the position on Bos taurus autosomes of variants linked with protein content in a Montbéliarde, b Normande, and c Holstein bulls

Fig. 3 Fig. 4

 34 Fig. 3 Sequence-based GWAS: -log 10 (P) plotted against the position on Bos taurus autosomes of variants linked with somatic cell score (SCS) in a Montbéliarde, b Normande, and c Holstein bulls

  a ID number associated with a group of QTL linked with milk production traits that had overlapping confidence intervals or less than 1 Mbp distance between the bounds of the confidence intervals b Individual QTL; milk yield (MY), fat content (FC), protein content (PC), fat yield (FY), protein yield (PY) c Confidence interval (CI) of the most significant individual QTL d ENSBTAG000000 e Effect of the variant, expressed in genetic standard deviation units f Extended confidence interval (CI) encompassing all the individual QTL detected in the same genomic region Table 6 QTL detected for udder traits in Montbéliarde (MON), Normande (NOR), and Holstein (HOL) bulls #: number of a ID number associated with a group of QTL linked with udder traits that had overlapping confidence intervals or less than 1 Mbp distance between the bounds of the confidence intervals b Individual QTL; somatic cell score (SCS), clinical mastitis score (CM), udder support (US), udder depth (UD), fore udder attachment (FUA), rear udder height (RUH), fore teat distance (FTD), udder balance (UB), teat orientation (TO), teat length (TL), and milking speed score (MSS) c confidence interval (CI) of the most significant individual QTL d ENSBTAG000000 e Effect of the variant, expressed in genetic standard deviation units f Extended confidence interval (CI) encompassing all the individual QTL detected in the same genomic region QTL ID a BTA Breed Traits (# QTL b ) CI of the TOP-CI c Variant with the most significant effect CI of the EXT-CI f

  with a group of QTL linked with milk production traits that had overlapping confidence intervals or less than 1 Mbp distance between the bounds of the confidence intervals b Individual QTL; milk yield (MY), fat content (FC), protein content (PC), fat yield (FY), protein yield (PY) c Candidate variants selected from sequence-based GWAS results

  Montbéliarde (MON), Normande (NOR), and Holstein (HOL) cows

Fig. 5

 5 Fig. 5 Validation GWAS: -log 10 (P) plotted against the position on Bos taurus chromosome 6 of variants linked with udder depth in a Montbéliarde and b Holstein cows, and c with somatic cell score (SCS) in Holstein cows; in black, 50K SNPs; in green, candidate variants; red circle indicates the variant with the most significant effect

Type of trait Trait and abbreviation Number of bulls with DYD Reliability of DYD mean (sd)

			MON	NOR	HOL	MON	NOR	HOL
	Milk production	Milk yield (kg) MY	2434	2175	6262	0.91 (0.09)	0.89 (0.11)	0.92 (0.05)
		Fat content (%) FC	2434	2175	6262	0.93 (0.08)	0.92 (0.10)	0.94 (0.04)
		Protein content (%) PC	2434	2175	6262	0.93 (0.08)	0.92 (0.10)	0.94 (0.04)
		Fat yield (kg) FY	2434	2175	6262	0.91 (0.09)	0.89 (0.11)	0.92 (0.05)
		Protein yield (kg) PY	2434	2175	6262	0.91 (0.09)	0.89 (0.11)	0.92 (0.05)
	Udder health	Clinical mastitis CM	1857	1427	4959	0.43 (0.21)	0.40 (0.22)	0.43 (0.21)
		Somatic cell score SCS	2438	2203	6318	0.87 (0.07)	0.85 (0.07)	0.88 (0.06)
	Udder morphology	Udder support US	2494	2180	6311	0.83 (0.07)	0.87 (0.06)	0.82 (0.08)
		Udder depth UD	2511	2020	6319	0.90 (0.05)	0.83 (0.07)	0.88 (0.06)
		Fore udder attachment FUA	2500	2164	5959	0.86 (0.07)	0.82 (0.07)	0.83 (0.08)
		Rear udder height RUH	2498	2147	6107	0.85 (0.07)	0.74 (0.10)	0.80 (0.09)
		Teat length TL	2515	-	6321	0.92 (0.05)	-	0.89 (0.05)
		Fore teat distance FTD	2509	2032	6319	0.89 (0.06)	0.86 (0.07)	0.88 (0.06)
		Udder balance UB	2478	2164	6275	0.77 (0.09)	0.81 (0.08)	0.81 (0.09)
		Teat orientation TO	2500	2175	6318	0.86 (0.06)	0.85 (0.06)	0.85 (0.07)
	Milking ease	Milking speed score MSS	2500	2164	6300	0.86 (0.07)	0.80 (0.08)	0.79 (0.09)

Table 2

 2

	Trait*	Montbéliarde			Normande				Holstein
		#QTL	TOT	Min	Max	#QTL	TOT	Min	Max	#QTL	TOT	Min	Max
	MY	1	1.6	1.6	1.6	1	2.0	2.0	2.0	3	10.0	1.1	7.8
	FC	5	19.7	1.8	11.1	3	23.2	1.7	18.4	8	37.0	0.5	26.8
	PC	6	17.2	1.2	7.0	5	20.0	2.0	8.0	11	27.7	0.7	8.2
	FY	1	1.1	1.1	1.1	1	3.0	3.0	3.0	4	13.8	0.6	10.9
	PY	0				0				4	3.7	0.3	1.7
	CM	0				0				0	0	0	0
	SCS	3	6.3	1.9	2.2	0				2	2.6	0.8	1.8
	US	1	4.3	4.3	4.3	1	2.5	2.5	2.5	1	1.6	1.6	1.6
	UD	2	4.9	1.8	3.1	0				4	6.3	0.8	1.7
	FUA	1	3.9	3.9	3.9	0				2	3.1	1.5	1.6
	RUH	1	3.2	3.2	3.2	1	5.2	5.2	5.2	2	3.9	1.0	2.9
	TL	0				-	-	-	-	1	3.4	3.4	3.4
	FTD	0				0				2	3.6	1.6	1.9
	UB	1	5.1	5.1	5.1	0				1	1.2	1.2	1.2
	TO	0				0				1	0.9	0.9	0.9
	MSS	2	3.9	1.9	2.0	0				2	2.2	1.0	1.2

Table 5 QTL detected for milk production traits in Montbéliarde (MON), Normande (NOR), and Holstein (HOL) bulls QTL ID a BTA Breed Traits (# QTL b ) CI of the TOP-CI c Variant with the most significant effect CI of the EXT-CI f Bounds (Mbp) # Genes Position (bp) Variant ID Functional annotation MAF -log 10 (P) b e

 5

	Functional annotation	TOP-CI: CI		EXT-CI: Extended
		of the QTL		CI
		with the most
		significant effect
		Number	%	Number	%
	Intergenic	6642	56.8	12,421	59.7
	Upstream	764	6.5	965	4.6
	Downstream	773	6.6	1039	5.0
	3' UTR	41	0.35	50	0.24
	5′ UTR	14	0.12	14	0.07
	Intronic	3286	28.1	6077	29.2
	Synonymous	105	0.90	133	0.64
	Non-coding transcript exon	4	0.03	6	0.03
	Splicing region	17	0.15	27	0.13
	Missense	50	0.43	66	0.32
	Total	11,696	100	20,798	100

Table 5 (continued) QTL ID a BTA Breed Traits (# QTL b ) CI of the TOP-CI c Variant with the most significant effect CI of the EXT-CI f Bounds (Mbp) # Genes Position (bp) Variant ID

 5

	Bounds (Mbp) # Genes
	sd
	MAF -log 10 (P) b e
	Functional annotation

#: number of

Table 7 GWAS validation results for production traits in Montbéliarde (MON), Normande (NOR), and Holstein (HOL) cows QTL ID a BTA Breed Traits (# QTL b ) Candidate variants c

 7

	1	2	2	0	5	3	4	8	1	1	1	3	0	0	2	1	0	0	0	3	11	3	2	4	24	4	8
	0.379 0.036 12.4-12.5	0.186 0.026 88.7-88.9	0.1787 0.026 88.7-88.9	0.272 0.034 71.9-72.2	-0.256 0.042 88.6-89.6 0.206 0.033 88.4-89.2	0.1789 0.03 88.4-89.2	-0.316 0.047 81.9-85.4 -0.718 0.114 23.4-23.7 -0.78 0.127 25.7-26.3 -0.859 0.147 28.1-29.2 -0.506 0.07 24.4-25.1 -0.283 0.038 62.6-62.7 -0.358 0.049 62.6-62.9 -0.272 0.042 62.5-62.9 -0.153 0.023 7.5-7.7 0.186 0.031 59.3-60.4	-0.199 0.032 60.5-61.5 0.1856 0.031 59.3-60.4	-0.192 0.028 25.6--28 0.326 0.053 33.2-35.3	-0.377 0.057 33.3-34.6 -0.302 0.051 38.5-38.9 -0.342 0.057 19.3-21.1 -0.167 0.023 49.1-51.5 0.2493 0.042 49.3-50.3	-0.347 0.051 48.3-51.2
	0.14 25	0.45 12.4	0.457 11.3	0.28 14.9	0.37 9.2	0.31 9.2	0.407 8.6	0.11 10.8	0.02 9.5	0.02 9.1	0.01 8.2	0.07 12.5	0.42 13.1	0.215 12.3	0.21 10.0	0.47 10.5	0.18 8.9	0.37 9.0	0.183 8.9	0.29 10.9	0.4 9.1	0.04 10.5	0.05 8.6	0.09 8.7	0.31 12.7	0.066 8.4	0.052 10.8
	12,441,072 rs109808447 TMTC2 Intron	88,800,994 rs110461240 ABCC9 intron	88,800,994 rs110461240 ABCC9 intron	72,022,013 rs136970036 Intergenic	88,744,985 rs110181141 Intergenic	88,881,928 rs474137839 Intergenic	88,697,293 rs110041776 GC intron	82,059,790 rs134625094 Intergenic	23,620,280 rs437752750 ME1 Intron	25,888,225 rs797068196 HEY2 upstream	28,197,098 rs455285205 38849 d intron	25,052,440 rs210030313 PLAG1 5′UTR /CHCHD7 upstream	62,679,331 rs109168890 Intergenic	62,714,882 rs109332098 Intergenic	62,694,032 rs109134926 Intergenic	7,565,865 rs445450010 Intergenic	60,151,690 rs109603247 Intergenic	60,523,834 rs134785404 Intergenic	60,151,690 rs109603247 Intergenic	27,082,437 rs209792391 Intergenic	34,288,737 rs382921722 Intergenic	34,310,163 rs482718265 Intergenic	38,725,076 rs133440951 Intergenic	20,158,547 rs43101108 Intergenic	49,421,004 rs208859984 KCNQ1 Intron	49,782,986 rs437410319 KCNQ1 Intron	49,299,622 rs210749543 NAP1L4 Intron
	1	2	2	0	1	3	4	3	1	1	1	3	0	0	2	1	0	0	0	2	11	3	2	4	4	4	8
	12.4-12.5	88.7-88.9	88.7-88.9	71.9-72.2	88.6-89.6	88.4-89.2	88.4-89.2	81.9-83.1	23.4-23.7	25.7-26.3	28.1-29.2	24.4-25.1	62.6-62.7	62.6-62.9	62.5-62.9	7.5-7.7	59.3-60.4	60.5-61.5	59.3-60.4	26.1-28.0	33.2-35.3	33.3-34.6	38.5-38.9	19.3-21.1	49.1-50.3	49.3-50.3	48.3-50.3
		UD(1)	FUA(1)	RUH(1)	UD(1)	SCS(1)	UD(1)	FTD(2)	SCS(1)	SCS(1)	SCS(1)	RUH(1)	FUA(1)	US(1)	US(1)	MSS(1)	MSS(1)	UD(1)	MSS(1)	FUA(3)	UB(1)	RUH(1)	TO(1)	MSS(1)	UD(3)	SCS(1)	UB(2)
		HOL	HOL	HOL	MON	HOL	HOL	HOL	MON	MON	MON	MON	MON	MON	NOR	HOL	HOL	MON	MON	HOL	MON	HOL	HOL	MON	HOL	HOL	HOL
		5	5	6	6	6	6	8	9	9	9	14	17	17	17	19	19	19	19	20	24	24	26	28	29	29	29
		42	42	43	44	45	45	46	47	48	49	50	51	51	52	53	54	55	55	56	57	58	59	60	61	61	61.

Acknowledgements

The authors gratefully acknowledge all their colleagues who participated in the BIG project; the French National Database (CTIG) and Valogene for

Funding

The bull genotypes were obtained by the CartoFine project funded by the French National Research Agency (ANR-05-GENANIMAL-007) and by APIS-GENE. About 15% of the cow genotypes originated from the LactoScan (ANR-08-GANI-034) and Amasgen (ANR-08-GENM-030-01) projects funded by the French National Research Agency and by APIS-GENE. This study was funded by INRAE under the BIG project.

Availability of data and materials

Phenotypes originated from the national database for genetic evaluation. Most of the cow genotypes originated from genomic selection programs that are managed by Valogene. All data belong to French farmers and cannot be disclosed without explicit authorization.

Supplementary information

Supplementary information accompanies this paper at https ://doi. org/10.1186/s1271 1-020-00575 -1.

Additional file 1: Table S1. Heritability of traits in Montbéliarde (MON), Normande (NOR), and Holstein (HOL) cattle.

providing phenotype and genotype data, respectively; and the contribution of the 1000 Bull Genomes Project.

Authors' contributions

MPS wrote the manuscript. TT and MPS led the analysis, synthesized and discussed the results, and designed the content of the article. MB managed sequence analyses of the 1000 Bull Genomes Project. PC, AB, DB, CH, and MPS developed computing programs. DB, CH, and MPS performed imputation analyses. TT, RL, CH, and MPS performed GWAS. SF and DB designed and managed the BIG project. All authors read and approved the final manuscript.

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

• fast, convenient online submission • thorough peer review by experienced researchers in your field

• rapid publication on acceptance

• support for research data, including large and complex data types • gold Open Access which fosters wider collaboration and increased citations maximum visibility for your research: over 100M website views per year

•

At BMC, research is always in progress.

Learn more biomedcentral.com/submissions

Ready to submit your research ? Choose BMC and benefit from:

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.