Goutam Sahana

Didier Boichard

Bernt Guldbrandtsen

Mogens Sandø Lund

Aurélien Capitan
email: aurelien.capitan@inra.fr

Md Mesbah-Uddin

Terhi Iso-Touru

Johanna Vilkki

Dirk - Jan

Supervisors Main supervisors

Fertility is an economically important trait in dairy cattle. Fertility is defined as the genetic ability of a cow to show oestrus and conceive after insemination, and to resume breeding after calving. This trait had a negative trend in recent years, partly due to intensive selection for production-related traits-which has a negative genetic correlation with fertility-using few bulls with high genetic merits for production traits. Before genomics, it was difficult to achieve faster genetic progress for fertility using pedigree-based breeding scheme due to low heritability of this trait. Genomic selection scheme provides a solution to this hard-to-select problem, especially for traits with low heritability. Identification of causal variants for recessive lethal mutations, when possible, and selection of a set of predictive markers that successfully track such causal variants, is vital for making mating-decision (e.g. to avoid atrisk mating) and for additional increase in genetic gains using genomic selection.

The overall aim of this PhD thesis is to identify causal variants for recessive lethal mutations and select a set of predictive markers that are in high linkage-disequilibrium with the causal variants for female fertility in dairy cattle. We addressed this broad aim under five articles/manuscripts that are presented through Chapter 2 to Chapter 6 in this thesis.

Chapter 2 describes a systematic approach of mapping recessive lethals in French Normande cattle using homozygous haplotype deficiency (HHD). This study shows the influence of sample size, quality of genotypes, quality of (genotype) phasing and imputation, age of haplotype (of interest), and last but not the least, multiple testing corrections, on discovery and replicability of HHD results. It also illustrates the importance of fine-mapping with pedigree and whole-genome sequence (WGS) data, (cross-species) integrative annotation to prioritize candidate mutation, and finally, large-scale genotyping of the candidate mutation, to validate or invalidate initial results.

Chapter 3 describes a high-resolution population-scale mapping of large chromosomal deletions from whole-genome sequences of 175 animals from three Nordic dairy breeds. This study employs three different approaches to validate identified deletions. Next, it describes population genetic properties and functional importance of these deletions. Finally, it illustrates deletion formation mechanisms based on the assembled sequence features at breakpoints. This study provides the basis for subsequent genetic studies performed during this PhD thesis.

Résumé in English iv

Chapter 4 deals with three main issues related to imputation of structural variants, in this case, large chromosomal deletions, e.g. availability of deletion genotypes, size of haplotype reference panel, and finally, imputation itself. To address the first two issues, this study describes a Gaussian mixture model-based approach where read-depth data from the variant call format (VCF) file is used to genotype a known deletion locus, without the need for raw sequence (BAM) file. Finally, it presents a pipeline for joint imputation of WGS variants along with large chromosomal deletions.

Chapter 5 describes genome-wide association studies for female fertility in three Nordic dairy cattle breeds using imputed WGS variants including large chromosomal deletions. This study is based on the analyses of eight fertility related traits using single-marker association, conditional and joint analyses. This study illustrates that inflation in association test-statistics could be seen even after correcting for population stratification using (genomic) principal components, and relatedness among the samples using genomic relationship matrices; however, this was known for traits with strong polygenic effects, among other factors. Finally, mapping of several new quantitative trait loci (QTL), along with the previously known ones, are reported in this study. This study also highlights the importance of including (imputed) large deletions for association mapping of fertility traits.

Chapter 6 describes prediction of genomic breeding values for fertility using SNP array-chip genotypes, selected QTL and large chromosomal deletion. Using genomic best linear unbiased prediction (GBLUP) method with one or several genomic-relationship matrices derived from a set of selected markers, this study reports higher prediction accuracy compared with previous report. This study also highlights the influence of selecting markers with best predictability, especially for a breed with small training population, in accuracy of genomic prediction. The results demonstrate that large deletions in general have a high predictive performance.

Finally, Chapter 7 provides a general discussion, conclusion, and perspective based on findings of this PhD study. v

Résumé på Dansk

Frugtbarhed er en økonomisk vigtig egenskab hos malkekvaeg. Frugtbarhed defineres som en koens evne til at viser brunsttegn, bliver draegtige efter insemination, og som kan genoptage avl efter kaelvning. Denne egenskab har de seneste år vaeret i tilbagegang på grund af kraftig selektion for produktionsrelaterede egenskaber, der har en negativ sammenhaeng med frugtbarheden. Det er svaert at selektere for (bedre) frugtbarhed ved hjaelp af stamtavlebaserede avlsprogrammer, da frugtbarhed har en lav arvelighed. Genomisk information kan muligvis løse dette problem. Identifikation af kausale markører, hvis muligt, og udvaelgelse af et saet praediktive markører, der med succes sporer de sande kausale alleler, er afgørende for at foretage parringsbeslutninger (f.eks. for at undgå parring mellem to baerere) og for at øge genetisk fremgang med genomisk

Acknowledgement

Firstly, I am grateful to the EGS-ABG consortium for giving me the opportunity in this double doctorate program, hosted by two leading universities in the field-Aarhus University and AgroParisTech. Last four years was a life-changing and eye-opening journey for me with lots of academic, research, cultural and international experiences that I could not imagine without this Erasmus Mundus endeavor. I am proud and feel lucky to be an EGS-ABG graduate. I am thankful to all EGS-ABG colleagues, past and present, for their help throughout this journey.

I am grateful to both of my PhD host institutes, QGG/Aarhus University and INRA/AgroParisTech, for offering this PhD. I am thankful to MBG/AU and GABI/INRA secretaries for their help regarding the administration. I have to mention one name, Karin Smedegaard (may she rest in peace). Her help and assistance started from receiving me from the Viborg train station, setting-up the MBG email, registering for social security, opening bank account … to visiting apartments and renting my first apartment in Viborg. I am grateful to her, will deeply miss her presence.

I was lucky to have both Goutam and Didier as my PhD supervisors. I am thankful to them for their help and guidance throughout my PhD. It was an excellent learning experience for me. I really enjoyed and appreciated their prompt feedback on my work (both research results and writing) with detailed and to-the-point suggestions, and constructive criticisms when needed. I am indebted to Goutam and Didier for their untiring efforts to educate me, to correct my mistakes, and finally to shape this dissertation in a presentable manner (just within one week!).

Specially, I want to thank Goutam under whose direct supervision I started my PhD journey and spent major part of the learning phase of my PhD. From day one, Goutam made it clear that I have to take charge of my PhD studies. No doubt, I was overwhelmed with this freedom and independence, but with his sincere supervision, I was able to navigate through my PhD, and gained enough confidence to pursue research topics from conception to dissemination. Before starting PhD, I had little experience working in Linux environment or with big data. I am indebted to Goutam for bearing with me, especially in the 1 st year of my PhD where he allowed me to attend several courses and spend time to learn programming/bioinformatics techniques (even though, I had no apparent progress in my research project). Without that, I could not imagine pursuing a PhD in quantitative genetics (let alone completing it!). I am thankful to Didier for his supervision during the second half of my PhD. His mentorship was pleasant, delightful and enlightening. Despite his busy schedule, I always had access and undivided attention from him. Every meeting and discussion I had with him, whether research or life in general, ended with enlightenment and wisdom, which greatly helped my research and personal life. Lastly, my stay in GABI was stress-free, and I was always in vacation mood-could be the Paris effect!!! I am also thankful to my co-supervisors Bernt, Mogens and Aurelien. I am grateful to Bernt from whom I learned the ins and outs of the Linux cluster in QGG, whole-genome sequence analysis, various aspects of genetics (to name few). Particularly, when I was struggling to interpret results in light of population/evolutionary genetics, I cherish the memories of those one-to-one discussions with Bernt (that often lasted several hours in late afternoon), which greatly helped me in finishing my first PhD paper. I am thankful to Mogens for bringing a broader, applied, and real-life problem solving perspective of the research to the supervision team, which was the North Star in my PhD that helped me not to fell prey of hair-splitting type analysis, and guided me to pursue topics that are relevant to livestock breeding and genetics. I am also grateful to Mogens for his keen interest in ensuring an allii round training for all QGG PhD fellows to prepare us for a better career in academia/research and/or industry. Finally, I am thankful to Aurelien from whom I learned a lot, especially, while working with the recessive lethal project. I admired his enthusiasm, friendly attitude and hands-on supervision approach; I really enjoyed working with him. Besides his guidance on mapping recessive lethals, I learned a lot on storytelling, and precise scientific writing-though it was lot of extra work to remove few of those "may be/could be" words (he would often say we have the data to check that...).

I am thankful to Grum, my QGG Buddy, for his help in social and work life. I also thank my office mates Lingzhao and Bingjie in QGG, and Rabia and Margarita in GABI/INRA for their pleasant company and stimulating work environment. I am grateful to all QGG/Aarhus University and G2B/INRA colleagues for their support throughout this journey.

Lastly, I am grateful to my family, parents, my wife Abeda, and son Abdullah, for their unconditional love, support and sacrifice. Without their amazing supports, this PhD would still be a dream. I dedicate this PhD Thesis to them.

Md Mesbah Uddin

Le chapitre 2 décrit une approche systématique de cartographie des variants létaux récessifs chez les bovins Normands français basée sur la recherche de déficit en haplotypes homozygotes (HHD). Cette étude montre l'influence de la taille de l'échantillon, de la qualité des génotypes, de la qualité du phasage des génotypes en haplotypes et de l'imputation, de l'âge de l'haplotype et enfin, de la définition des seuils de signification prenant en compte les tests multiples, sur la découverte et la reproductibilité des résultats de HHD. Elle illustre également l'importance de la cartographie fine avec les données de généalogie et de séquence de génome entier (WGS), l'annotation intégrative (entre espèces) pour hiérarchiser les mutations candidates et, enfin, le génotypage à grande échelle de la mutation candidate, pour valider ou invalider les mutations initiales. Cette étude met en évidence une mutation létale dans le gène CAD dont la fréquence dans la population est de 3%.

Le chapitre 3 décrit une cartographie à haute résolution de grandes délétions chromosomiques de séquences du génome dans une population de 175 animaux appartenant à trois races laitières nordiques. Cette étude utilise trois approches différentes pour valider les résultats de la viii cartographie. Le chapitre décrit les propriétés génétiques des populations et l'importance fonctionnelle des délétions identifiées. Enfin, il illustre les mécanismes de formation des délétions basés sur les caractéristiques des séquences d'ADN aux points de cassure. Cette étude fournit les bases pour les études génétiques ultérieures réalisées au cours de cette thèse sur le délétions.

Le chapitre 4 traite de trois questions liées à l'imputation de variants structuraux, ici de délétions chromosomiques importantes: la disponibilité des génotypes de délétion, la taille du panel de référence d'haplotypes et, enfin, l'imputation elle-même. Pour aborder les deux premières questions, cette étude décrit une approche basée sur un modèle de mélange gaussien dans laquelle les données de profondeur de lecture provenant de fichiers au format VCF (variant call format) sont utilisées pour génotyper un locus de délétion connu, en l'absence d'information sur la séquence brute. Enfin, il présente un pipeline pour l'imputation conjointe de variants WGS et de grandes délétions chromosomiques.

Le chapitre 5 décrit des études d'association pangénomiques de la fertilité femelle dans trois races de bovins laitiers nordiques à l'aide de variants WGS imputés et de grandes délétions chromosomiques. Cette étude concerne huit caractères de fertilité et utilise des analyses d'association mono-marqueur, conditionnelles et conjointes. Cette étude montre qu'une surestimation, ou « inflation », des statistiques de test peut être observée même après correction pour la stratification de la population à l'aide de composantes principales génomiques et pour les structures familiales à l'aide de matrices de relations génomiques. Ce biais était connu pour les caractères très polygéniques. Enfin, cette étude présente plusieurs locus de traits quantitatifs (QTL) nouveaux et confirme plusieurs autres déjà connus. Elle souligne également l'importance d'inclure les grandes délétions (imputées) pour la cartographie par association des caractères de fertilité.

Le chapitre 6 décrit la prédiction des valeurs génomiques de fertilité (ou indice de fertilité) à l'aide de génotypes à puces SNP, de QTL sélectionnés et de délétions chromosomiques importantes. En utilisant la méthode de meilleure prédiction linéaire sans biais génomique (GBLUP) avec une ou plusieurs matrices de relations génomiques dérivées d'un ensemble de marqueurs sélectionnés, cette étude rapporte une précision de prédiction améliorée. Cette étude met également en évidence l'influence de la sélection des marqueurs les plus prédictifs, en particulier pour une race ayant une population d'apprentissage réduite, sur la précision des ix prédictions génomiques. Enfin, les résultats démontrent que les grandes délétions ont en général un pouvoir prédictif élevé.

Enfin, le chapitre 7 propose une discussion générale, une conclusion et des perspective basées sur cette thèse.

x

Table of Contents

Table of Contents

Female fertility is an economically important trait in dairy cattle. Due to low heritability (0.01 to 0.04) of this trait, it is relatively hard to improve using traditional pedigree-based selection scheme. Fertility is negatively correlated with production traits such as milk and protein yield [START_REF] Pryce | Genetic aspects of common health disorders and measures of fertility in Holstein Friesian dairy cattle[END_REF][START_REF] Rauw | Undesirable side effects of selection for high production efficiency in farm animals: a review[END_REF][START_REF] Roxström | Genetic and Environmental Correlations Among Female Fertility Traits and Milk Production in Different Parities of Swedish Red and White Dairy Cattle[END_REF]. Therefore, an intensive selection for production traits, over the years, using few bulls with high breeding values resulted in a decrease in cows' fertility. The indirect response to selection for production traits explains only half of the reduction seen in fertility. Factors such as higher inbreeding rate [START_REF] Howard | Invited review: Inbreeding in the genomics era: Inbreeding, inbreeding depression, and management of genomic variability[END_REF]) also contributes to this decreasing trend, primarily due to embryonic lethality (VanRaden et al., 2011, Cole et al., 2016). However, selection based on genetic markers, with or without pedigree information, in a scheme known as 'genomic selection' [START_REF] Meuwissen | Prediction of total genetic value using genome-wide dense marker maps[END_REF], provides an opportunity to overcome this negative trend in female fertility (Garcia-Ruiz et al., 2016). For genomic selection, large reference population were assembled for major dairy breeds with high quality of phenotypes and genotypes. Due to their large size and their family structure (artificial insemination bulls with many daughters with fertility records), these reference populations allow for high fertility evaluation reliability, in spite of this low heritability. Noteworthy, these reference populations are also a unique resource for identifying the causal factors for recessive genetic defects and mapping quantitative trait loci (QTL).

Recessive lethals are responsible for a substantial economic loss in dairy cattle (Cole et al., 2016). To identify causal factors/QTL for recessive lethals, two approaches could be used.

First, searching for homozygous haplotype deficiency (HHD) in the population. This approach does not require phenotype. Here, when homozygotes for a common haplotype are absent or less frequent in the population than expected, may indicate embryonic lethality. This HHD approach has been shown to be very successful in cattle where genotyping is widely used and several discoveries were made in last few years. For example, seven recessive haplotypes in Holstein were discovered using HHD, namely, HH1 to HH5 (VanRaden et al., 2011, Fritz et al., 2013, Sahana et al., 2013, Daetwyler et al., 2014, McClure et al., 2014), HH6 (Fritz et al., 2018) and HH7 [START_REF] Hoze | Short Communication: A splice site mutation in CENPU is associated with recessive embryonic lethality in Holstein cattle[END_REF]. The second approach is performing genome-wide association studies (GWAS) for fertility traits. For example, in Nordic Red Dairy cattle, two large chromosomal deletions, one affecting fertility [START_REF] Schulman | Mapping of fertility traits in Finnish Ayrshire by genome-wide association analysis[END_REF], Kadri et al., 2014) and the other causing stillbirth (Sahana et al., 2016), were mapped using GWAS approach. In recent years, a third approach, next-generation sequencing (NGS)-based reverse-genetics also becoming a popular choice with lot of success both in humans [START_REF] Macarthur | A systematic survey of loss-of-function variants in human protein-coding genes[END_REF][START_REF] Rivas | Human genomics. Effect of predicted protein-truncating genetic variants on the human transcriptome[END_REF] and in animals (Charlier et al., 2016, Li et al., 2016, Michot et al., 2016). However, in this PhD thesis, we used former two approaches.

Although identification of causal variants is very challenging, the identification of recessive deleterious alleles is very important to manage their frequencies in a population. In genomic selection, use of causal variants information could theoretically provide higher accuracy of evaluation (van den Berg et al., 2016), and persistence of prediction accuracy over generations (Meuwissen and Goddard, 2010). In mating plans, this information can also be used to avoid at-risk mating and thus can reduce embryonic mortality in the population (Cole et al., 2016).

As whole-genome sequencing (WGS) provides a nearly complete list of variants, it may facilitate to identify causal variants. This is one of the major advantage of using WGS data over standardized genotyping panels-that include only a small proportion of the variants and are unlikely to include the causal ones. NGS is rapidly becoming the primary focus to characterize the genetic basis of quantitative traits. WGS could then be used to fine-map the causal mutations. This information could subsequently be implemented in the breeding plan to avoid carrier-to-carrier mating, thus improving fertility in dairy cattle.

The organization of this PhD thesis is as follows. In Chapter 1, I have introduced the topic of this thesis with relevant backgrounds followed by the overall aim of the PhD and five study specific objectives. In subsequent five chapters (Chapter 2 to Chapter 6), I have presented the main results and discussions of this PhD thesis under five separate articles. Finally, in Chapter 7, I have presented general discussions and future perspectives based on this PhD study.

Large-scale sequencing and genotyping

This PhD study utilized WGS and large-scale genotyping data that was generated/gathered through two primary international collaborations, namely, the 1000 Bull Genomes Project (1KBGP) (Daetwyler et al., 2014) and the EuroGenomics consortium (Lund et al., 2011, Boichard et al., 2018). The 1KBGP is a large-scale WGS consortium that provided a (nearly) full catalog of small sequence variants (e.g. SNPs and indels) for all modern cattle breeds, which is a valuable resource for identifying candidate mutation. The 1KBGP also provides a large reference panel for imputation of WGS variants into the SNP array-typed cattle populations. The EuroGenomics consortium was initiated with a goal of enlarging reference population for genomic selection. This consortium also designed a low-density custom SNP array, called EuroG10k, for large-scale genotyping. EuroG10k has one common part that contain SNPs from BovineLD, Bovine50k and BovineHD BeadChips to assure better imputation, and several private parts that contain various classes WGS markers, such as, QTL detected through GWAS, structural variants, and putative causal variants selected using NGSbased reverse-genetics approach (Boichard et al., 2018).

Genetic Markers

Today, most genetic studies, in humans (Genomes Project et al., 2015), mouse [START_REF] Li | Global genetic analysis in mice unveils central role for cilia in congenital heart disease[END_REF][START_REF] Bult | Mouse Genome Database (MGD)[END_REF], and major agricultural species, such as cattle (Daetwyler et al., 2014, Bouwman et al., 2018), sheep [START_REF] Jiang | The sheep genome illuminates biology of the rumen and lipid metabolism[END_REF], pig [START_REF] Li | Genomic analyses identify distinct patterns of selection in domesticated pigs and Tibetan wild boars[END_REF], chicken (Wong et al., 2004), wheat [START_REF] He | Exome sequencing highlights the role of wild-relative introgression in shaping the adaptive landscape of the wheat genome[END_REF][START_REF] Pont | Tracing the ancestry of modern bread wheats[END_REF], rice [START_REF] Zhao | Pan-genome analysis highlights the extent of genomic variation in cultivated and wild rice[END_REF], and many other species, are based on genetic markers derived from whole-genome sequencing (WGS). This revolution in genomics is due to a rapid advancement of WGS techniques coupled with gradual decline in both sequencing and genotyping cost [START_REF] Goodwin | Coming of age: ten years of next-generation sequencing technologies[END_REF]. In most of these studies, these markers include small sequence variants such as single-nucleotide polymorphism (SNP) and small (usually <50 bp) insertion and/or deletion (indel). Large (>50 bp) DNA polymorphisms known as structural variants (SVs) are studied less frequently, primarily because genotype calling is not straightforward and requires complex procedures. These SVs include deletions, duplications, insertions, inversions, and translocations. Unbalanced SVs, i.e. fewer or more copies of some sequences in an organism compared with the reference genome, are called copy-number variants (CNVs) (Alkan et al., 2011). Figure 1.1 illustrates some of these DNA polymorphisms and Figure 1.2 shows a generic pipeline for mapping and genotyping WGS variants.

In short, different methods can be used to call SV genotypes. Some methods rely on WGS data and use different kinds of information: variation in read depth, for insertions and deletions, unexpected read pairs (too close or too far to each other), split reads (when two parts of the same read map to different regions of the genome assembly). Some other methods rely on SNP arrays signals: different intensities may reflect copy number variants [START_REF] Kadri | Prediction of a deletion copy number variant by a dense SNP panel[END_REF]; Mendelian incompatibilities for close markers are often the consequences of deletions [START_REF] Capitan | A 3.7 Mb deletion encompassing ZEB2 causes a novel polled and multisystemic syndrome in the progeny of a somatic mosaic bull[END_REF]; finally, properly designed tests may be used to test breakpoints of structural variants (Boichard et al., 2018). This latter approach is very powerful as it can generate a very large number of reliable genotypes at a given variant and is a method of choice for targeted SV. The focus of several studies of this PhD study was on large chromosomal deletions due to its potential as causal variants especially for fitness related traits. Furthermore, there were fewer studies where large deletions or any other types of SVs were investigated beyond mapping their genomic locations and estimating their population frequencies. However, mapping SVs in itself is a challenging task (Bickhart and Liu, 2014). Moreover, SVs studies also suffers from high false discovery rates, and often requires the use of multiple SV callers to get a consensus Here, for breakpoints, SV call with small confidence interval is given preference over other calls, when there is a ~70% reciprocal overlap between the calls. 2011), numerous lethals were discovered in several dairy cattle breeds (Fritz et al., 2013, Sahana et al., 2013, Pausch et al., 2015, Kipp et al., 2016, Schwarzenbacher et al., 2016, Michot et al., 2017, Fritz et al., 2018).

Search for HHD is a vital tool for mapping genetic defects in situations where direct phenotypes are not available, such as for embryonic lethals. This search is quite straightforward. In most routine genomic evaluation procedures, genotypes are imputed and therefore phased. These phased data are used to estimate haplotype frequencies, expected and observed number of homozygous haplotypes, and to infer the status of each individual carrier. In the targeted region, candidate variants and genes can be identified by comparing the carrier-status at the candidate variant for sequenced individuals with their carrier-status at the haplotype. Among the usually limited number of candidate variants, the best one is selected based on its functional annotation, such as predicted functional effect, sequence conservation, and cross-species annotation. In many situations, the prioritized variant locate within the coding sequence of a gene and its annotation is accurate. Finally, this candidate variant (and possibly several others) can be added on the SNP chip used for routine genotyping in genomic selection program, providing large numbers animals with direct genotypes. This information is especially powerful to validate (or invalidate) the candidate variant.

HHD method is especially efficient when the haplotype frequency is moderate to high and linkage disequilibrium (LD) with the causal variant is strong. As an example, a few thousands genotypes are required to map an embryonic lethal at a frequency of 0.1, when 0 homozygote is observed instead several tens (1%) are expected. However, the expected number of homozygotes becomes very limited when the frequency is low (e.g. 1/10,000 when the frequency is 1%) and very large numbers, e.g. hundreds of thousands, of genotyped animals are often required to identify lethal haplotypes with low to moderate frequencies. Due to the availability of large-scale genotyping data in genomic selection, search for HHD is a routine practice for dairy cattle breeds in many countries. Furthermore, as the size of genotyped cohort increases over time, the power of detecting rare and low frequency causal variants with HHD also increases. For example, [START_REF] Hoze | Short Communication: A splice site mutation in CENPU is associated with recessive embryonic lethality in Holstein cattle[END_REF] recently mapped a candidate mutation in CENPU gene with an allele frequency of 0.8% in French Holstein that was tagged by Holstein haplotype 7 (HH7) with a population frequency of 1.1%. In many situations, these lethal variants were found to be 'untolerated' mutations in the coding sequence, in genomic regions highly conserved in the evolution, and results in a truncated, incomplete, and usually nonfunctional protein product. For instance, several of these mutations responsible for early embryonic mortality affect cell division (Fritz et al., 2018) or DNA base synthesis pathways (Michot et al., 2017, Mesbah-Uddin et al., 2019), i.e. basic mechanisms of life.

Genotype Imputation

Genotype imputation is a critical step in most modern genetic studies. Indeed genomic information deeply varies across individuals from low density SNP chip to WGS at various depths, i.e. from several thousand to several tens of millions markers. Inferring the complete information for all individuals is a prerequisite for subsequent genetic analysis such as genomewide association studies or genomic evaluation. Imputation is a statistical approach to predict genotypes for un-typed markers based on a densely genotyped haplotype reference panel relative to the imputation target population. Because these methods use haplotype information, mostly they also include the genotype phasing process when genotyped are unphased.

There are several tools available for genotype phasing and imputation with very different approaches (Marchini andHowie, 2010, van Leeuwen et al., 2015). The methods derived from human genetics, such as Beagle (Browning and Browning, 2016), Minimac (Das et al., 2016), and IMPUTE2 [START_REF] Marchini | A new multipoint method for genome-wide association studies by imputation of genotypes[END_REF] are the oldest. These tools rely on LD at the population level and use Bayesian algorithms. More recently, algorithm more suited to livestock population were developed, taking advantage of the small effective population size and the large half-sib family size, i.e. the long-range LD and the small number of haplotypes segregating in the populations. These methods build the library of the haplotypes observed in the reference population and then use this library to fill the missing information. In addition, they efficiently use pedigree information. FImpute (Sargolzaei et al., 2014) is the most widely used software for phasing and imputation in cattle populations. The relative efficiency of imputation strongly varies according to the size of the population, to its structure, extent of LD, and marker density. A combination of methods for phasing and then imputation is sometimes found optimal (Delaneau et al., 2014). FImpute predicts (usually) the most likely genotypes; however, FImpute does not provide imputation accuracy information as an output, contrary to Beagle, Minimac and IMPUTE2. Bayesian methods provide expected genotypes (i.e. allele dosages) which can take into account the uncertainty of the imputed genotypes and therefore, more appropriate than 'best-guess' genotypes for downstream analysis like association studies.

Genome-wide association study (GWAS)

Genome-wide association study (GWAS) is a method of choice for identifying quantitative trait loci (QTL) in virtually all species, due to the advent of high throughput genotyping and sequencing techniques. In GWAS, a dense array of genetic markers, either directly genotyped or imputed, is used to capture a substantial proportion of common genetic variation of the given species [START_REF] Mccarthy | Genomewide association studies for complex traits: consensus, uncertainty and challenges[END_REF]. GWAS exploits LD pattern at population level, i.e. the nonrandom association of alleles at different loci [START_REF] Slatkin | Linkage disequilibrium--understanding the evolutionary past and mapping the medical future[END_REF], to identify association of markers or haplotypes with the trait of interest (and QTL, if any). In ideal situation, these chromosome segments usually carry identical markers or haplotypes, and hence are expected to carry QTL alleles [START_REF] Hayes | Overview of Statistical Methods for Genome-Wide Association Studies (GWAS)[END_REF].

Over the years, GWAS has proven to be an efficient approach for gene mapping both in monogenic and polygenic traits. For example, Charlier et al. (2008) first reported fine-mapping of causal variants for five monogenic traits in cattle. Similarly, earlier studies by [START_REF] Barendse | A validated whole-genome association study of efficient food conversion in cattle[END_REF], [START_REF] Kolbehdari | A whole-genome scan to map quantitative trait loci for conformation and functional traits in Canadian Holstein bulls[END_REF], [START_REF] Daetwyler | A genome scan to detect quantitative trait loci for economically important traits in Holstein cattle using two methods and a dense single nucleotide polymorphism map[END_REF] and [START_REF] Druet | Fine mapping of quantitative trait loci affecting female fertility in dairy cattle on BTA03 using a dense single-nucleotide polymorphism map[END_REF] also showed success in mapping variants associated with polygenetic traits. Following those early studies, several hundreds of GWASs were performed in cattle over the last decade.

There are numerous statistical methods for association mapping. Today, single-marker association using mixed linear model (MLM) is a method of choice for most GWAS. In MLM analysis, false discoveries could be limited by accounting for both population stratification (usually by using few principal components (PCs) of the genomic relationship matrix) and relationships among individuals (by using pedigree or genomic kinship matrix) [START_REF] Zhang | Mixed linear model approach adapted for genome-wide association studies[END_REF]). Henderson's matrix notation of the MLM is as follows:

𝑦𝑦 = 𝑋𝑋𝑋𝑋 + 𝑍𝑍𝑍𝑍 + 𝑒𝑒
where y is the vector of observed phenotypes; u is a vector of random additive genetic effects;

X and Z are the known design matrices; and e is a vector of residual effects. β is a vector containing the fixed effects: in addition to the effect of the genetic marker of interest, it includes other fixed effects such as the population mean, sex, age, PCs, etc. The u and e vectors are assumed to be normally distributed with a null mean and a variance of:

𝑣𝑣𝑣𝑣𝑣𝑣 � 𝑍𝑍 𝑒𝑒 � = � 𝐺𝐺 0 0 𝑅𝑅 �,
where 𝐺𝐺 = 2𝐾𝐾𝜎𝜎 𝑎𝑎 2 with σa 2 as the additive genetic variance and K as the kinship matrix. If homogeneous variance is assumed for the residual effect, R = Iσe 2 , where σe 2 is the residual variance (it includes all other variances that is not captured by σa 2 ), and I is an identify matrix.

The proportion of the total variance explained by the additive genetic variance is defined as

narrow sense heritability, ℎ 2 = 𝜎𝜎 𝑎𝑎 2 𝜎𝜎 𝑎𝑎 2 +𝜎𝜎 𝑒𝑒 2 .
False positive association. Two major sources of false positives in GWAS is population stratification and hidden relatedness in the mapping population. These biases cause spurious association in GWAS. However, in MLM, correction for population structure and relatedness could be done using first few PCs and genomic relationship matrix.

Multiple-testing correction. Typically, hundreds of thousands to several millions of hypothesis testing are performed to identify markers associated with a trait. Therefore, stringent p-value threshold is required to account for these multiple testing and reduce false discoveries (type 1 error). Bonferroni correction is most commonly used multiple testing adjustment: here, p-value is divided by the total number of tests performed. In GWAS, a commonly used p-value threshold is <5×10 -8 , at 5% level of significance and ~1 million independent tests (i.e. 0.05/10 6 ).

Conditional GWAS. Single-marker GWAS often provides hundreds of loci with significant p-value (given the study had enough power to capture true QTL, i.e. able to reject null hypothesis when alternative is true). Conditional association analysis could then be performed for a given locus or for a chromosome to identify independent GWAS signal(s). Here, marker with the best association signal is used as a cofactor for subsequent round of GWAS.

Meta-analysis. Meta-analysis is a powerful tool to boost QTL detection power by combining GWAS summary statistics from multiple populations [START_REF] Willer | METAL: fast and efficient meta-analysis of genomewide association scans[END_REF]; examples of this approach in cattle include multi-breed meta-analysis for cattle stature (Bouwman et al., 2018) and milk protein composition in dairy cattle [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]. Following the method by [START_REF] Bolormaa | A multi-trait, meta-analysis for detecting pleiotropic polymorphisms for stature, fatness and reproduction in beef cattle[END_REF], multi-trait meta-analysis is also used in cattle to identify QTL affecting multiple traits, i.e. pleotropic QTL; example of this approach in cattle includes meta-analysis for mammary gland morphology [START_REF] Pausch | A multi-trait meta-analysis with imputed sequence variants reveals twelve QTL for mammary gland morphology in Fleckvieh cattle[END_REF].

Genomic Prediction

Early prediction of genetic merits of an organism is of vital importance for any breeding program concerning agricultural species. Availability of large-scale SNP array genotypes and WGS data provides tremendous advantage on selecting breeding population based on their breeding values estimated from genomic information, with or without pedigree information, rather than selecting solely based on performances and pedigrees. Genetic gain (yearly), ∆𝐺𝐺 = 𝑖𝑖𝑖𝑖𝜎𝜎 𝑎𝑎 𝐿𝐿 , depends on four parameters: selection intensity (i), selection accuracy (r), generation interval (L), and additive genetic standard deviation (𝜎𝜎 𝑎𝑎 ). Genomic selection can affect at least three of these parameters. In its simplest form, substantial reduction in generation time is possible since genomic evaluation can be performed at birth. In cattle, traditional selection was late due to progeny testing of bulls and performance testing of cows. [START_REF] Schaeffer | Strategy for applying genome-wide selection in dairy cattle[END_REF] showed that genomic selection could double genetic gain by substantial reduction of generation interval in dairy cattle. However, this situation is not general to all species, especially those with shorter generation interval, which explains why genomic selection is so popular in cattle. Genomic selection can also affect selection accuracy if the reference population is large, and selection intensity if genotyping is cheap, allowing large-scale evaluation of many candidates.

Nevertheless, genomic selection has brought a revolutionary change in selection and breeding for all major agricultural species (Georges et al., 2019).

The accuracy of prediction is a crucial component for any successful genomic selection program. Numerous studies were performed over the last few years. To improve prediction accuracy, different approaches were tested:

-some studies focused on the prediction methods [START_REF] Meuwissen | Prediction of total genetic value using genome-wide dense marker maps[END_REF][START_REF] De Los Campos | Prediction of complex human traits using the genomic best linear unbiased predictor[END_REF][START_REF] Kemper | Improved precision of QTL mapping using a nonlinear Bayesian method in a multi-breed population leads to greater accuracy of across-breed genomic predictions[END_REF][START_REF] Moser | Simultaneous discovery, estimation and prediction analysis of complex traits using a bayesian mixture model[END_REF]: different priors were used, from a constant Gaussian prior to all SNP (GBLUP) to a prior reflecting an oligogenic determinism;

-few focused on training population characteristics, such as size, genetic relationship with the candidates, and genetic parameters of the trait, such as heritability and genetic correlation (Calus et al., 2018, van den Berg et al., 2019); -others focused on selecting better predictors [START_REF] Sørensen | Genomic Feature Models[END_REF], Brøndum et al., 2015[START_REF] Vanraden | Selecting sequence variants to improve genomic predictions for dairy cattle[END_REF], Hay et al., 2018), i.e. not only generic SNP but variants with a higher probability of biological effects. These predictors were selected on the basis of their functional annotation or on the effect on the phenotype, as estimated in independent populations.

However, there are many opportunities yet to be explored to improve prediction accuracy further. Besides SNPs and indels, inclusion of other classes of genetic markers, such as CNVs, also needs further exploration. Feature selection using reverse-genetic approach could also be very useful for prediction purposes, as shown by Charlier et al. (2016), [START_REF] Michot | A reverse genetic approach identifies an ancestral frameshift mutation in RP1 causing recessive progressive retinal degeneration in European cattle breeds[END_REF] and Boichard et al. (2018). Last but not the least, automated feature selection using machine learning approaches can also be considered, such as using deep learning approaches as shown by [START_REF] Bellot | Can Deep Learning Improve Genomic Prediction of Complex Human Traits? Genetics[END_REF] for predicting complex traits in humans, although the success so far was limited.

Aim and objectives of this PhD study

The overall aim of this PhD thesis was to identify causal variants for deleterious recessive mutations and select a set of predictive markers for female fertility in dairy cattle. These results support recessive embryonic lethality with nearly complete penetrance, as was previously reported in CAD mutants in several Eukaryote species. The only homozygous cow had an extremely poor udder conformation, suggesting a potential role of CAD in udder development, but no effect was detected when comparing daughter yield deviations of 250 heterozygous bulls with that of 2,912 homozygotes for the ancestral allele. Together, this study showed the importance of large-scale screening for HHD with hundreds of thousands of animals, validating results with independent dataset, and considering unexpected live homozygotes, to avoid both false positive and false negative discoveries. These discoveries will be used primarily in mating decisions to avoid at-risk mating. In addition, we recommend including CAD in the breeding objectives of Normande cattle.

Key words: embryonic lethal, CAD, homozygous haplotype deficiency, dairy cattle, large-scale genotyping

Introduction

Inbreeding and genetic drift has increased substantially in dairy cattle over the last decades due to the extensive use of a few elite artificial insemination sires [START_REF] Weigel | Controlling Inbreeding in Modern Breeding Programs[END_REF][START_REF] Bjelland | Evaluation of inbreeding depression in Holstein cattle using whole-genome SNP markers and alternative measures of genomic inbreeding[END_REF]. This situation can increase the frequency of some deleterious alleles in the population, and thus unmask recessive genetic defects affecting various fitness traits (Cole et al., 2016).

Homozygosity mapping is an efficient approach to map such recessive deleterious alleles with few inbred live cases (Charlier et al., 2008). However, it requires the prior detection of live animals showing distinctive features, which is not possible for all types of genetic defects (e.g. when homozygotes die very early in the gestation period or when they display symptoms that can be confounded with non-genetic diseases). In 2011, VanRaden et al. proposed a new approach for identifying recessive genetic defects from large-scale screening of homozygous haplotype deficiency (HHD). Embryonic lethality or increased perinatal mortality is suspected when observed homozygotes are significantly less frequent than expected based on population Chapter 2 21 frequency. This approach enabled the detection of a number of haplotypes with both complete and incomplete penetrance in several breeds (Sahana et al., 2013, Pausch et al., 2015, Schwarzenbacher et al., 2016a), though it required tens or hundreds of thousands of genotyped individuals for rarer disorders.

However, some of these haplotypes could not be confirmed in independent studies, possibly due to false positive results or population structure within a breed. For example, haplotypes HH2 in Holstein and BH1

in Brown Swiss, showed significant HHD in the USA (VanRaden et al., 2011) but not in Europe (Fritz et al., 2013, Segelke et al., 2016, Fritz et al., 2018). For this reason, and since no phenotypic association and candidate mutation have been reported so far (McClure et al., 2014), they have been recently removed from the list of genetic defects considered by the National Association of Animal Breeders (NAAB) of the USA (Olivier Bulot, personal communication).

Previous studies also highlighted the importance of large-scale genotyping to validate candidate causal mutations (Schwarzenbacher et al., 2016a, Schwarzenbacher et al., 2016b, Michot et al., 2017, Fritz et al., 2018). Doing so, our team recently excluded a candidate variant in the SHBG gene for MH1 haplotype in Montbeliarde cattle (Fritz et al., 2013) and validated a new causal variant in the PFAS gene (Michot et al., 2017).

In 2013, Fritz et al. reported six haplotypes (NH1 to NH6) showing significant HHD in 11,466 Normande cattle, including four also showing some reduction in conception rate in at-risk versus control mating. In this large-scale screen, we were interested in validating these haplotypes in Normande cattle, detecting novel loci, and when possible, identifying causal variants.

Methods

Study Samples and Genotype Data

We analyzed 77,815 Normande cattle from the French genomic evaluation database, genotyped between the year 2008 and 2018 with different Illumina SNP-chips (Illumina Inc., San Diego, CA, USA), such as, 777K

BovineHD BeadChip, 50K BovineSNP50 BeadChip [START_REF] Matukumalli | Development and characterization of a high density SNP genotyping assay for cattle[END_REF], BovineLD BeadChip (Boichard et al., 2012a), and EuroG10K BeadChip (Boichard et al., 2018). In the total dataset, 27,632 animals had both parents genotyped, and 49,690 had sire and maternal grandsire (MGS) genotyped. The others, i.e. 353 and 140 animals with no or a single parent genotyped, respectively, were included as sires or MGS only. Since this dataset included all the animals studied by Fritz et al. (2013), we also performed a specific analysis on this subpopulation. This dataset comprised of 2,303 animals with both parents genotyped, and 9,163 animals with sire and MGS (but not dam) genotyped, as well as 111 sires or MGS with no or a single parent genotyped.

The quality control of the genotype data was carried out following the French genomic evaluation pipeline.

This includes: (i) filtering of low-quality genotypes, (ii) verifying the pedigree for Mendelian inconsistency, and (iii) phasing and imputing genotypes (Sargolzaei et al., 2014) from low to medium-density (for details see Boichard et al. (2012b)). After quality control, 43,801 autosomal SNPs from the BovineSNP50 BeadChip were kept for subsequent analysis. The genomic coordinates presented in this study refer to bovine genome assembly UMD3.1 [START_REF] Zimin | A whole-genome assembly of the domestic cow, Bos taurus[END_REF].

Screening for Homozygous Haplotype Deficiency

Following the method used by VanRaden et al. (2011) and Fritz et al. (2013), we screened the 29 bovine autosomes, with sliding windows of 20 consecutive markers, for homozygous haplotype deficiency (HHD).

We counted occurrences of every haplotypes, and calculated population frequency from the maternal chromosomal phases. In addition, we estimated the expected number of homozygotes using within-family transmission probability, such as, from carrier sires and dams to offspring, or from carrier sires and MGS (with un-genotyped dam) to offspring. For subsequent analyses, we considered haplotypes with population frequency greater than 1.0%, representing on average 17 different haplotypes for each 20-marker window (range 4 to 34 haplotypes). In total, we performed 717,747 tests to identify HHD, and therefore considered haplotypes with P < 1.39 × 10 -8 for further analysis, after correcting for multiple testing at 1% level of significance. The expected number of homozygotes (lambda) assuming neutrality was computed from sire, dam or MGS genotypes and frequency in the population as in Fritz et al. (2013). We calculated the probability of observing "q" homozygotes with expectation of "lambda", using Poisson distribution "ppois(q, lambda)" function in R software version 3.4.3 (R Core Team, 2017). After multiple testing correction, we considered haplotypes for further analysis when the ratio of observed to expected homozygotes was less than 0.25.

Comparison of Haplotype Status for NH1-NH6 between the 2013 and Present Studies

Haplotype status as computed by Fritz et al. (2013) were available for 2,364 AI bulls out of the 11,466 animals initially studied. We compared carrier status between the two studies for these animals.

Analysis of the Survival of Animals Homozygous for NH2

Information on the date of birth, date of death and cause (slaughter or natural death) was extracted from the French genomic evaluation database for 228 females homozygous for NH2 and 8,871 non-homozygous paternal half-sisters. In both groups, we evaluated the proportion of animals that died of natural causes at one, two and three years of age.

Evaluation of Haplotype Effect on Fertility Traits

We evaluated the haplotype effect on two 'binary' female fertility traits, conception rate (CR) and non-return rate at 56 days (NRR56), separately in heifers and lactating cows, using phenotypic records from the French genomic evaluation database (period January 2000 to May 2018). The CR is a measure of success ('1') or failure ('0') of each insemination, assessed by a calving after a compatible gestation length. NRR56 measures whether a female is re-inseminated within 56 days after a previous insemination. It is coded '0' if at least one insemination is recorded within this period and '1' otherwise. As both phenotypes are routinely evaluated in the French national evaluation system, non-genetic effects estimates (such as, effects of herd × year, year x month of insemination, parity, days in milk, AI technician, use of sexed semen, etc.) were available and were used to pre-adjust the phenotypes of these two traits. The mating class was coded '1' for at-risk mating (between carrier bulls and daughters of carrier bulls), and '0' for non-risk mating (all other combinations of genotypes regarding the bulls and the sires of the dams). We analyzed the phenotypic effect of the haplotype between at-risk mating vs non-risk mating using the following fixed effect model:

𝑌𝑌 𝑖𝑖𝑖𝑖𝑖𝑖 = µ + 𝑠𝑠𝑠𝑠𝑣𝑣𝑒𝑒 𝑖𝑖 + 𝑚𝑚𝑣𝑣𝑚𝑚𝑠𝑠𝑚𝑚𝑚𝑚_𝑐𝑐𝑐𝑐𝑣𝑣𝑠𝑠𝑠𝑠 𝑖𝑖 + 𝑒𝑒 𝑖𝑖𝑖𝑖𝑖𝑖
where, Yijk represents the phenotype of interest (i.e. adjusted CR or NRR56), µ is the overall phenotypic mean, sirej is the fixed sire effect, mating-classk is the fixed effect of the mating status, and eijk is the random residual error. The analysis was performed using the GLM procedure of SAS software v9.4 'Proc GLM'

(SAS Institute Inc., NC, USA). In addition, we calculated the expected decrease in fertility, assuming full lethality under homozygosity, using the formula:

1 4 ( 1 2-𝑓𝑓 ℎ𝑎𝑎𝑎𝑎_𝑘𝑘
)µ, where fhap_k is the frequency of the haplotypek, and µ is the phenotypic mean for CR or NRR56 (as in Michot et al. (2017)).

Fine-Mapping of the Normande Haplotype 7 (NH7) Locus

To fine map the causal mutation associated with the newly detected NH7 haplotype, we followed the approach used by [START_REF] Sonstegard | Identification of a nonsense mutation in CWC15 associated with decreased reproductive efficiency in Jersey cattle[END_REF] and Fritz et al. (2018). Briefly, this consisted in identifying phenotypically normal animals presenting a run-of-homozygosity (ROH) on either side of the candidate haplotype among the inbred descendant of the bull ("Nonnic", NORFRAM002977029292), the most ancient and influential carrier of NH7 in our dataset. Such asymptomatic individuals are very unlikely to be homozygous for the causal mutation, and thus regions with long ROH can be ruled-out for causal variant screening.

Analyzing Whole-Genome Sequence (WGS) Data for Candidate Mutations

We analyzed WGS of 2,333 animals (44 Normande cattle and 2,289 animals from other cattle breeds) from the 1000 Bull Genomes Project (Run6; 1KBGP) (Daetwyler et al., 2014) for identifying putative causal variant. Among the WGS animals, only two Normande bulls carried NH7 haplotype and the rest did not. We performed Pearson's correlation between haplotype status and allele dosage for all bi-allelic single nucleotide polymorphisms (SNPs) and small insertions and deletions (indels) within NH7 (±500 kb), to identify highly associated variants. Both, haplotype status and allele dosage were coded as '0', '1', or '2', when the animal carried either zero, one or two copies of alternative allele/candidate haplotype, respectively.

We presented squared values of Pearson's correlation (r 2 ) in a scatterplot.

We also analyzed structural variant (SV) calls on 389 whole-genome sequenced French cattle, which included 44 WGS Normande cattle and 345 animals from other breeds. The SVs were identified using three different software, namely BreakDancer v1.3.6 [START_REF] Chen | BreakDancer: an algorithm for highresolution mapping of genomic structural variation[END_REF], Delly v0.6.1 [START_REF] Rausch | DELLY: structural variant discovery by integrated paired-end and split-read analysis[END_REF], and Pindel v2.5 [START_REF] Ye | Pindel: a pattern growth approach to detect break points of large deletions and medium sized insertions from paired-end short reads[END_REF], and consensus SV calls from these tools were included for the analysis (for details see Letaief et al. (2017)). For candidate SVs, we focused our search within NH7 and ±500 kb surrounding regions. We checked for the concordance between the carrier status of NH7 and SVs.

Integrative Annotation of Candidate Mutations and Nearby Genes

We annotated all strongly correlated variants (r > 0.8) using Variant Effect Predictor (VEP v87) software (McLaren et al., 2016), which includes annotation of sequence ontology, SIFT score [START_REF] Kumar | Predicting the effects of coding non-synonymous variants on protein function using the SIFT algorithm[END_REF] for missense variants, overlapping gene name, etc. We retrieved three gene-level constraint scores: (i) the ratio of non-synonymous (dN) to synonymous (dS) substitution rate (dN/dS scores) for cow-human 1-to-1 orthologs, using BioMart (Kinsella et al., 2011[START_REF] Zerbino | Ensembl 2018[END_REF], (ii) Residual Variation Intolerance Score (RVIS) [START_REF] Petrovski | Genic intolerance to functional variation and the interpretation of personal genomes[END_REF] and (iii) missense Z-score [START_REF] Lek | Analysis of proteincoding genetic variation in 60,706 humans[END_REF] for human orthologs of cattle genes. We converted genomic coordinates of variants from "UMD3.1/bosTau6" to "GRCh37/hg19", for retrieving variant-level conservation scores using LiftOver software (http://genome.ucsc.edu/cgibin/hgLiftOver) [START_REF] Casper | The UCSC Genome Browser database: 2018 update[END_REF]. We retrieved evolutionary conservation statistics, such as, Genomic

Evolutionary Rate Profiling (GERP++) "rejected substitutions" scores for 35-mammalian alignments [START_REF] Davydov | Identifying a high fraction of the human genome to be under selective constraint using GERP++[END_REF], PhyloP and PhastCons scores for 100-vertebrate alignments [START_REF] Pollard | Detection of nonneutral substitution rates on mammalian phylogenies[END_REF], using UCSC Table Browser (http://genome.ucsc.edu/cgi-bin/hgTables) (Karolchik et al., 2004).

Multiple-Sequence Alignment

We performed multiple-sequence alignment among 24 different species, representing six groups, e.g.

vertebrates, metazoan, plants, fungi, protists, and bacteria. We extracted protein sequences from either Ensembl Genomes [START_REF] Kersey | Ensembl Genomes 2018: an integrated omics infrastructure for non-vertebrate species[END_REF] or UniProt (The UniProt, 2017). We used CLUSTAL Omega version-1.2.4 [START_REF] Sievers | Fast, scalable generation of high-quality protein multiple sequence alignments using Clustal Omega[END_REF] for multiple-sequence alignment. For the candidate gene, we also retrieved functional annotation of protein domains from NCBI's Conserved Domain Database (CDD) [START_REF] Marchler-Bauer | CDD/SPARCLE: functional classification of proteins via subfamily domain architectures[END_REF] [accessible from "https://www.ncbi.nlm.nih.gov/cdd"].

Large-Scale Genotyping of the CAD Candidate Mutation in French Cattle Breeds

A molecular test was designed for genotyping a candidate missense mutation in the CAD gene (g.72399397A>G on Chromosome 11, CAD p.Tyr452Cys) with the EuroG10K BeadChip (Boichard et al., 2018) routinely used for genomic evaluation. Genotypes of 381,916 animals, from 20 French cattle breeds, were available for this variant.

Genotyping of the CAD Mutation by PCR and Sanger Sequencing

We used PCR and Sanger sequencing to genotype the same mutation in a single individual found to be homozygous for the NH7 haplotype based on 50K genotypes, and a non-carrier animal. DNA was obtained from LABOGENA DNA (Jouy-en-Josas, France). We designed primers using Primer3 software [START_REF] Untergasser | Primer3--new capabilities and interfaces[END_REF]. The forward and reverse primers are 5′-AGGGGATCACAGAATTGGCA-3′ and 5′-GCCTCTGTCTCCTTGGGATT-3′, respectively. We performed PCR using "GoTaq Flexi DNA

Polymerase" (Promega, Madison, WI, USA) in Mastercycler Pro (Eppendorf, Montesson, France) thermocycler, following manufacturer's instructions and sent amplicons to Eurofins MWG (Hilden, Germany) for sequencing. We next used novoSNP software [START_REF] Weckx | novoSNP, a novel computational tool for sequence variation discovery[END_REF] for variant calling and visualization.

Phenotypic Examination of the Homozygous Cow and Paternal Half-Sisters

Clinical examination of the NH7 homozygote cow, named "Genetique", was not possible since this animal was slaughtered at the beginning of her first lactation due to poor udder conformation (personal communication from the breeder). However, own performance records for 34 morphological traits were available for this cow and its 2,239 paternal half-sisters from the French National databases (see Supplementary Figure S1). Individual performances were ranked for each trait and the position of "Genetique" was expressed in term of percentile.

Effect of the CAD Mutation on Udder Conformation Traits

Since "Genetique" displayed extremely poor udder conformation, we evaluated the effect of the CAD missense mutation, on four udder traits, namely, teat direction (TD), teat placement (TP), udder cleft (UC), and udder depth (UD). These traits are initially recorded in a scale of 1 (low) to 9 (high) (see ICAR (2018) for details), before being pre-adjusted for non-genetic factors and mean-centered in the French national evaluation system. Daughter yield deviations (DYD) were available for 250 heterozygous carriers and 2,912 non-carriers of the CAD mutation. The phenotypic effect of the mutation was estimated using the following fixed effect model:

𝑌𝑌 𝑖𝑖𝑖𝑖𝑖𝑖 = µ + 𝑠𝑠𝑠𝑠𝑣𝑣𝑒𝑒 𝑖𝑖 + 𝐶𝐶𝐶𝐶𝐶𝐶_𝑚𝑚𝑒𝑒𝑚𝑚𝑔𝑔𝑚𝑚𝑦𝑦𝑔𝑔𝑒𝑒 𝑖𝑖 + 𝑒𝑒 𝑖𝑖𝑖𝑖𝑖𝑖
Here, Yijk represents the phenotype of interest (i.e. DYD for TD, TP, UC, or UD), µ is the overall phenotypic mean, sirej is the fixed effect of the sire of the bull, CAD_genotypek is the fixed effect of the genotype (k = "AG" or "AA"), and eijk is the random residual error. The analysis was performed using the GLM procedure of SAS software v9.4 (Proc GLM)' (SAS Institute Inc., NC, USA).

Results and Discussion

Homozygous Haplotype Deficiency in Normande Cattle

In a sliding window approach, we analyzed 77,815 Normande cattle phased and imputed for 43,801 autosomal markers, to identify homozygous haplotype deficiency (HHD). Surprisingly, none of the six Normande haplotypes (NH1 to NH6) previously reported by Fritz et al. (2013) showed significant HHD in this follow-up study after correction for multiple testing (Table 2.1).

To elucidate the causes of these discrepancies, we reanalyzed the mapping population of 11,466 Normande cattle from Fritz et al. (2013) using the current French phasing and imputation dataset. Contrary to Fritz et al. (2013), the current dataset do not comprise the markers of the Illumina BovineSNP50v1 BeadChip, which have been removed from the subsequent versions of the chip for technical reasons. This concerns markers ARS-BFGL-NGS-97976 and BTA-58084-no-rs within NH1, Hapmap59512-rs29027076 within NH3, BTB-01218649 within NH5 and BTA-37438-no-rs and ARS-BFGL-NGS-29141 within NH6 haplotype. The removal of these problematic markers, which likely caused local phasing and imputation errors, explains the substantial changes in haplotype frequencies (Table 2.1) and carrier status (Table S2) observed between the two studies. Insufficient correction for multiple testing (arbitrary threshold of P < 1 × 10 -4 instead of Bonferroni correction) was also another reason for the false discovery. Among the six haplotypes, only NH2 was significant for HHD after Bonferroni correction (Table S1). 1.7×10 -10 1 Start and End positions represent first and last 50K-marker of the haplotype; positions correspond to bovine genome assembly UMD3.1 2 Information within the parenthesis and Normande haplotype NH1 to NH6 are from Fritz et al. (2013) † P-values are from Poisson distribution, calculated using R software's "ppois(q=Observed, lambda=Expected)" function * Genotypes for markers ARS-BFGL-NGS-97976 and BTA-58084-no-rs within NH1, Hapmap59512-rs29027076 within NH3, BTB-01218649 within NH5 and BTA-37438-no-rs and ARS-BFGL-NGS-29141 within NH6 haplotype were considered in Fritz et al. (2013) but not in the present study.

Finally, it is worth mentioning that the population analyzed by Fritz et al. (2013) consisted in animals aged one year or more whereas the present dataset includes numerous young animals genotyped in their first months of life, which is expected to reduce the power of detection when mapping loci responsible for juvenile mortality. To verify if NH2 is associated with postnatal mortality, and potentially explain why it is not detected in the present study, we analyzed the survival of 228 homozygous females and 8,871 nonhomozygous paternal half-sisters. We found no difference between the two groups at one, two and three years of age (Supplementary Table S3). In this context, we assume that NH2 was a false positive due to a lower power of detection in the analysis of Fritz et al. (2013) as compared with that of the current study.

Indeed the present dataset is sevenfold larger than the previous one, and the frequency of NH2 more than doubled in the population, which increased the probability of observing homozygotes.

Importantly, we detected two novel haplotypes with significant deficit of homozygotes (P < 1.39×10 -8 ). One, named NH7, is located on chromosome 11 between positions 70,750,209 and 72,476,622 bp, and the other, named NH8, on chromosome 24 between positions 52,606,336 and 54,542,635 bp. Haplotype frequencies in Normande breed were 3.7% and 1.9% for NH7 and NH8, respectively. In our dataset, we observed one homozygote for NH7 and seven homozygotes for NH8, while the expectations under neutrality were 57 and 40, respectively (Table 2.1). The oldest carriers in our dataset were born in 1977 for NH7 ("Newgate" NORFRAM002277007498 and "Nonnic" NORFRAM002977029292, two sons of "Valhalla" NORFRAM007668031767 who was not genotyped) and in 2003 for NH8 ("Uvray" NORFRAM002951401444). NH8 is a relatively new haplotype in Normande cattle, created by a recombination in the maternal gamete that was transmitted to "Uvray", with fewer opportunities to be observed in the homozygous state in the population. However, we considered both haplotypes for analyzing their effect on CR and NRR56. 2 Expected effects on CR and NRR56 were calculated using the formula:

1 4 � 1 2-𝑓𝑓 ℎ𝑎𝑎𝑎𝑎_𝑘𝑘
� µ ; here fhap_k is the frequency of NH7 or NH8 in the population, and µ is the phenotypic mean of CR or NRR56.

The NH7 but not the NH8 Haplotype Compromises Conception Rate in Normande Cattle

To elucidate whether NH7 and NH8 have negative effect on CR and NRR56, we analyzed ~10 million and ~6 million mating records, respectively, between genotyped bulls and daughters of genotyped bulls (Table 2.2).

For NH7, we observed a marked decrease in CR of 5.53% in heifers and 5.56% in cows, which is close to the expectation (7.30% and 6.01%, respectively). We also observed a moderate drop in NRR56 of 1.76% in heifers and 1.52% in cows, much lower than the 9.75 and 8.55% expected values in both groups assuming complete lethality in homozygotes. The re-insemination within 56 days (NRR56) is an indirect measure of embryonic losses within 35 days of gestation, since the average estrus cycle in cow is 21 days [START_REF] Reece | Dukes' physiology of domestic animals[END_REF]. Our analysis shows that on average only 18% of the expected embryonic losses occurred within 35 days of gestation, while the majority of the failure may occurred later in the gestation. Here, we suspect a gradual embryonic loss over the length of gestation, rather than abrupt mass loss in the first month after fertilization.

We also analyzed the effect of NH8 on CR and NRR56, and found that this haplotype does not affect these two traits in Normande cattle (Table 2.2). It is worth mentioning that detection power was limited by the reduced number of at risk mating (0.06% vs. >0.5% of total mating for NH7) in our dataset because of the low frequency (1.9%) and the recent spreading of this haplotype in the population. Another possibility is that NH8 cause postnatal mortality, possibly with incomplete penetrance. Seven homozygotes, all inbred descendants of "Uvray", have been observed for this haplotype. These were all alive and aged between 11 months and 3 years at the time of this study. Despite NH8's significant HHD, it is necessary to accumulate more carrier-to-carrier mating, and inheritance in more generations, to be confident about the effect and HHD. Considering this and the fact that juvenile mortality is out of the scope of our study, we focused the subsequent analyses on NH7.

The Homozygote for NH7 has Particularly Poor Udder Morphology

We found a single cow, named "Genetique", homozygous for NH7 haplotype and for a 20-Mb surrounding segment (~4.2 Mb upstream and 14.8 Mb downstream of the source haplotype) due to recent inbreeding. This individual deserves special attention. If it is homozygous at the causal mutation, penetrance is not complete but the animal may be affected by severe abnormality. Alternatively, it could be homozygous for the haplotype but not for the mutation, by carrying an ancestral version of the haplotype before the mutation event along with the lethal haplotype. The breeder informed us that this animal had normal health but was slaughtered during its first lactation due to very poor udder morphology. To verify this, we analyzed own performance records for 34 morphological traits on "Genetique" and its 2,239 paternal half-sisters. We confirmed that it was in the normal range, except for four udder traits, e.g. teat placement, teat direction, udder depth, and udder cleft, for which it ranked amongst the lowest 1, 5, 2, and 4% of its siblings, respectively (Figure 2.1, and Supplementary Figure S1). Furthermore, "Genetique" ranked amongst the lowest 1% on a linear combination of the scores of the four traits (Figure 2.1).

Fine-Mapping of the NH7 Locus in a ~1.01Mb Genomic Region

To fine-map the NH7 locus, we analyzed pedigree and found that among genotyped animals, half-brothers "Newgate" (NORFRAM002277007498) and "Nonnic" (NORFRAM002977029292) were the two most ancient carriers of the haplotype. Among them, however, "Nonnic" is the ancestor of many bulls with high breeding value. For example, "Elixir" (NORFRAM 005389014161), a descendent of "Nonnic", is the sixth most influential bull of the Normande breed with 4.7% contribution in the population [START_REF] Idele | Indicators of genetic variability -Normande[END_REF]. Therefore, we considered only phenotypically normal descendants of "Nonnic" (and "Elixir") to fine-map the NH7 locus.

We identified two cows, "0520" (NORFRAF003708210520) and "Fanion" (NORFRAF001447822444), presenting a run-of-homozygosity on either side of the candidate haplotype (Figure 2.2). Both inherited NH7 paternally from "Nonnic" (through "Elixir") and NH7 recombinant haplotypes maternally. In dairy cattle, females were rarely genotyped and therefore, it is difficult to track maternally inherited recombinant haplotypes [START_REF] Adams | Identification of a nonsense mutation in APAF1 that is likely causal for a decrease in reproductive efficiency in Holstein dairy cattle[END_REF]. Nevertheless, it is likely, from pedigree and sire-haplotype analyses, that "Nimbus" (Figure 2.2a) received the recombinant haplotype from "Elixir" and passed it to "0520". On the other hand, "Fanion" received a portion of the recombinant haplotype from "Nonnic", though the origin and source of the remaining portion is inconclusive due to un-genotyped female ancestors (Figure 2.2a).

It is noteworthy that both cows have normal health and udder morphology. In addition, during this study, "0520" and "Fanion" were in second and fourth lactation, respectively. Such phenotypically normal animals are very unlikely to be homozygous for genomic region harboring the causal mutation, unless the mutation is not fully penetrant. Under this assumption, the critical interval is defined by the two recombination points located: (i) between marker rs29022293 (Chr11:71,465,910) and rs29010798 (Chr11:71,520,448) in "Fanion", and (ii) between marker rs110176879 (Chr11:72,159,929) and rs109299742 (Chr11:72,476,622), in "0520". This enabled us to fine-map the locus from a ~1.73Mb initial region to a ~1.01Mb mutationcritical region (Figure 2.2b).

Identification of Causal Variant from WGS Data

To identify putative causal variants, we considered bi-allelic WGS variants (SNPs and indels), and large SVs.

We used SV calls on 389 WGS French cattle (2 NH7 carriers and 387 non-carriers) to identify SVs overlapping the NH7 haplotype and/or ±500 kb surrounding regions. We identified two deletions and one inversion within this interval, but none of them was present in the genomes of the NH7 carriers (Supplementary Table S4).

Then, we performed Pearson's correlation between NH7 carrier status of 2,333 animals from the 1KBGP (2 carriers and 2,331 non-carriers) and allele dosages for 33,640 bi-allelic WGS variants within the haplotype (±500 kb). Here, we identified eight strongly correlated (Pearson's correlation, r > 0.8) variants within the given interval (Figure 2.3). All these variants have high confidence call (avg. phred-scaled quality=162; min. 97), with sufficient mapping quality (avg. MQ=57; min. MQ=33) and read coverage (avg. coverage = 11.75x; min. 8.95x). Interestingly, five out of the eight variants were within our fine-mapped critical region, which could include the causal variant. However, it is known from studies that the best candidates do not always have perfect correlation with the causal haplotype (Michot et al., 2017, Fritz et al., 2018), which could be due to various artifacts (e.g. errors in phasing and imputation, errors in sequencing and WGS genotype calls, incomplete linkage disequilibrium with causal variant, etc.). We therefore considered all the eight variants for subsequent gene-and variant-level prioritization.

Identification of Causal Variant from WGS Data

To identify putative causal variants, we considered bi-allelic WGS variants (SNPs and indels), and large SVs. We used SV calls on 389 WGS French cattle (2 NH7 carriers and 387 non-carriers) to identify SVs overlapping the NH7 haplotype and/or ±500 kb surrounding regions. We identified two deletions and one inversion within this interval, but none of them was present in the genomes of the NH7 carriers (Supplementary Table S4).

Then, we performed Pearson's correlation between NH7 carrier status of 2,333 animals from the 1KBGP (2 carriers and 2,331 non-carriers) and allele dosages for 33,640 bi-allelic WGS variants within the haplotype (±500 kb). Here, we identified eight strongly correlated (Pearson's correlation, r > 0.8) variants within the given interval (Figure 2.3). All these variants have high confidence call (avg. phred-scaled quality=162; min. 97), with sufficient mapping quality (avg. MQ=57; min. MQ=33) and read coverage (avg. coverage = 11.75x; min. 8.95x). Interestingly, five out of the eight variants were within our fine-mapped critical region, which could include the causal variant. However, it is known from studies that the best candidates do not always have perfect correlation with the causal haplotype (Michot et al., 2017, Fritz et al., 2018), which could be due to various artifacts (e.g. errors in phasing and imputation, errors in sequencing and WGS genotype calls, incomplete linkage disequilibrium with causal variant, etc.). We therefore considered all the eight variants for subsequent gene-and variant-level prioritization. (Chr11:72,524,223) to rs110443631 (Chr11:78,468,627). The recombination occurred between marker rs29022293 (Chr11:71,465,910) and rs29010798 (Chr11:71,520,448) in "Fanion", and between marker rs110176879 (Chr11:72,159,929) and rs109299742 (Chr11:72,476,622), in "0520". The resulting critical region is located between marker rs29022293 and rs109299742, i.e. between chromosomal interval of 71,465,910 and 72,476,622 bp. Here, "NH7" denotes "Normande haplotype 7"; "?" denotes un-genotyped animal.

Figure 2.3. Identification of putative causal variants within NH7 locus (±500kb) using Pearson's correlation between haplotype status vs allele dosage of WGS variants. Each dot in the scatterplot indicates a WGS variant. Perfectly correlated (r 2 =1) CAD missense mutation is indicated by an arrow (with annotation).

We annotated these eight WGS variants using the Variant Effect Predictor (VEP) software tool (version 87) from Ensembl (McLaren et al., 2016). There were one intergenic variant, two intronic variants in PLB1, three intronic variants in BRE, and one intronic and one missense variants in CAD gene (Table 2.3). Next, to prioritize the causal variant from the list, we assessed functional impact of these variants using three genelevel-, e.g. dN/dS, RVIS, and Z-score, and four variant-level-, e.g. SIFT, GERP, PhyloP and PhastCons, conservation scores.

To identify selective constraint on genes, we first analyzed the dN/dS ratio of cow-human 1-to-1 orthologs. Though all the three genes are somewhat constrained (dN/dS < 1), the selective constraint is stronger on CAD and BRE (dN/dS < 0.1). We later focused on RVIS and Z-score for human homologs of those genes, and found that CAD is the most constraint gene among the three, in terms of tolerance to functional variants (Table 2.3).

We next focused on variant-level conservation scores, e.g. SIFT, GERP, PhyloP, and PhastCons, to elucidate the strength of evolutionary constraint on those positions. SIFT score represents conservation at proteinlevel and predict whether an amino acid substitution is deleterious (SIFT < 0.05) or tolerated (SIFT > 0.05) for protein function. SIFT predicted the missense mutation CAD g.72399397T>C (p.Tyr452Cys) as "deleterious" for the protein function (Table 2.3). In contrast, GERP, PhyloP, and PhastCons represent conservation at nucleotide-level, generated from multi-species whole-genome sequence alignments, where higher positive scores indicate stronger evolutionary constraint. We were able to retrieve conservation scores for seven positions, using "LiftOver" and "UCSC Table Browser", while the remaining one did not have any corresponding position from UMD3.1/bosTau6 to GRCh37/hg19. Interestingly, all three classifiers predicted CAD g.72399397T>C as the most constrained position, where damaging mutations are purged from the population by strong purifying selection (Table 2.3). 1 Chromosomal positions correspond to bovine genome assembly UMD3.1 2 dN/dS scores for human-cow 1-to-1 orthologs using BioMart (Kinsella et al., 2011) from Ensembl [START_REF] Zerbino | Ensembl 2018[END_REF] 3 Residual Variation Intolerance Score (RVIS) is retrieved from RVIS v4 (based on ExAC release 2.0; last accessed on 10 September 2018 from: http://genic-intolerance.org/) [START_REF] Petrovski | Genic intolerance to functional variation and the interpretation of personal genomes[END_REF]. Positive and negative RVIS scores represent tolerant and intolerant genes, respectively. RVIS "x%": it represents that the respective gene is amongst the "x" percentile of most intolerant human genes. 4 Missense Z-scores were retrieved from http://exac.broadinstitute.org/ [START_REF] Lek | Analysis of proteincoding genetic variation in 60,706 humans[END_REF] Integrating gene-and variant-level conservation scores is an effective approach to prioritize causal variants from human genomes [START_REF] Petrovski | Genic intolerance to functional variation and the interpretation of personal genomes[END_REF]. Although bovine genome annotations are somewhat limited to sequence ontology and SIFT scores, incorporating cross-species annotations from well-studied genes or sequence homologs in humans improved our prioritization of the causal variant. Finally, the perfect correlation with NH7, localization within the critical interval, along with strong gene-and variant-level evolutionary conservation, taken together, strongly support CAD g.72399397T>C (p.Tyr452Cys) as the causal variant tagging the (nearly) lethal haplotype in Normande cattle.

The Tyrosine Residue at CAD p.Tyr452Cys is Fully Conserved from Humans to Bacteria

Pyrimidines are important building blocks of DNA and RNA, as well as precursors for uridine-diphosphate (UDP) sugars that play a vital role in post-translational modification, such as, UDP-dependent glycosylation of proteins [START_REF] Jones | Pyrimidine nucleotide biosynthesis in animals: genes, enzymes, and regulation of UMP biosynthesis[END_REF][START_REF] Fairbanks | Importance of ribonucleotide availability to proliferating T-lymphocytes from healthy humans. Disproportionate expansion of pyrimidine pools and contrasting effects of de novo synthesis inhibitors[END_REF]. Pyrimidines can either be synthesized de novo from small precursors (e.g. glutamine, ATP, and HCO3 -), such as, in S-phage of proliferating cells [START_REF] Mitchell | De novo pyrimidine nucleotide biosynthesis in synchronized rat hepatoma (HTC) cells and mouse embryo fibroblast (3T3) cells[END_REF], or through salvaging of nucleosides (e.g. uridine), such as, in resting or matured cells [START_REF] Jones | Pyrimidine nucleotide biosynthesis in animals: genes, enzymes, and regulation of UMP biosynthesis[END_REF][START_REF] Fairbanks | Importance of ribonucleotide availability to proliferating T-lymphocytes from healthy humans. Disproportionate expansion of pyrimidine pools and contrasting effects of de novo synthesis inhibitors[END_REF]. In mammals, the first three reactions in the de novo biosynthesis of pyrimidine are catalyzed by CAD-a tri-functional cytosolic protein composed of carbamoyl-phosphate synthetase 2 (CPSase-2), aspartate transcarbamylase (ATCase), and dihydroorotase (DHOase) domains [START_REF] Jones | Pyrimidine nucleotide biosynthesis in animals: genes, enzymes, and regulation of UMP biosynthesis[END_REF][START_REF] Moreno-Morcillo | Structural Insight into the Core of CAD, the Multifunctional Protein Leading De Novo Pyrimidine Biosynthesis[END_REF]). An illustration of bovine CAD protein domains, from NCBI's Conserved Domain Database (CDD), is presented in Figure 2.4a.

Our candidate CAD mutation is located in CPSase-2 domain (Figure 2.4). We suspect that the wild-type allele (p.Tyr452 residue) is vital for CPSase-2 domain's activity, and therefore should be under strong evolutionary constraint. To assess the selective constraint on CAD, we performed multi-species sequence alignment among 24-species using Clustal Omega software [START_REF] Sievers | Fast, scalable generation of high-quality protein multiple sequence alignments using Clustal Omega[END_REF]. We found that this T-to-C transition mutation (CAD g.72399397T>C) produces a tyrosine-to-cysteine substitution (p.Tyr452Cys) on

CPSase-2 domain, replacing a fully conserved amino acid residue among the 24-species analyzed (Figure 2.4b) distributed in all life kingdoms.

Large-Scale Genotyping Confirmed Recessive Inheritance of CAD p.Tyr452Cys in Normande Cattle

We genotyped 33,323 Normande cattle for Chr11 g.72399397T>C (CAD p.Tyr452Cys) mutation (corresponding test in the Chip: g.72399397A>G in the TOP format) using Illumina EuroG10K BeadChip (Boichard et al., 2018) (Table 2.4). In this dataset, the NH7 haplotype and g.72399397G allele frequencies were 2.58 and 2.92%, respectively (Table 2.5). We did not observe homozygotes for g.72399397G allele, though 28 were expected under neutrality (P < 7×10 -13 ), which largely corroborate with the observed recessive-lethal inheritance of NH7 (Table 2.1). This CAD mutation was also strongly associated with NH7 haplotype status; Pearson's correlation between g.72399397G allele dosages and haplotype status is, r = However, the association was not perfect: 237 heterozygotes for the CAD mutation do not carry NH7 haplotype and eight heterozygous for NH7 do not carry the CAD mutation (Table 2.5). Such imperfect association can arise from recombination between NH7 and other haplotypes, along with other reasons, such as, genotyping, phasing or imputation errors. Further investigations revealed that approximately 67%

(158/237) of the heterozygotes for the CAD mutation but not for NH7 carry the same recombinant NH7 haplotype with a long stretch of source haplotype downstream. In addition, the eight animals heterozygous for CAD mutation do not carry the original NH7 but a recombinant identical-by-state (IBS) portion of the haplotype.

Similarly, previous studies showed that haplotypes are not a perfect proxy for causal mutations (Michot et al., 2017, Fritz et al., 2018), and there were instances where two different versions of the same haplotype segregate in the population, i.e. one with-and other without-causal mutation (Kipp et al., 2016[START_REF] Menzi | A transposable element insertion in APOB causes cholesterol deficiency in Holstein cattle[END_REF]). Yet, haplotype is a good substitute to avoid carrier-to-carrier mating, when causal mutation is unknown (VanRaden et al., 2011, Sahana et al., 2013), or when test results are confidential due to patents (Cole et al., 2016).

We next analyzed 348,593 additional animals (from 19 other French cattle breeds), genotyped for genomic selection, to elucidate whether this mutation is also observed in other populations. We found that the CAD mutation is only segregating in Normande cattle (Table 2.4).

Surviving CAD Homozygotes Exhibit Broad Range of Symptoms Depending on the Species and Nature of Mutation

For verification, we genotyped "Genetique", the cow homozygous for a 20-Mb region including NH7, through PCR and Sanger sequencing and found that it was also homozygous for our candidate mutation (Figure 2.5). Despite being alive and apparently healthy and fertile, this animal had very poor udder morphology (scored amongst the lowest 1% of animals with same sire as discussed earlier; Figure 2.1). This is consistent with previous observation in model organisms and humans where deleterious mutation in CAD lead to embryonic lethality with incomplete penetrance [START_REF] Norby | The biochemical genetics of rudimentary mutants of Drosophila melanogaster. I. Aspartate carbamoyltransferase levels in complementing and non-complementing strains[END_REF][START_REF] Falk | Sex-linked auxotrophic and putative auxotrophic mutants of Drosophila melanogaster[END_REF][START_REF] Willer | Analysis of the Zebrafish perplexed mutation reveals tissuespecific roles for de novo pyrimidine synthesis during development[END_REF][START_REF] Franks | C. elegans pharyngeal morphogenesis requires both de novo synthesis of pyrimidines and synthesis of heparan sulfate proteoglycans[END_REF][START_REF] Cox | Novel role for carbamoyl phosphate synthetase 2 in cranial sensory circuit formation[END_REF][START_REF] Ng | Biallelic mutations in CAD, impair de novo pyrimidine biosynthesis and decrease glycosylation precursors[END_REF]. Noteworthy, studies in human and Drosophila demonstrated that uridine supplementation reduces the severity of the symptoms by salvage pathway of pyrimidine biosynthesis [START_REF] Falk | Sex-linked auxotrophic and putative auxotrophic mutants of Drosophila melanogaster[END_REF][START_REF] Ng | Biallelic mutations in CAD, impair de novo pyrimidine biosynthesis and decrease glycosylation precursors[END_REF][START_REF] Koch | CAD mutations and uridine-responsive epileptic encephalopathy[END_REF]. Depending on the species and nature of the mutation, surviving homozygotes display a broad range of clinical manifestations including defective development of embryonic ectoderm and ectodermal derivatives. Defective pectoral fins in zebrafish [START_REF] Willer | Analysis of the Zebrafish perplexed mutation reveals tissuespecific roles for de novo pyrimidine synthesis during development[END_REF], nerve cells in humans [START_REF] Ng | Biallelic mutations in CAD, impair de novo pyrimidine biosynthesis and decrease glycosylation precursors[END_REF][START_REF] Koch | CAD mutations and uridine-responsive epileptic encephalopathy[END_REF], and wings in fruit flies [START_REF] Norby | A specific nutritional requirement for pyrimidines in rudimentary mutants of Drosophila melanogaster[END_REF][START_REF] Norby | The biochemical genetics of rudimentary mutants of Drosophila melanogaster. I. Aspartate carbamoyltransferase levels in complementing and non-complementing strains[END_REF][START_REF] Falk | Sex-linked auxotrophic and putative auxotrophic mutants of Drosophila melanogaster[END_REF], are among such examples. The poor udder morphology observed in the homozygote cow could be added to this list since the mammary gland is also of ectodermal origin [START_REF] Hens | Key stages of mammary gland development: molecular mechanisms involved in the formation of the embryonic mammary gland[END_REF][START_REF] Cowin | Molecular mechanisms guiding embryonic mammary gland development[END_REF][START_REF] Macias | Mammary gland development[END_REF].

To test whether the CAD mutation had an effect on the development of the mammary gland, we analyzed the performances of the daughters of 250 heterozygous carriers and 2,912 non-carriers for four traits (teat direction, teat placement, udder cleft, and udder depth). We did not find significant difference between DYD of heterozygous vs wild homozygous bulls (Figure 2.6). In addition, there was no reported udder related QTL in Normande cattle, either within NH7 interval or on chromosome 11 [START_REF] Boichard | Detection of genes influencing economic traits in three French dairy cattle breeds[END_REF], Marete et al., 2018). This is in line with a possible recessive effect of the mutation on udder morphology in surviving homozygotes, although more live homozygous cows are needed to confirm this assumption.

Conclusion

Here, we perform a systematic screening for recessive lethal mutations from large-scale haplotype scan to single-variant resolution mapping in Normande cattle. We report the identification of a deleterious substitution in CAD (p.Tyr452Cys) responsible for embryonic lethality with incomplete penetrance. We illustrate the influence of sample size, age of haplotype, and correction for multiple testing, on discovery and replicability of HHD results. Furthermore, we highlight the importance of large-scale genotyping of candidate variants, and phenotyping of homozygotes, to validate or invalidate them. Though such approaches require collection of genomic data on hundreds of thousands of animals over several years, thanks to genomic selection, this is a routine practice in dairy cattle breeding. Finally, we added a test for CAD (p.Tyr452Cys) to the EuroG10K BeadChip, which is routinely used for genomic evaluation in France. We recommend using CAD information to avoid at-risk mating. Because the economic weight of embryonic lethal defects is limited, we also suggest including it in the breeding objective as proposed by others, instead of eliminating all carriers.

Author Contributions: MMU, DB and AC conceived the study. MMU, CH, AB, PM, RL, MB, SF, DB and AC contributed in collection and/or analysis of the data. MMU drafted the manuscript. GS, DB and AC critically revised the manuscript. AC coordinated the study. All authors read and approved the final manuscript.

Furthermore, a ~525 KB deletion on chromosome 23 is reported to be associated with stillbirth in Nordic Red Cattle (Sahana et al., 2016).

Earlier SV studies on cattle were mostly SNP-array based, such as array-comparative genomic hybridization [START_REF] Liu | Analysis of copy number variations among diverse cattle breeds[END_REF], 50K BovineSNP50 BeadChip [START_REF] Hou | Genomic characteristics of cattle copy number variations[END_REF] or 777K BovineHD BeadChip (BovineHD chip) [START_REF] Xu | Population-genetic properties of differentiated copy number variations in cattle[END_REF] based. But, using these approaches a substantial portion of the genome could not be explored and breakpoint resolution is still an issue (Alkan et al., 2011). However, whole-genome sequence (WGS) based techniques could improve resolution as well as power to capture SVs in a wide size and frequency spectrum (Alkan et al., 2011). For example, majority of the novel SVs in humans [START_REF] Mills | Mapping copy number variation by populationscale genome sequencing[END_REF][START_REF] Sudmant | An integrated map of structural variation in 2,504 human genomes[END_REF] and mouse [START_REF] Yalcin | Sequence-based characterization of structural variation in the mouse genome[END_REF] were detected using WGS approaches. Besides, breakpoint sequences could also be assembled with high accuracy from sequencing reads [START_REF] Chen | TIGRA: a targeted iterative graph routing assembler for breakpoint assembly[END_REF], which are necessary for elucidating the mechanisms underlying SV formation [START_REF] Carvalho | Mechanisms underlying structural variant formation in genomic disorders[END_REF].

In the advent of next-generation sequencing (NGS) techniques, hundreds of cattle (bull or cow) genomes were re-sequenced in collaborative initiative such as 1000 Bull Genomes Project (1KBGP) (Daetwyler et al., 2014) (and other independent projects, e.g. [START_REF] Jansen | Assessment of the genomic variation in a cattle population by re-sequencing of key animals at low to medium coverage[END_REF], Brondum et al., 2014)) to build a comprehensive database of sequence variants, mainly SNPs and indels. This NGS data provides a unique opportunity to study SVs in cattle. However, few studies (Boussaha et al., 2015, Chen et al., 2017) utilized these (and/or other) NGS resources so far for studying SVs in cattle.

Therefore, in this study we scanned the whole genome sequences of 175 cattle from three dairy breeds, namely Holstein, Jersey, and Nordic Red Cattle, to discover large deletions segregating in the population, and analyze their population-genetic properties. In particular, we focused on understating the population diversity, stratification, and plausible functional effects. We also explored the probable mechanisms of SV formation for a set of breakpoint-resolved deletions.

Materials and methods

Animal samples and ethics

This study was performed on whole-genome sequences (WGS) of 175 dairy cattle from three breeds, e.g. 67

Holstein, 27 Jersey, and 81 Nordic Red Cattle. The sample included 7 Holstein cows and 168 bulls from these three breeds -144 animals from Run 5 of 1K bull genomes project (1KBGP) (Daetwyler et al., 2014) and 31 animals from Nordic sequence data [START_REF] Brondum | Strategies for imputation to whole genome sequence using a single or multi-breed reference population in cattle[END_REF]. Genome sequences were generated using Illumina paired-end sequencing to an average coverage of 10-fold.

Here, we did not include any experimentation on animals and only dealt with analysis-ready WGS data; hence, no ethical approval was required.

following: 1) an initial denaturation at 95

Analysis of population genetic properties

The population genetic properties of deletions, among the three breeds, were studied in terms of population diversity, population structure, and population differentiation. Population diversity was calculated using "VariantsPerSampleAnnotator" from Genome STRiP software, which provides distribution of variants across samples and populations. We performed principal component analysis (PCA) using PLINK (v1.90p) software [START_REF] Purcell | PLINK: a tool set for whole-genome association and population-based linkage analyses[END_REF] to distinguish three cattle breeds (details in Supplementary data). We calculated VST [START_REF] Redon | Global variation in copy number in the human genome[END_REF]) -a population stratification measure of structural variants (highly correlated with Wright's fixation index, FST [START_REF] Wright | Evolution in Mendelian Populations[END_REF]), for each deletion locus using variant allele frequency (VAF) and genotypes from pairwise comparison of one breed with the rest, such as Holstein vs Jersey+Nordic Red Cattle, and vice versa.

Functional annotation and enrichment analysis

Functional annotation of deletions were performed using Variant Effect Predictor (VEP-87) software (McLaren et al., 2016), and enrichment of protein domains (InterPro [START_REF] Finn | InterPro in 2017-beyond protein family and domain annotations[END_REF] and Pfam [START_REF] Finn | The Pfam protein families database: towards a more sustainable future[END_REF]) and pathways (KEGG [START_REF] Kanehisa | KEGG: new perspectives on genomes, pathways, diseases and drugs[END_REF]) were analyzed using STRING-v10 database [START_REF] Szklarczyk | STRING v10: protein-protein interaction networks, integrated over the tree of life[END_REF].

Selective constraints on genes were measured from the ratio of nonsynonymous (dN) to synonymous (dS) substitution rate, i.e. dN/dS ratio, between cow-mouse 1-to-1 orthologs downloaded from Ensembl database [START_REF] Yates | Ensembl 2016[END_REF] (release 87, last accessed on February 21, 2017) using BioMart (Kinsella et al., 2011).

Here we analyzed whether dN/dS of genes overlapping deletions are higher (i.e. less constrained) than that of mouse lethal genes (from

Data manipulation, visualization and statistical analysis

All statistical analyses and plots were generated in RStudio software (RStudio Team, 2016) running R software version 3.3.2 (R Core Team, 2016), unless mentioned otherwise. BEDTools (v2.26.0) software (Quinlan and Hall, 2010) is used for identifying the overlap between deletion calls and other genomic features, such as, UMD3.1 assembly gaps (from UCSC Table Browser), CNVs from dbVar database [START_REF] Lappalainen | DbVar and DGVa: public archives for genomic structural variation[END_REF], three known deletions from (Charlier et al., 2012, Kadri et al., 2014, Sahana et al., 2016), QTL from QTLdb. VCFtools (v0.1.15) software and PLINK (v1.90p) software were used for analyzing the VCF file.

Results and discussion

Discovery and genotyping of deletions

Deletion discovery and genotyping were carried out using Genome STRiP. After filtering, we report 8,480 large deletions with genotypes in 67 Holsten, 27 Jersey, and 81 Nordic Red Cattle. The deletion size ranged from 199 bp to 773 KB with a mean of 4.5 KB (median = 1 KB), which is approximately 10 times smaller compared to 184 deletion-CNVs (mean = 44.5 KB, median = 7.7 KB) reported in a recent 777K BovineHD BeadChip (BovineHD chip) based study [START_REF] Xu | Population-genetic properties of differentiated copy number variations in cattle[END_REF], reflecting the resolution of our sequence-based calls. Only 18% of the deletion calls have overlap (≥1 base pair) with previously reported bovine deletions (or CNV-loss) in the dbVar database (accessed on January 27, 2017), while remaining 82% are novel.

However, ~72% of our deletion regions remained unique when compared to all CNVs (gain or loss) and copy number variable regions (CNVRs) in the database. Interestingly, majority (~80%) of these overlapping regions are from an earlier WGS-based study, where genome sequences of 27 Holstein, 17 Montbéliarde and 18 Normande bulls were analyzed (Boussaha et al., 2015). Nonetheless, we were able to broaden the accessible deletion size-range, more importantly towards smaller one unascertainable by usual SNP-array based approaches. We also report high quality genotypes for all the 8,480 deletions. Apparently, there are more low frequency variants than that of high frequency one, and the frequency distribution is very similar to humans [START_REF] Mills | Mapping copy number variation by populationscale genome sequencing[END_REF]) (Figure 3.1).

Previous NGS-based studies on cattle were mostly limited to SV discovery, while copy-number states (genotypes) were inferred using BovineHD chip (Chen et al., 2017), or custom SNPs array (Boussaha et al., 2015). However, in this study we estimated the copy-number at each deletion locus (per sample) from readdepth within the region using a constrained Gaussian mixture model with three classes, e.g. copy-number zero (i.e. homozygous deletion), one (i.e. heterozygous deletion) and two (i.e. homozygous reference). It is known from human studies that majority of the (bi-allelic) common SVs segregate on specific SNP haplotypes [START_REF] Mccarroll | Integrated detection and populationgenetic analysis of SNPs and copy number variation[END_REF], Conrad et al., 2010), which could be imputed with high accuracy (Handsaker et al., 2011, Handsaker et al., 2015). Thus, this approach has the potential, albeit with large reference, for accurate haplotype phasing and imputation of SVs to large cohorts of low-density chip-typed animals with no additional cost.

Validation of deletions

We validated the results using three approaches: (i) using BovineHD chip intensity data, (ii) breakpoint assembly and alignment, and (iii) PCR + sequencing of amplicons.

BovineHD chip intensity.

We used 777K BovineHD chip intensity data of 26 Holstein animals, both chiptyped and sequenced, to validate the deletion calls using Genome STRiP' IntensityRankSum test. We had partial power to investigate all deletions due to the sparsity of the array-probes (one probe per ~3.5 KB), and were underpowered to accurately verify small deletions (e.g. overlapping one-probe). Furthermore, we could only test a deletion for which at least one of the 26 samples had non-reference allele. Therefore, an estimate of false discovery rate (FDR) for the deletion call-set is provided here from the overall p-value distribution.

In this approach, we were able to interrogate ~8.3% of the total call, majority of which contain a single arrayprobe within the region (Table 3.1 and Table S1). We found that deletions overlapping only one array-probe had higher FDR (11.3%) compared to two or more probes. And finally, we showed that our deletion call-set had an overall FDR of 8.8%, which is within our chosen threshold of FDR ≤ 10%. Targeted breakpoint assembly. We next validated three randomly chosen set of ten-deletions each by assembling breakpoint sequences using TIGRA [START_REF] Chen | TIGRA: a targeted iterative graph routing assembler for breakpoint assembly[END_REF]: 10 deletions ≤ 500 bp, 10 deletions > 500 bp but ≤1KB, and 10 deletions with VAF ≤ 0.10. Out of the thirty, we successfully resolved breakpoints of 26 deletions (~87% success rate) using a combination of TIGRA and BLAT search [START_REF] Kent | BLAT--the BLAST-like alignment tool[END_REF] (Table S2). Additionally, we assembled breakpoints of three previously reported deletions that were also segregating in our study populations, such as a ~662 KB deletion on chromosome 12 encompassing RNASEH2B, GUCY1B2, and FAM124A, a ~3 KB deletion on chromosome 21 encompassing FANCI, and a ~525 KB deletion on chromosome 23 encompassing BTBD9, GLO1, and DNAH8. Breakpoints of the former two deletions were previously reported (Charlier et al., 2012, Kadri et al., 2014), which exactly matched with our predicted breakpoint sequences (Figures S1abc andS2abc). While for the later, we resolved breakpoint sequences in this study (Figure 3.2ab). Overall, the success rate of our deletion-breakpoint assembly was better than the reported success rate of TIGRA on similar sized read-length [START_REF] Chen | TIGRA: a targeted iterative graph routing assembler for breakpoint assembly[END_REF].

And Genome STRiP's breakpoint predictions were on average within 20 bp of the validated breakpoint, which is within the tool's reported estimate of 1-20 bp (Handsaker et al., 2011). PCR and amplicon sequencing. We then experimentally validated breakpoints for "Chr23:12,291,761-12,817,087" deletion, previously reported to be associated with stillbirth in Nordic Red Cattle (Sahana et al., 53 2016). Four animals were used for PCR validation: 2 carriers and 2 non-carriers. The breakpoint spanning PCR products of 359 bp were only observed in the carrier animals, while no amplicon was seen for noncarriers (Figure 3.3a). The 359 bp amplicon was then sequenced, and exact breakpoint sequences were observed (Figure 3.3b), thus confirming the breakpoint for this deletion.

Population genetic properties of deletions

Population diversity. We explored population diversity among the three dairy cattle breeds from perindividual deletion-heterozygosity and homozygosity. We found that individuals from Nordic Red Cattle exhibits 3.5% and 6.4% higher deletion-heterozygosity than in Holstein and Jersey, respectively. Median numbers of heterozygote-deletion were 1272, 1229 and 1196, in Nordic Red Cattle, Holstein, and Jersey, Chapter 3 54 respectively (Figure 3.4a). In contrast to heterozygosity, Jersey animals showed highest levels of deletionhomozygosity, followed by Holstein (Figure 3.4b). Similar estimates of genetic diversity were also reported for these breeds in a SNP heterozygosity and runs-of-homozygosity analysis -where Jersey exhibited lowest (genome-wide) average nucleotide diversity (and higher number/size of runs-of-homozygosity) followed by Holstein and Nordic Red Cattle [START_REF] Zhang | In each window, 600 conditioning haplotypes were generated, among which, 400 were Hamming distance conditioning (--states 400) and the remaining 200 were chosen randomly (--states-random 200) to improve the mixing of the MCMC. However, both the 50k and 777k dataset were already phased using Beagle. We took the Beagle haplotypes and GPs, and performed phasing (and calling)[END_REF]. These differences could be understood from the current effective population size (Ne) of these breeds, e.g. Ne of Jersey, Holstein, and Nordic Red Cattle are 73, 99 and 106, respectively [START_REF] Bovine Hapmap | Genomewide survey of SNP variation uncovers the genetic structure of cattle breeds[END_REF]; this entailed higher diversity in Nordic Red Cattle, and Holstein, than in Jersey. From singletons estimate it is also evident that Nordic Red Cattle has more rare deletions compared to Holstein and Jersey (Figure 3.4c); partly could be due to incorrect ascertainment of rare ones. Principal Component Analysis (PCA). We performed Principal Component Analysis (PCA) using the deletion genotypes of the samples. Around 6K deletions with VAF between 0.02 and 0.90 were used in the analysis. For comparison, we also performed PCA on ~168K bi-allelic SNPs randomly selected from 29 autosomes of the same individuals. The first two principal components (PCs) from both deletion and SNPbased PCA clearly distinguished the three breeds, and jointly explained 20% and 16.2% of the variance, respectively (Figure 3.5ab). In addition, PC3 and PC4 recapitulated substructures within Nordic Red cattle (Figure S3), and first five PCs cumulatively explained 33.6% (with the deletions) and 28.4% (with the SNPs) of the variance (Figure S4). Our deletion results agree with the known population structure of the three gains and losses in mammalian evolution [START_REF] Niimura | Extensive gains and losses of olfactory receptor genes in mammalian evolution[END_REF]. And population specific copy-number variations of olfactory receptor genes were also reported in human (deletions) [START_REF] Van Ziffle | Homozygous deletion of six olfactory receptor genes in a subset of individuals with Beta-thalassemia[END_REF] and cattle (gains) [START_REF] Lee | Analysis of cattle olfactory subgenome: the first detail study on the characteristics of the complete olfactory receptor repertoire of a ruminant[END_REF]. Nevertheless, this is the first report, to our knowledge, of homozygous deletion of olfactory receptor genes in cattle.

Figure 3.7. Difference between dN/dS ratios of mouse-lethal and deletion-overlapped genes in cattle. Cow genes for which

one-to-one mouse orthologs available were considered for a one-sided Wilcoxon rank-sum test. Mouse lethal genes are from [START_REF] Dickinson | High-throughput discovery of novel developmental phenotypes[END_REF].

QTL Enrichment. We next explored the enrichment (or depletion) of quantitative trait loci (QTL) on deleted regions (at least 1 bp overlap with deletion). We retrieved ~24K autosomal QTL from QTLdb reported to be associated with any of the six trait classes, e.g. "Health", "Reproduction", "Milk", "Exterior", "Production" and "Meat and Carcass". The association of deletions with diseases, fitness or fertility related traits is well evident [START_REF] Weischenfeldt | Phenotypic impact of genomic structural variation: insights from and for human disease[END_REF]. Hence, we suspected enrichment of fitness and fertility related traits for our deletions. As expected, health (2 fold) and reproduction (1.5 fold) related QTL were significantly enriched, while other trait classes were highly depleted (Table 3.2). Higher enrichment of health related QTL could be driven by immune-system genes, which were also highly enriched in our dataset (discussed earlier). [START_REF] Hu | Developmental progress and current status of the Animal QTLdb[END_REF]. QTL from autosomes of Holsteins, Jersey, Nordic Red Cattle, and Ayrshires were considered for Fisher's exact test (two-sided).

Deletion formation mechanisms

Finally we explored the probable mechanisms of deletion formation. There are two key mechanisms of structural variants formation (for detail see review (Hastings et al., 2009, Carvalho and[START_REF] Carvalho | Mechanisms underlying structural variant formation in genomic disorders[END_REF]); for example, recurrent SVs often result from non-allelic homologous recombination (NAHR) between large low-copy repeats (LCRs), and thus, contain extensive sequence homology provided by LCRs, such as segmental duplicates, at the flanking regions [START_REF] Carvalho | Mechanisms underlying structural variant formation in genomic disorders[END_REF]. In contrast, non-recurrent SVs often form either by microhomology-mediated end joining (MMEJ) or non-homologus end joining (NHEJ), which requires limited to no sequence homology, and thus could be characterized by microhomologies or simple blunt ends at the breakpoint junction [START_REF] Hastings | Mechanisms of change in gene copy number[END_REF].

Breakpoint information is crucial for understanding the mechanism, and therefore, we analyzed 29 breakpoint resolved deletions from our validation set. We found that 24 of 29 deletions contain microhomology ranging from 2-31 bp at the breakpoint, and two of which also contain insertions (S2 Table ).

In addition, 4 deletions exhibited non-reference insertion at breakpoint junctions, and one deletion with no apparent homology. However, the number of breakpoint sequences analyzed here were not a robust representation of our deletion call-set (less than 0.5% deletions), though selected randomly (for validation), we were able to demonstrate that majority of deletions contain microhomoloy at breakpoint, followed by few insertions, and rarely with no homology. Our results largely agree with the trend reported for large deletions in humans, e.g. 70.8% deletions exhibited microhomology/homology and 16.1% insertions at the breakpoint [START_REF] Mills | Mapping copy number variation by populationscale genome sequencing[END_REF].

Limitations

This study only focused on identifying deletions in cattle because of their potential relevance to loss-offunction and embryonic lethality. However, we had limited success to identify small deletions, such as <200 bp due to reduced sensitivity of the SV caller. It is also not a comprehensive list of deletions for these samples, since we could have missed many true deletions due to sensitivity, coverage, or stringent filtering (among other reasons). Furthermore, the short read length (~100 bp) in our WGS dataset also made it difficult to resolve breakpoints from regions of long repeats.

Conclusions

Loss-of-function variants are responsible for a substantial yearly-economic loss in dairy industry, where a limited number of elite sires are in extensive use for rapid genetic gains. Mapping of such variants is essential for effective breeding planning and genomic selection. Here we showed an NGS-based analytical framework suitable for population-scale mapping of large deletions in cattle, leveraging the available WGSs. Here we described population-genetic, functional, and evolutionary properties of discovered deletions. We identified and confirmed a ~525 KB deletion on chromosome 23, causing stillbirth in Nordic Red Cattle. We demonstrated that Nordic Red Cattle had higher population diversity than Holstein and Jersey, and deletiongenotype could recapitulate genetic structure of these breeds. Natural gene knockouts are enriched for immune-related and olfactory receptor genes. We also showed that deletions are significantly enriched for health and fertility related QTL, while depleted for production related QTL. Our population genetic and functional analysis showed promise for inclusion of SVs in genomic studies in dairy cattle. This deletion catalog will facilitate discovery, genotyping, and imputation of deletions in large cohorts of animals, and subsequent studies for gene mapping and genomic prediction of breeding values.

Estimating Deletion Genotype Likelihoods from VCF file

We estimated genotype likelihoods for the selected deletions on 597 additional animals from the 1KBGP (Run-6) using RD data from VCF file.

Calculation of expected RD. Studies showed that GC bias, i.e., RD variation due to differences in the percentage of guanine (G) and cytosine (C) bases in the region, can confound the inference of copy-number from RD (Abyzov et al., 2011, Benjamini and[START_REF] Benjamini | Summarizing and correcting the GC content bias in high-throughput sequencing[END_REF]. For the genome assembly UMD3.1, we first calculated the GC% in bins of 100 bp after excluding potential copy-number variable regions, such as sex chromosomes, mitochondrial sequences, unplaced contigs, assembly gaps, repeat sequences, and CNVs and SVs from DGVa database (last accessed on January 22, 2019, from: https://www.ebi.ac.uk/dgva). Next, for each animal, we extracted RD data from the corresponding genomic intervals from VCF files, and calculated the expected RD (𝜇𝜇 𝐺𝐺𝐺𝐺 ) and variance (𝑣𝑣𝑣𝑣𝑣𝑣 𝐺𝐺𝐺𝐺 ) for different GC% (for details, see Supplementary Methods).

There was substantial variations in RD against different GC% between samples (Supplementary Figure S1).

To account for this heterogeneous GC% bias for a given deletion locus, the expected RD was assumed to be the average RDs from its GC% bin, instead of the overall genomic coverage. For example, when interrogating a 10 Kb deletion locus with 30% GC content, we used the genomic average over all 30% GC bins as the expected RD for each corresponding animal.

Gaussian mixture model. To compute genotype likelihoods for a deletion, we first extracted the SNPs (say n SNPs in total) within the deleted segment and retrieved read depth data from the allelic-depth (AD) tag of the VCF file. Here, we only considered bi-allelic SNPs with QUAL ≥30 and a thinning of 10 bp. Next, we fitted a Gaussian mixture model (GMM) to the data, assuming a linear relationship between the observed RD within the locus and (unobserved) copy-number (CN) of the locus, such as,

𝑔𝑔(𝑅𝑅𝐶𝐶 𝑖𝑖 ) = � 𝑤𝑤 𝑖𝑖 𝑁𝑁(𝑅𝑅𝐶𝐶 𝑖𝑖 | 𝑚𝑚𝑍𝑍 𝑖𝑖 , 𝑣𝑣𝑣𝑣𝑣𝑣 𝑖𝑖 ) 2 𝑖𝑖=0
Here, RD is a vector of n data points corresponding to the n SNPs; k (0 to 2) is an indicator variable for the Gaussian component; 𝑤𝑤 𝑖𝑖 is the relative weight for CN classes (e.g., 𝑤𝑤 0 , 𝑤𝑤 1 , 𝑤𝑤 2 for CN0, CN1 and CN2, respectively); 𝑚𝑚𝑍𝑍 𝑖𝑖 and 𝑣𝑣𝑣𝑣𝑣𝑣 𝑖𝑖 are the expectation and variance of the corresponding Gaussian distribution.

The GMM parameters were estimated using an expectation-maximization (EM) algorithm as follows:

Initialization of the EM. Assuming a deletion locus, we constrained GMM to fit exactly three components with a fixed 𝑚𝑚𝑍𝑍 𝑖𝑖 for each k such that 𝑚𝑚𝑍𝑍 0 = 0 for CN0, 𝑚𝑚𝑍𝑍 1 = 𝜇𝜇 𝐺𝐺𝐺𝐺 2 for CN1, and 𝑚𝑚𝑍𝑍 2 = 𝜇𝜇 𝐺𝐺𝐺𝐺 for CN2. We initialized 𝑤𝑤 𝑖𝑖 with three equal weights of 1/3, and 𝑣𝑣𝑣𝑣𝑣𝑣 𝑖𝑖 with the value of 0.2, 𝑔𝑔𝑎𝑎𝑖𝑖 𝐺𝐺𝐺𝐺 2 , and 𝑣𝑣𝑣𝑣𝑣𝑣 𝐺𝐺𝐺𝐺 , for CN0, CN1 and CN2, respectively. E-step. Given the parameters, the expected values for the latent variable 𝑍𝑍 𝑖𝑖𝑖𝑖 were estimated as:

𝑍𝑍 𝑖𝑖𝑖𝑖 = 𝑤𝑤 𝑖𝑖 𝑁𝑁(𝑅𝑅𝐶𝐶 𝑖𝑖 | 𝑚𝑚𝑍𝑍 𝑖𝑖 , 𝑣𝑣𝑣𝑣𝑣𝑣 𝑖𝑖 ) ∑ 𝑤𝑤 𝑖𝑖 𝑁𝑁(𝑅𝑅𝐶𝐶 𝑖𝑖 | 𝑚𝑚𝑍𝑍 𝑖𝑖 , 𝑣𝑣𝑣𝑣𝑣𝑣 𝑖𝑖 ) 2 𝑖𝑖=0
M-step. The membership weights and variances were updated using the 𝑍𝑍 𝑖𝑖𝑖𝑖 values from the E-step (while keeping 𝑚𝑚𝑍𝑍 𝑖𝑖 fixed):

𝑤𝑤 𝑖𝑖 = 1 𝑚𝑚 � 𝑍𝑍 𝑖𝑖𝑖𝑖 𝑎𝑎 𝑖𝑖=1 𝑣𝑣𝑣𝑣𝑣𝑣 𝑖𝑖 = ∑ 𝑍𝑍 𝑖𝑖𝑖𝑖 (𝑅𝑅𝐶𝐶 𝑖𝑖 -𝑚𝑚𝑍𝑍 𝑖𝑖 ) 2 𝑎𝑎 𝑖𝑖=1 ∑ 𝑍𝑍 𝑖𝑖𝑖𝑖 𝑎𝑎 𝑖𝑖=1
For a given deletion, parameters 𝑤𝑤 𝑖𝑖 and 𝑣𝑣𝑣𝑣𝑣𝑣 𝑖𝑖 were estimated iteratively from the observed RD data of each animal using the EM algorithm until convergence. The scripts used to prepare the RD data from the VCF file and to estimate deletion PL could be accessed from:

"https://github.com/MMesbahU/ImputeDelPipeline/tree/master/read_Depth_genotyping".

Refining the genotypes. Next, to refine the deletion genotypes along with the WGS SNPs and indels, we combined PL estimates for all these variants and performed chromosome-wise phasing using Beagle v4.r1274 software (Browning and Browning, 2016) (hereafter referred to as Beagle; for details see "Phasing").

Evaluation of the deletion genotyping pipeline. We evaluated the pipeline for Brachyspina-associated deletion located on chromosome 21 between the positions 21,184,869 and 21,188,202 bp (Charlier et al., 2012). In this dataset, we had 113 Holstein with known Brachyspina carrier-status. There were 13 carriers and 100 non-carriers of the deletion. Using the GMM approach described before, we first estimated deletion PL using the SNP RD from the VCF file. After excluding variants within the deletion, we extracted PL for WGS SNPs and indels within ±1 Mb of the deletion. Next, to refine the deletion genotypes, we performed phasing using the surrounding haplotypes using Beagle (with same parameters as discussed in the "Phasing" step). Finally, the genotype concordance (true vs. Beagle's best guess genotypes) was manually inspected.

Phasing

Following the method used by Delaneau et al. (2014), we phased our genotype dataset using a combination of software Beagle and SHAPEIT v2.r837 (Delaneau et al., 2013a) (hereafter referred to as SHAPEIT). First, providing the WGS PLs as input, we ran Beagle for 10 burn-in and 15 sampling iterations with a window size of 12,000 markers and an overlap of 2,000 markers between consecutive windows. From this step, we obtained posterior genotype probabilities (GP) for each bi-allelic variant for the three possible genotypes, e.g., homozygous reference, heterozygous (or hemizygous for deletion), and homozygous for alternate allele (or homozygous deletion). We fixed the genotype for a variant as known when Beagle GP ≥0.99, and for the remaining variants, we re-called the genotype using SHAPEIT. To initialize phasing and calling in SHAPEIT, we used the haplotypes generated in the Beagle step. After this initialization, SHAPEIT was run for 12 pruning stages of four iterations each, followed by 20 main MCMC sampling iterations. We used a window size of 0. 175 animals. For the remaining 597 animals, we estimated PL by fitting the GMM with three components to the RD data extracted from the VCF file. We then used surrounding haplotypes to refine deletion genotypes in animals with ambiguous (or missing) PL estimates.

We evaluated the deletion genotyping pipeline for classifying animals for their carrier status for the Brachyspina-associated deletion (Charlier et al., 2012) in a panel of 113 animals with known carrier status including 13 hemizygotes and 100 homozygotes for the reference allele. We found that our two-step approach-estimating PL using WGS SNP RD data from the VCF file followed by a genotype refinement step using Beagle-accurately classified deletion genotypes in all the animals (13 hemizygotes and 100 homozygotes for the reference allele).

The WGS Reference Panel of 772 Animals

Following previous reports [START_REF] Brøndum | Strategies for imputation to whole genome sequence using a single or multi-breed reference population in cattle[END_REF][START_REF] Pausch | Evaluation of the accuracy of imputed sequence variant genotypes and their utility for causal variant detection in cattle[END_REF], where the authors showed that a multi-breed reference panel provides better downstream imputation performance in cattle compared with a single-breed reference panel, we built a multi-breed WGS reference panel of 772 animals that included 92 JER, 132 RDC and 548 HOL. We applied the GMM pipeline to estimate the genotype likelihoods for all the deletions in the additional animals included from the 1KBGP. Next, to refine the genotypes and get posterior genotype probability (GP), we performed phasing by chromosome using Beagle with the PL estimates for all SNPs, indels and deletions as inputs.

However, studies in humans have shown that an initial phasing of WGS genotypes using Beagle followed by fixing the genotype with high posterior probability (e.g., Beagle GP >0.995) and re-calling the remaining markers genotypes using SHAPEIT "-call" provides further improvement in the reference panel (Delaneau et al., 2014, Genomes Project et al., 2015). Besides, SHAPEIT can also improve haplotype phasing by leveraging on the haplotype sharing within and between populations with heterogeneous ancestry (Delaneau et al., 2013b). Here, both our WGS reference and SNP array-typed target populations consist of animals with diverse ancestry (see Fig. 3 from Bouwman et al. (2018) for examples of haplotype diversity in JER, HOL, and RDC; see [START_REF] Zhang | Genome-wide association studies and genomic prediction of breeding values for calving performance and body conformation traits in Holstein cattle[END_REF] for a detailed account in RDC). In addition, on average ~2.7% of the total WGS genotypes in our dataset had Beagle GP <0.99 (Table S1). Furthermore, there were 7.3% (1.1 of 14.7 million) SNPs, 18.3% (22.4 of 122.1 k) indels and 66.6% (3.9 of 5.8 k) deletions where more than 10% genotypes were below this threshold (Table S1). Therefore, we expected that by using SHAPEIT we could further improve the reference panel by improving (low-quality) genotype calls and haplotype phasing. Based on this assumption, we next refined the reference panel using SHAPEIT, where all genotypes were rephased, and re-called when Beagle GP was less than 0.99.

To evaluate the effect of this two-step refinement approach on imputation accuracy we prepared RefBeagle and RefBeagle_SHAPEIT of the chromosome 29. First, we performed leave-out-trials on 112 deletions with MAC ≥5 from this chromosome and evaluated the imputation accuracy using (squared) Pearson's correlation (r 2 ) between true vs. imputed dosages of 0, 1, or 2 for homozygous reference, hemizygous or homozygous deletion, respectively (see Methods "Cross validation"). After the trials, five deletions became homozygous for the reference allele. These were excluded from subsequent analysis. Overall, the average gain in accuracy was r 2 = 0.081 (95% CI: 0.053, 0.109) with RefBeagle_SHAPEIT vs. RefBeagle (P = 7.9 × 10 -8 ; two-tailed paired t-test). The performance gain was even higher (r 2 = 0.105 [CI: 0.069, 0.141], P = 1.7 × 10 -7 ) for variants where more than 10% of the genotypes were under the GP threshold (Table 4.1). Furthermore, 76 vs. 90% of the deletions achieved r 2 >0.80 using the reference RefBeagle vs. RefBeagle_SHAPEIT, respectively (Figure 4.1a).

Next, we evaluated the downstream performances of the two reference panels to impute jointly a target population of ~13k animals with 777k genotypes to WGS SNPs, indels and deletions of chromosome 29 using Minimac. As expected, the reference panel RefBeagle_SHAPEIT showed better downstream imputation performance for all variant types compared to that of RefBeagle (Figure 4.1b, 4.2ab, and Table 4.1). The average gains in accuracy in terms of Minimac r 2 were 0.032 (CI: 0.012, 0.053; P = 0.002) for deletions, 0.020 (CI: 0.019, 0.020; P <1.0 × 10 -16 ) for SNPs and 0.034 (CI: 0.031, 0.037; P <1.0 × 10 -16 ) for indels (Table 4.1). However, after Beagle step, ~8% of the SNPs, 20% of the indels and 70% of the deletions had more than 10% of genotypes below the GP threshold (Table S1). Indeed, for these markers we observed a clear difference between before vs. after SHAPEIT step. For example, the average gain in accuracy after the SHAPEIT step was 0.041 (CI: 0.011, 0.070; P = 0.007) for deletions, 0.058 (CI: 0.056, 0.059; P <1.0 × 10 -

16

) for SNPs and 0.076 (CI: 0.065, 0.087; P <1.0 × 10 -16 ) for indels (Table 4.1). Besides, 87 vs. 93% of the deletions (Figure 4.1b), 87 vs. 89% of the SNPs (Figure 4.2a), and 78 vs. 83% of the indels (Figure 4.2b) achieved r 2 >0.80 using the reference panel RefBeagle vs. RefBeagle_SHAPEIT, respectively. It is worth noting that for the majority of the imputed variants, allele frequencies differed between the two reference panels.

Nevertheless, the average r 2 in different 5% MAF bins were also higher with RefBeagle_SHAPEIT vs. RefBeagle (Figure S2 and Table S2). Thus, it is likely that SHAPEIT "-call" step improved genotype calls, providing an improved WGS reference panel for the chromosome. These results are in agreement with the reported performance gains in humans (Delaneau et al., 2014). reference panel phased using Beagle followed by re-phasing using SHAPEIT.

Imputation of WGS SNPs, Indels and Deletions to the Chip-Typed Cattle Population

Our imputation target was a population of ~13k animals (1,640 JER, 6,371 HOL and 4,955 RDC) that were genotyped using either the 50k or the 777k BeadChip. We prepared this dataset using the same pipeline: two-step phasing followed by imputation from the 50k to the 777k genotypes using Minimac (see Methods "Phasing" and "Imputation"). Next, for the 29 bovine autosomes, we jointly imputed this target 777k dataset to WGS SNPs, indels and deletions using the reference panel RefBeagle_SHAPEIT comprised of 772 animals. In Figure 4.4a-c, we present the MAF vs. Minimac r 2 estimates for the imputed WGS markers. The average imputation accuracy for SNPs, indels and deletions were 0.928 (SD 0.207), 0.900 (0.217) and 0.904 (0.220), respectively. Here, irrespective of the marker type, the common variants (MAF >5%) were well imputed with an average r 2 >0.90 (Table 4.2). Furthermore, low frequency variants (1%≤ MAF ≤5%) were also well imputed with an average r 2 value of 0.889 (0.250) for SNPs, 0.874 (0.245) for indels and 0.927 (0.154) for deletions. As for the rare variants (MAF <1%), we observed a sharp reduction in imputation accuracy compared to that of low frequency variants: differences of 0.361, 0.325 and 0.214 were seen for SNPs, indels and deletions, respectively (Table 4.2). Similar low imputation accuracies for rare variants were also reported previously [START_REF] Brøndum | Strategies for imputation to whole genome sequence using a single or multi-breed reference population in cattle[END_REF][START_REF] Pausch | Evaluation of the accuracy of imputed sequence variant genotypes and their utility for causal variant detection in cattle[END_REF]. Rare variants are often novel in the population, have low linkage-disequilibrium with neighboring markers and poorly represented in the reference panel, if at all.

For these variants, we further stratified the r 2 into MAF bins of 0.1% (Table 4.2). We found that both the SNPs and indels with MAF >0.5% had an imputation accuracy r 2 >0.50. The reduction was mainly due to variants in the bottom half (in terms of MAF) of this category. In contrast to SNPs and indels, rare deletions with MAF = 0.2% were also imputed with high accuracy (r 2 = 0.633), and for the lowest 0.1% MAF bin (i.e.

0< MAF ≤0.1%), r 2 was 0.313 (Table 4.2). Nevertheless, compared to the previous report on this population (see Figure S2 from Wu et al. ( 2016)), we observed substantial gains in accuracy in a wide range of MAF bins (Figure 4.3, Table 4.2 andS3). These results also showed that, given a suitable reference panel, Minimac could jointly impute various classes of bi-allelic sequence variants-in a frequency spectrum range from rare to common-to the existing chip-typed population with high accuracy.

Principal Component Analysis

Next, to assess the usability of this dataset in population-genetic studies, we performed principal component analysis using a random sample of the imputed variants with MAF >1%. Here we verified that the three classes of imputed markers were able to distinguish the three population present in our dataset in a similar way. Indeed, with the top two principal components from SNPs, indels and deletions, we were able to illustrate the population structure present in our dataset (Figure 4.4).

Conclusion

In this study, we presented strategies to incorporate large chromosomal deletions in population genetic studies, along with SNPs and indels. First, we showed an approach to estimate deletion genotype likelihood using the read depth data from the VCF file. Second, using a two-step genotype refinement, we built WGS reference panels that included SNPs, indels, and deletions. Finally, we imputed our 777k population to full sequence. We found that taking genotype uncertainty in the reference panel into account leads to substantial gains in accuracy. These results highlighted the feasibility of building a WGS reference panel comprising different classes of sequence variants and jointly imputing them to the existing SNP chip-typed large cohorts of animals. In addition, we showed that, as WGS resources, such as the 1KBGP, both include animals with shallow and deep sequencing depth, refining the low-confidence genotypes yielded a substantial gain in imputation accuracy. Furthermore, such an approach also provided accurate genotypes for rare variantswhich often include putative causal variants. This, in turn, is expected to improve downstream association signals.

This study was focused only on the imputations of deletions along with SNPs and indels. With such an accurate imputation procedure even for rare and low frequency variants, it becomes possible to explore their potential biological effects by genome-wide association study and to study their predictive ability in genomic evaluation. This requires further investigation and will be presented in another study. It is likely that the strategies and pipelines demonstrated here will facilitate joint imputation of CNVs, SVs, SNPs and indels, and subsequent population genetic studies in cattle. Nevertheless, this imputation resource will be used in association mapping and estimation of (genomic) breeding values in Nordic dairy cattle breeds.

. Imputation accuracy of WGS SNPs, indels and deletions on the 29 bovine autosomes. Each dot represents an

imputed marker with MAF >0. There are (a) 14,070,960 SNPs, (b) 122,054 indels and (c) 5,730 deletions. Here, the imputation accuracy were greater than 0.80 for 91% SNPs, 85% indels and 87% deletions. The solid line in each scatterplot represents the average Minimac r 2 in MAF bins of 1% (see Table S3 for corresponding r 2 values with standard deviations). For the imputed variants, the Pearson's correlation coefficients, r, between the MAF in the reference panel vs. the target population were 0.98 for SNPs, 0.97 for indels and 0.98 for deletions. The observed average differences in MAF, between target vs. reference, were 0.0023 (two-tailed t-test P <1.0×10 -16 ) for SNPs, 0.0036 (P = 6.7×10 -12 ) for indels and 0.0009 (P = 0.73) for deletions. to improve GWAS signals even with a fixed size of the mapping population. The first one is to improve accuracy of imputation, and the second is to include other classes of genetic markers, besides SNPs and indels, such as chromosomal deletions, duplications, insertions, to name a few. In earlier studies, we attempted to address both factors and reported substantial gains in accuracy of imputation for WGS SNPs, indels and large chromosomal deletions (thereafter named deletions) in three Nordic dairy cattle breeds, namely, HOL, JER, and RDC (presented in Chapter 3 and 4 of this thesis). In this study, we performed GWAS for eight fertility traits in these three breeds using imputed WGS dataset comprise of SNP, indel and deletion.

Methods

Genotypes and phenotypes

We used the imputed WGS dataset comprised of SNPs, indels and deletions from Chapter 4, where we performed a two-step imputation from the Illumina BovineSNP50 to BovineHD BeadChip (777k, hereafter) genotypes, followed by joint imputation from the 777k genotypes to WGS SNPs, indels and deletions. The phasing was performed using a combination of Beagle v4.r1274 (Browning and Browning, 2016) and SHAPEIT v2.r837 (Delaneau et al., 2013) software following the method reported in (Delaneau et al., 2014), and joint imputation was performed using Minimac3 v.2.0.1 software (Das et al., 2016). Overall, average imputation accuracies (Minimac r 2 ) were 0.93 for SNPs, 0.90 for indels and 0.90 for deletions. From this imputation dataset (which had several filters for low quality markers including deviation from Hardy-Weinberg proportion, P-value ≤1 × 10 -5 ), we considered markers with minor allele frequency (MAF) >1% and Minimac r 2 value ≥0.1 for GWAS. In the final dataset, we had 5,596 HOL, 1,215 JER and 4,507 RDC bulls with both genotypes and phenotypes (de-regressed estimated breeding values). At this MAF and Minimac r 2 , the genotype dataset comprised 12,174,287 (12,076,004 SNPs;94,581 indels;3,702 deletions) markers in HOL,10,135,228 (10,057,293 SNPs;74,829 indels;3,106 deletions) markers in JER and 12,855,372 (12,748,285 SNPs,102,843 indels,4,244 deletions) markers in RDC. For each breed, the phenotype dataset comprised of de-regressed breeding values for eight fertility traits: namely, number of inseminations per conception in heifers (AISh) or in cows (AISc), interval (number of days) from calving to first insemination in cows (ICF), interval (number of days) from first to last insemination in heifers (IFLh) or in cows (IFLc), fertility index (FI), and non-return rate in heifers (NRRh) or in cows (NRRc) (for details see (NAV, 2013)). Summary of these phenotypes was presented in Table 5.1.

Genome-wide association study

We performed single-marker association analysis using following mixed linear model:

𝒚𝒚 = 𝟏𝟏 𝑎𝑎 𝜇𝜇 + 𝑏𝑏𝒙𝒙 + � 𝑐𝑐 𝑖𝑖 𝑷𝑷𝑷𝑷 𝒊𝒊 10 𝑖𝑖=1 + 𝒈𝒈 + 𝜺𝜺
where y is a vector of phenotypic records for fertility traits, 1n is a vector of ones, µ is the mean of the given trait, b is the additive allele substitution effect of the candidate marker to be tested, x is a vector of allele dosages (coded as 0, 1 or 2), PCi is a vector of the principal component i (first 10 PCs were considered) and ci is the corresponding coefficients, g is a vector of polygenic effect with 𝒈𝒈 ~𝑁𝑁�0, 𝑮𝑮𝜎𝜎 𝑔𝑔 2 �, where G is the genomic relationship matrix calculated using the 777 k markers without the test chromosome (leave-onechromosome-out (LOCO) approach), and 𝜺𝜺 is the vector of residual effect with 𝜺𝜺 ~𝑁𝑁(0, 𝑰𝑰𝜎𝜎 𝜀𝜀 2 ), where I is an identity matrix. Differences in reliabilities of de-regressed proofs were neglected because of the large number of daughters. The analysis was performed using GCTA software (version 1.92.1 beta6) (Yang et al., 2011a).

To declare an association with the trait, we used a genome-wide significant threshold of P-value <5 × 10 -8 , considering 1 million independent test at 5% level of significance.

We estimated the genomic inflation, lambda, for the single-marker GWAS P-values, using GenABEL software "estlamda" function [START_REF] Aulchenko | GenABEL: an R library for genome-wide association analysis[END_REF].

Conditional and joint association analysis

Next, to identify (additional) independent signal(s) from each significant GWAS locus, we performed stepwise model selection using GCTA-COJO analysis [START_REF] Yang | Conditional and joint multiple-SNP analysis of GWAS summary statistics identifies additional variants influencing complex traits[END_REF]. Here, we used a LD threshold of 0.80 to identify independent marker(s), and assumed complete independence for markers ±10 Mb away from a conditional marker. After each iteration the marker with lowest P-value was selected when it passed our GWAS threshold of <5 × 10 -8 and was included to conditional marker list to perform next iteration. When there were multiple markers with same (lowest) P-value, the first marker based on chromosome position was arbitrarily chosen. Finally, joint association analysis was perform on markers selected from the same chromosome.

Identification of candidate genes

We used bovine reference genome assembly UMD3.1 [START_REF] Zimin | A whole-genome assembly of the domestic cow, Bos taurus[END_REF][START_REF] Elsik | Bovine Genome Database: new tools for gleaning function from the Bos taurus genome[END_REF] for all the genomic coordinates presented in this study. To identify candidate genes within a QTL, we used Ensembl (release 94) webserver [START_REF] Zerbino | Ensembl 2018[END_REF] and manually inspected 1 Mb genomic region centering at the top marker (i.e. top marker ±500 kb surrounding region). From Ensembl, we retrieved annotations for genes and known phenotypes associated with those genes. However, when multiple genes were found to be located within the interval, for reporting, we preferred genes with fertility related phenotypes in other species (if available).

Results and Discussion

We performed single-marker association analysis for eight female fertility traits (Table 5.1) in Holsten, Jersey and Nordic Red Dairy cattle using imputed WGS SNPs, indels and deletions. Association test statistics were inflated (Table 5.2) when we only accounted for relatedness using genomic relationship matrix (GRM) derived from the 777 k markers in a LOCO approach. Therefore, following [START_REF] Jiang | A resource-efficient tool for mixed model association analysis of large-scale data[END_REF], we adjusted for both population stratification, using the first ten principal components (PCs) (calculated from a random set of 10% markers selected from the 29 autosomes) as fixed effects, and relationships among animals using the LOCO approach. Genomic inflation rates, lambda, before and after fitting the PCs were presented in Table 5.2. We observed a decrease in lambda values after the inclusion of PCs as cofactor, although it was still very far from the standard value of "1" under the assumption of the distribution of test statistics under the null-hypothesis. However, when a trait is under polygenic inheritance, i.e. when multiple (small effect) causal variants underlie the trait variation (also known as infinitesimal model), a substantial inflation of GWAS test statistics is expected even in the absence of population structure (Yang et al., 2011b). The extent of this inflation is governed by several factors, namely, heritability of the trait, study sample size, pattern of linkage-disequilibrium, and number of causal variants (for detail see Yang et al. (2011b)). The genomic inflation seen in this study despite the correction for population stratification and relatedness could be explained by the nature of polygenic inheritance of these fertility traits we analyzed, e.g. several non-zero effect variants in Bayesian mixture models as reported by Brøndum et al. (2015). The phenotypes used here are de-regressed breeding values for bulls with high reliability and are equivalent to a high heritability trait.

The long-range LD pattern seen in the study population [START_REF] De Roos | Linkage disequilibrium and persistence of phase in Holstein-Friesian, Jersey and Angus cattle[END_REF], van den Berg et al., 2016) could also explain inflated lambda values observed in this study. Among the three breeds, HOL had the highest sample size and therefore, highest power to identify QTL, was also the breed where we observed highest lambda values. However, the adjustment for both population stratification and relatedness could be very conservative; nonetheless, we expect that the test statistics will be robust against spurious associations and false discoveries.

The summary of single-marker associations, for eight female fertility traits in HOL, JER and RDC, is presented in Table 5.3 and corresponding Manhattan plots are presented in Figure 5.1, 5.2 and 5.3, respectively.

Deletions associated with female fertility traits

In HOL, we identified three deletions that had significant genome-wide association (P-value <5×10 -8 ) with ICF, FI and NRRc, respectively (Table 5.4). We identified three candidate genes, ACSS1, CFAP61 and HELLS, within these deletion loci (±500 kb), that are known to be associated with embryonic lethality, premature death, and male infertility in mouse (Table 5.5).

In RDC, 13 deletions reached genome-wide significance; associated phenotypes include seven fertility traits (Table 5.4). We did not find significant association between deletions and ICF. Among these 13 deletions, strongest association (p-value = 3.5×10 -42 ) was seen for a previously known deletion on chromosome 12 located between position 20,100,648 and 20,763,119 bp (esv4015629) that cause embryonic lethality in cattle (Kadri et al., 2014). The esv4015629 had significant association with all fertility traits analyzed, except ICF as expected. Kadri et al. (2014) also did not find association between esv4015629 and ICF. Several deletion loci showed association with multiple fertility traits. For example, FI, IFLc and AISc are the three highly correlated traits (correlation of 0.97 and 0. 91 between FI vs. IFLc, and FI vs. AISc, respectively (NAV, 2013)) that had seven associated deletion loci in common (Table 5.4). Several genes within these loci have known fertility related phenotypes, including embryonic lethality, in mouse (Table 5.5).

In JER, no deletion locus reached genome-wide significance. The top five associations included esv4014584 (Chr13:50733796-50735009; associated with ICF, P-value = 3.2×10 -6 ), esv4018559 (Chr6:33504998-33506227; associated with AISc, P-value = 3.4×10 -6 ), esv4019122 (Chr6:33738269-33738743; associated with AISc, P-value = 3.4×10 -6 ), esv4017388 (Chr1:145570019-145570346; associated with NRRc, P-value = 8.5×10 -6 ), and esv4017421 (Chr3:7140763-7141120; associated with NRRh, P-value = 1.8×10 -5 ). Among these five deletions, only esv4017388 located within a QTL associated with NRRc (Figure 5.1h). Candidate genes for this locus include ADARB1 (ENSBTAG00000017486; 52 kb upstream of esv4017388) and POFUT2 (ENSBTAG00000007818; 467 kb downstream of esv4017388). In mouse, ADARB1 causes postnatal lethality (MGI: 891999), and POFUT2 causes embryonic lethality during organogenesis as well as lethality between implantation and somite formation (MGI: 1916863).

Fertility index (FI) associated SNPs and indels

For SNPs and indels, we only presented FI associated independent GWAS loci. The seven FI component traits need further analyses to identify independent as well as pleotropic loci. COJO and multi-trait Meta analyses for these traits are in progress.

FI associated loci in HOL.

In HOL, 30,384 SNPs and 178 indels in total had reached GWAS significant Pvalue threshold for the eight fertility traits (Table 5.3 and Figure 5.1). After COJO analysis, we identified 37 FI associated independent GWAS loci out of these significant markers (Table 5.6). In Table 5.7, we presented the candidate genes for these loci. Top three FI associated loci are rs434006863 in chromosome 13 (P-value = 1.5×10 -22 ), rs380439408 in chromosome 24 (P-value = 5.4×10 -9 for this JER population. However, in this study, we did not find any significant association with FI for JER. Two factors, namely, sample size and control for genomic inflation, could explain this situation. For JER, we had ~4 times fewer samples compared with HOL and RDC, and QTL detection power was therefore lower; however, sample size between this two studies did not changed significantly (1,225 vs. 1,211 in this study). Difference between studies can be explained by the way to account for genomic inflation, much more stringent in this study: lambda of 1.88 in Hoglund et al. (2015b) vs. 1.27 in this study.

FI associated loci in RDC. Finally, in RDC, we identified 23,481 SNPs and 189 indels that had genomewide significant association with fertility traits (Table 5.3). After COJO analysis, we reported 23 FI associated independent GWAS loci, nine of which are located in chromosome 12 (Table 5.8). The LD pattern (correlation) of these nine lead markers is presented in Table 5.9. Functional annotations of the genes located within these QTL region (top markers ±500 kb) is presented in Table 5.10. Interestingly, the candidate genes within these loci also had known fertility related functions in cattle, mouse, and zebrafish (Table 5.10).

Conclusion

In this study, we identified several novel QTL for female fertility in Nordic dairy cattle. With our imputed data, we were able to identify a previously known deletion QTL in RDC; however, GWAS signals for this QTL were stronger with our imputed data than in the previous report. Other deletion loci also reached GWAS threshold but were not among the top signals of those loci. Those deletions could be in strong LD with the QTL. Interestingly, majority of the candidate genes within the reported QTL regions had established fertility related function in mouse and zebrafish. QTL identified in this study will be used in subsequent study to elucidate the effect of these markers in prediction accuracy for female fertility. Traits: number of inseminations per conception in cows (AISc) and in heifers (AISh); interval (number of days) from calving to first insemination (ICF); interval from first to last insemination in cows (IFLc) or in heifers (IFLh); fertility index (FI); non-return rate in cows (NRRc) or in heifers (NRRh). When a deletion was associated with multiple trait, then GWAS summary for trait with strongest signal (bold face) was presented in the table. No significant deletion was found in JER breed. 2.147 0.364 3.7×10 -9 Chr: Chromosome; Position: base pair positions in UMD3.1; A1/A2: Effect allele/other allele; Freq: frequency of A1 allele; b: allele substitution effect, se: standard error of b estimates; P: p-value from single marker GWAS. Subscript 'J' indicates values from conditional and joint analysis.

In RDC, we achieved highest prediction accuracy with a two-component model of 50k and FI associated markers (r = 0.578) followed by a three-component model comprising 50k, FI and deletion markers (r = 0.572). Prediction accuracy with 50k markers was the third best from the nine scenarios we considered in this study (Figure 6.1c). Our results are consistent with the previous report (Brøndum et al., 2015). Interestingly, we observed a relative gain of 10% in prediction reliability compared with that of Brøndum et al. (2015), although we overestimated bulls' genetic merit (Table 6.3).

Overall, predictability of imputed deletions alone was rather impressive considering the relative size of this marker-set (Table 6.1). For example, in RDC, even with ~4 times more markers in FI dataset, prediction accuracy was ~14% lower compared with that of using deletion markers only. Furthermore, we achieved ~80% of 50k markers' accuracy using 10 times fewer markers in the deletion dataset (0.451 with deletion vs. 0.566 with 50k; ~4k deletion markers vs. ~41k in 50k dataset).

Conclusion

In this study, we illustrated the potential of large chromosomal deletions as well as imputed WGS markers in predicting genomic breeding values in dairy cattle. For female fertility, we reported a substantial gain in prediction accuracy/reliability compared with previous reports in the same population. We showed that in a two-component model, including deletions along with 50k markers did not improve prediction accuracy compared with that of one-component 50k model. However, we also showed that they did not reduce accuracy, and this could mean that these imputed deletions did not add noise to our prediction model and could be useful for predicting fitness related traits. However, direct genotyping of these imputed deletions and WGS markers is necessary to explore the full potential of these markers as predictor. Here, 50k: Illumina BovineSNP50 BeadChip markers. DEL: large chromosomal deletions (size >50 bp) that are not in the GWAS marker list. GWAS: markers selected from single-marker GWAS of fertility index (FI) and five component traits, e.g. number of inseminations per conception in heifers (AISh) or cows (AISc), interval (number of days) from calving to first insemination in cows (ICF), interval (number of days) from first to last insemination in heifers (IFLh) or cows (IFLc). HOL: Holstein; JER: Jersey; RDC: Nordic Red Dairy cattle; GRM/G: genomic relationship matrix. . Here, we reported squared prediction accuracy, i.e. prediction reliability, which could be scaled using average reliability of FI in the test population (HOL: 0.645, JER: 0.587, and RDC: 0.654).

Recessive Lethals

Poor fertility is generally related to embryonic mortality. Embryonic mortality can be caused by any lethal mutation homozygous in an embryo and affecting any essential gene, not necessarily always related to reproduction. Early embryonic mortalities are not observed and therefore, usual mapping methods cannot be used because biological material from dead foetus is not available. Nevertheless, when the number of genotypes is high enough in a population, such defects can be indirectly observed by a deficit in some genotypes. This approach is much more powerful with haplotypes that are more informative than single SNP and better reflect identical by descent (IBD) chromosomal regions. This method proposed by VanRaden et al.

(2011) has been successfully applied in several breeds by different authors. In Normande breed, we have been able to map such a haplotype with a strong deficit in homozygous state and have identified the causal variant in CAD gene (Chapter 2). This study was rich in conclusions: 1) several previously described haplotypes in this breed were not confirmed, emphasizing the importance of a high quality imputation and phasing process and of a carefully adjusted statistical threshold accounting for the huge number of comparisons performed;

2) like other recessive lethals (Michot et al., 2017, Fritz et al., 2018), this mutation affected a very basic metabolic pathway (DNA base synthesis) and induced rapid mortality, when the compensatory metabolic stocks of maternal origin are depleted;

3) because of their strong effects, such mutations are always recent and therefore usually breed specific, but they can reach moderate frequencies (3-10%) because of the limited effective population size of our modern dairy breeds and strong genetic drift, as well as less efficient counter selection due to their recessive determinism; 4) because these defects are recessive, they are difficult to map with conventional GWAS, except when the frequency is already high and when the proportion of at risk mating becomes important. From this point of view, haplotype deficiency and GWAS are two complementary approaches to identify deleterious recessive mutations. 5) This haplotype deficiency study is a by-product of the large-scale genotyping for genomic selection. It requires very large numbers of genotypes and one can anticipate that the number of detected loci will increase over time, simply because detection power increases. It illustrates one out of many applications that can be derived from this large genotype data accumulating over time. Once discovered, such a mutation can be readily Chapter 7 121 added to SNP genotyping chip, disseminating this information to the whole population at nearly no additional cost. 6) Finally, from a practical point of view, because these defects are recessive and frequencies are usually below 10%, their economic impact is limited and does not justify a drastic elimination of the carriers. Instead, including this information in mating plans to avoid mating between carriers (or potentially carriers) is recommended. In addition, a careful monitoring of its frequency is important to avoid any further dissemination of the mutation and ensure a gradual elimination.

Use of Whole-Genome Sequence Variants

The markers used in majority of the genomic prediction studies were SNPs and (in fewer studies) indels. Furthermore, to our knowledge, no studies systematically explored the use of structural variants in the cattle genome for mapping or prediction purposes. This situation does not reflect the relative importance of SNP and SVs. SVs are less numerous than SNPs but the number of affected bases in the genome is of the same magnitude and each SV is much more likely to modify biological processes than any single SNP, especially when they affect one or even several genes. This situation of limited study is explained by the relative ease of access to the information for these two classes of variants. SNP data are readily available from SNP chips, without additional bioinformatics analyses. In contrast, with the noticeable but limited exception of SV breakpoints designed to be tested like SNPs by custom array (Boichard et al., 2018), most SVs are difficult to call.

Since our interest was to identify lethal mutations in dairy cattle, we choose to focus on large (>50 bp) chromosomal deletions besides SNPs and indels (Chapter 3-6). This is because deletion of an essential gene or its gene-regulatory region will certainly lead to loss-of-function (LOF) and can lead to lethality or a substantial reduction in fitness for animals homozygous for the deletion. Technological advances over the past decades for whole genome sequencing and availability of large WGS data in dairy cattle offer the possibility of testing different class of genomic variants for severe recessive genetic disorders. SVs can be called from sequence data with different algorithms (variation in read depth, discordant mapping of single reads or read pairs) but this remains tedious on a large scale and with a strong proportion of false positives. They can also be called from genotyping fluorescence signals and sometimes from Mendelian incompatibilities, but only for large SVs. However, SV genotype calling as well as imputation are not trivial and many challenges remain regarding the use of other type of genetics [START_REF] Meadows | Dissecting evolution and disease using comparative vertebrate genomics[END_REF]Lindblad-Toh, 2017, Van de Peer et al., 2017). While looking for candidate mutations or candidate genes, this evolutionary conservation across species could be a powerful tool. Genes that are essential for survival and fitness are usually under strong purifying selection compared with nonessential genes [START_REF] Hart | High-Resolution CRISPR Screens Reveal Fitness Genes and Genotype-Specific Cancer Liabilities[END_REF], Wang et al., 2015), have fewer paralogs and involve in more protein-protein interactions (Blomen et al., 2015[START_REF] Hart | High-Resolution CRISPR Screens Reveal Fitness Genes and Genotype-Specific Cancer Liabilities[END_REF]. These genes contain more singletons and very few deleterious variants at rare or low frequency (Blomen et al., 2015[START_REF] Hart | High-Resolution CRISPR Screens Reveal Fitness Genes and Genotype-Specific Cancer Liabilities[END_REF]. Including CAD gene, these features were common for majority of the previously detected lethal genes in dairy cattle (OMIA, 2018).

Conclusions

-Causal factors for recessive lethals could be detected using a combination of HHD screen and subsequent validation -Positional mapping and genotyping of large chromosomal deletions using WGS data provides opportunity for incorporating deletions in genomic studies -Copy-number status, such as copy-number loss or gain, could be inferred using SNP read depth data from the VCF file -Large deletions could be imputed with high accuracy along with SNPs and indels -Better imputation could provide stronger GWAS signal -Imputed WGS selected markers could improve within-breed prediction accuracy

Perspectives

-Detection of lethal haplotype, as well as lethal mutation will facilitate management of lethals in the population by avoiding at-risk mating. A first step is to include these markers on the SNP chip used for genomic selection and genotypes will be available to the whole genotyped population, making it possible to use it in selection and mating decisions.

-GWAS at the sequence level is supposed to include the causal variants. However, many association studies are limited to SNPs. Because SV are more likely to be causal than a SNP, it is important to include them in the studies, even if it increases genotype calling complexity.

-A balance between genetic gains and inbreeding should be maintained for sustainability in dairy cattle breeding. Rapid inbreeding trend favours both genetic drift and homozygotes proportions-two risk factors for fertility.

Chapter 7 127 -Breeding objectives should include both production and fitness traits to maintain a steady progress for all traits -Direct genotyping of imputed deletions and potential QTL is necessary to explore the full potential of these markers as predictor. It is possible to simplify SV genotyping process by SNP-like tests based on the SV breakpoints. Then large-scale information is readily available.

-Use of fewer predictive markers should be given preference for genomic prediction in small breeds, because there is less information to estimate each marker effect.

-Lastly, automated selection of predictive markers using deep learning should be explored

Figure 1 . 1 .

 11 Figure 1.1. An illustration of various classes of whole-genome sequence variants. (Adopted from presentation slides by Mesbah-Uddin, M., GenSAP Annual Meeting 2017, http://mbg.au.dk/fileadmin/gensap/events/2017/Session1/Md_Mesbah_Uddin_GenSAP2017.pdf)

  the read-depth variation to estimate SV genotypes, one could only study deletions (copynumber loss) or duplications (copy-number gain). An illustration of SNP read-depth based genotyping of large deletion is presented in Figure 1.4.

Figure 1 . 3 .

 13 Figure 1.3. Illustration of merging SV calls generated using multiple tools. Shaded-box represents confidence interval around the SV breakpoints and the straight horizontal line represents start and end of the breakpoint.

Figure 1 . 4 .

 14 Figure 1.4. Graphical presentation of SNP read-depth based genotyping of large deletion. CN: Copy-number; RD: read-depth; µj and varj represent mean and variance of j th Gaussian. (adopted from Mesbah-Uddin et al. (2018))

Figure 2

 2 Figure 2.1. Ranking of animals (n=2,240) for four udder traits (%). Black dots indicate the relative percentile of "Genetique" (the cow homozygous for NH7 haplotype) amongst the 2,239 paternal half-sisters. The percentile was calculated as: # 𝑜𝑜𝑓𝑓 𝑎𝑎𝑎𝑎𝑖𝑖𝑎𝑎𝑎𝑎𝑎𝑎𝑎𝑎 𝑤𝑤𝑖𝑖𝑤𝑤ℎ (≤) 𝑜𝑜𝑓𝑓 𝑎𝑎 𝑔𝑔𝑖𝑖𝑔𝑔𝑔𝑔𝑎𝑎 𝑎𝑎𝑠𝑠𝑜𝑜𝑖𝑖𝑔𝑔 𝑤𝑤𝑜𝑜𝑤𝑤𝑎𝑎𝑎𝑎 𝑎𝑎𝑎𝑎𝑖𝑖𝑎𝑎𝑎𝑎𝑎𝑎𝑎𝑎

Figure 2

 2 Figure 2.2. Fine-mapping at NH7 locus. (a) Pedigrees of the two inbred descendants of "Nonnic" (NORFRAM002977029292), "0520" (NORFRAF003708210520) and "Fanion" (NORFRAF001447822444), both carry a copy of the 26-marker long original NH7 haplotype and a recombinant NH7 haplotype facilitating fine-mapping at the locus. (b) Diagram showing original NH7 and recombinant haplotype located between the marker rs41610536 (Chr11:70,750,209) and rs109299742 (Chr11:72,476,622), along with 100 markers on both sides of the haplotype: upstream markers are from rs41569373 (Chr11:65,246,700) to rs41572406 (Chr11:70,703,633), and downstream markers are from rs109666167(Chr11:72,524,223) to rs110443631 (Chr11:78,468,627). The recombination occurred between marker rs29022293(Chr11:71,465,910) and rs29010798 (Chr11:71,520,448) in "Fanion", and between marker rs110176879 (Chr11:72,159,929) and rs109299742 (Chr11:72,476,622), in "0520". The resulting critical region is located between marker rs29022293 and rs109299742, i.e. between chromosomal interval of71,465,910 and 72,476,622 bp. Here, "NH7" denotes "Normande haplotype 7"; "?" denotes un-genotyped animal.

Figure 2 . 4 .

 24 Figure 2.4. Protein domains and conserved sequences in CAD. (a) Schematic representation of Bovine CAD protein domains from NCBI's Conserved Domain Database (CDD) (domain architecture ID 11486011; accessed from webserver "https://www.ncbi.nlm.nih.gov/cdd"). The diagram represents four functional protein domains, such as: (i) a small subunit of the glutamine-dependent form of carbamoyl-phosphate synthase II (GATase1_CPSase from 178-354; accession: cd01744), (ii) a large subunit of carbamoyl-phosphate synthase II (CPSase-II from 392-1440; accession: TIGR01369)-which catalyzes the first step in de novo pyrimidine biosynthesis, (iii) a subunit of dihydroorotase (DHOase from 1460-1806; accession: cd01316)-which catalyzes the third step in de novo pyrimidine biosynthesis, and (iv) a subunit of aspartate/ornithine carbamoyltransferase (OTCase_N from 1925-2065; accession: pfam02729)-which catalyzes the second step in de novo pyrimidine biosynthesis. (b) Multi-species sequence alignment representing a 61-residue long region in CPSase-II domain of CAD protein (centering at BTA11 p.Y452C) among nine vertebrates, five metazoans, two plants, three fungi, three protists, and two bacterial species. An arrow highlights the position of the CAD mutation, in bovine and orthologus mutation in humans.An asterisk (*) at the bottom of the alignment represents residues conserved in all 24-species (shaded vertical boxes); a colon (:) represents conservation among residues with strong similarity; a period (.) represents conservation among residues with weak similarity. BTA: Bos taurus autosome; HSA: Homo sapiens autosome.

Figure 2 . 5 .

 25 Figure 2.5. Genotype of NH7 homozygote at Chr11 g.72399397T>C mutation. The CAD mutation is indicated by the shaded vertical box.

Figure 2 . 6 .

 26 Figure 2.6. Effect of CAD mutation on four-udder traits. * Here, P-values are from Student's t-test on 250-animals with AGgenotype (carriers) vs 2,912 animals with AA-genotypes (non-carriers).

Figure 3 . 1 .

 31 Figure 3.1. Number of ascertained deletions relative to variant allele count. Here, variant allele frequency (VAF) is expressed in terms of variant allele count. Deletions down to an allele count of 1 (VAF=0.0026 and 0.0032, in cattle and humans, respectively) are also represented here. Human deletion calls by Mills et al. (Mills et al., 2011) were downloaded from 1K Genomes Project FTP server (ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/pilot_data/paper_data_sets/companion_papers/mapping_structural_variation).

Figure 3

 3 Figure 3.2. A ~525-KB deletion on chromosome 23 discovered using Genome STRiP (a), and resolved breakpoint sequences from TIGRA and BLAT search (b). Shaded bases are a 5-bp microhomology at breakpoint junction. (This figure was drawn and modified using Inkscape version 0.91.)

Figure 3 . 3 .

 33 Figure 3.3. Experimental validation of the ~525 KB deletion on chromosome 23. (a) PCR amplification across (left) and within the deletion (right) for two carrier (∆/+) and two homozygous wild-type (+/+) animals. Water: negative control, M: molecular weight marker (GeneRuler™ 100bp DNA ladder, Fermentas). (b) Sequence trace of the 359 bp amplicon bridging the breakpoint. Shaded bases are a 5-bp microhomology at breakpoint junction. (This figure was drawn and modified using Inkscape version 0.91.)

Figure 3 . 4 .

 34 Figure 3.4. Population diversity. (a) Heterozygous deletions per genome. (b) Homozygous deletions per genome. (c) Singletons per genome. Only high-confidence genotype calls are included. The y-axis in (a), (b), and (c) represents the number of heterozygous, homozygous, and singleton deletions per genome, respectively. HOL: Holstein; JER: Jersey; RDC: Nordic Red Cattle.

Figure 4 . 1 .

 41 Figure 4.1. Imputation accuracy for 107 deletions on chromosome 29 using the reference panels RefBeagle vs. RefBeagle_SHAPEIT. (a) Leave-out-trials. For each deletion, 77 trials were performed: 76 trials of 10 animals each and 1 trial of 12 animals. Imputation accuracy was calculated using the Pearson's correlation, r 2 , between the true and imputed alternative allele dosages. (b) Downstream performance. Using the two reference panels, a target population of ~13k animals with the 777k genotypes was imputed for the deletions (along with other WGS markers) of chromosome 29 using Minimac. Each dot represent one deletion; dashed lines indicate r 2 value of 0.80. RefBeagle: reference panel phased using Beagle; RefBeagleSHAPEIT:

Figure 4 . 2 .

 42 Figure 4.2. Downstream imputation performance between RefBeagle vs. RefBeagle_SHAPEIT on test chromosome 29. (a) SNPs. In total, 326,838 SNPs were imputed from this chromosome; each dot represents a SNP. (b) INDELs. There were 2,723 imputed indels in this chromosome; each dot represents an indel. The x-and y-axis represents corresponding Minimac r 2 value; dashed lines indicate Minimac r 2 value of 0.80. RefBeagle: reference panel phased using Beagle; RefBeagleSHAPEIT: reference panel phased using Beagle followed by re-phasing using SHAPEIT.

Figure

  Figure 4.3. Imputation accuracy of WGS SNPs, indels and deletions on the 29 bovine autosomes. Each dot represents an

Figure 4 . 4 .

 44 Figure 4.4. Principal component analysis using imputed WGS SNPs, indels and deletion. Here, the analyses were performed on (a) ~137k SNPs, (b) ~12k indels and ~4.5k deletions using PLINK v2. These variants were randomly selected from the imputed WGS variants with MAF >1%. PC: principal component; HOL: Holstein; JER: Jersey; RDC: Nordic Red Dairy cattle.

Figure 5 . 1 .

 51 Figure 5.1. Manhattan plots of single-marker GWAS in Holstein cattle. Each dot represent one WGS marker; dashed horizontal line is GWAS significance threshold of 5×10 -8 .

Figure 5 . 2 .

 52 Figure 5.2. Manhattan plots of single-marker GWAS in Jersey cattle. Each dot represent one WGS marker;dashed horizontal line is GWAS significance threshold of 5×10-8 .

Figure 5 . 3 .

 53 Figure 5.3. Manhattan plots of single-marker GWAS in Nordic Red Dairy cattle. Each dot represent one WGS marker; dashed horizontal line is GWAS significance threshold of 5×10 -8 .

(

  Figure1.3), or requires population scale sequencing data to boost the discovery signal. Small SV discovery population size, which is typical for WGS dataset in most research labs, is a further constraint on subsequent imputation and genetic study. One solution to small sample size problem could be the use of shared genomic resources available through different

collaboration, such as the 1000 Bull Genomes Project where variant call format (VCF) files are shared among the collaborators. It is important to note, however, that BAM files are not always available for these resources and we need methods adapted to VCF files. These files include read-depth data supporting each SNP genotype call and this information can be used to estimate SV genotypes, as shown in Chapter 4 of this thesis. Of course, the success of this approach will heavily depend on the sequencing read-depth, size of SV region, and more importantly availability of SNPs within the interrogated region. Since, this technique analyses

A missense mutation (p.Tyr452Cys) in the CAD gene compromises reproductive success in French Normande cattle Md Mesbah-Uddin, 1

	2.1 Abstract
	We scanned the genome of 77,815 Normande cattle with different Illumina SNP-chips to map recessive
	We addressed embryonic lethal mutations using homozygous haplotype deficiency (HHD). We detected two novel
	this broad aim under five separate studies that I have presented in Chapter 2 to Chapter 6 in haplotypes on chromosomes 11 and 24, but did not confirm six previously reported haplotypes. The one on
	Chapter 2. chromosome 11 showed a marked reduction in conception rates and moderate decrease in non-return rate in this thesis with the objectives of each study listed below:-at-risk versus control mating, supporting late embryonic mortality. After fine-mapping and analyzing whole-
	1. To validate (or invalidate) six previously reported lethal haplotypes in Normande genome sequences we prioritized a missense mutation in CAD (g.72399397T>C; p.Tyr452Cys)-a gene
	cattle, to identify new loci, and when possible, to identify causal variants in this dairy encoding a protein essential for de novo pyrimidine biosynthesis-as a candidate causal variant. This
	breed (Chapter 2) transition mutation replaces a tyrosine residue, which is perfectly conserved among living organisms, by a
	cysteine residue in the carbamoyl-phosphate synthetase 2 domain of the protein. A single animal was 2. To map large chromosomal deletions from whole-genome sequences of 67 Holstein, confirmed to be homozygous for the mutation based on Sanger sequencing. However, large-scale genotyping
	27 Jersey, and 81 Nordic Red Dairy cattle, and to analyse their (identified deletions) of the candidate variant with Illumina EuroG10k Beadchip revealed an absence of live homozygote in a
	population-genetic properties (Chapter 3) panel of 33,323 Normande animals and an absence of carriers in 348,593 animals from 19 other cattle breeds.
	3. To estimate deletion genotype-likelihood using SNP read-depth data from VCF file, to ,2* Chris Hoze, 2,3 Pauline Michot, 2,3 Anne Barbat, 2 Rachel Lefebvre, 2 Mekki
	build a haplotype reference panel that include SNPs, indels, and deletions, and to Boussaha, 2 Goutam Sahana, 1 Sébastien Fritz, 2,3 Didier Boichard, 2 and Aurélien Capitan 2,3*
	jointly impute these variants into the existing SNP array-typed animals (Chapter 4)
	4. To perform genome-wide association studies with imputed SNPs, indels and deletions
	for eight female fertility traits in Holstein, Jersey and Nordic Red Dairy cattle
	(Chapter 5)
	5. To investigate the effect of imputed WGS selected markers on accuracy of genomic
	prediction for female fertility in Holstein, Jersey and Nordic Red Dairy cattle (Chapter
	6)

Table 2 .1. List of Normande haplotypes with deficit of homozygotes (haplotypes NH1-NH6 from Fritz et al. (2013); NH7-NH8 from the present study)

 2

	Haplotype	Chromosome	Start	End	Frequency		No. of Homozygotes
	ID		Position 1	Position	(%)	Expected Observed	P-value †
	NH1 *	24	38,086,180	39,153,166	1.9 (1.8) 2	9 (12)	6 (0)	0.21 (5.3×10 -4 )
	NH2	1	145,682,206 146,833,973	6.3 (3.8)	394 (49)	314 (14)	1.7×10 -5
								(5.7×10 -7 )
	NH3 *	4	92,261,682	93,828,702	4.6 (5.9)	189 (41)	159 (10)	0.01 (1.3×10 -6 )
	NH4	6	37,704,254	38,869,785	4.0 (5.2)	172 (38)	108 (12)	1.08×10 -7
								(2.5×10 -5 )
	NH5 *	7	3,661,458	4,592,227	18.0 (1.9)	2,764	2,735	0.29 (6×10 -7 )
						(58)	(20)
	NH6 *	15	59,789,154	61,044,897	12.7 (1.9)	1,525	1,503	0.29 (3×10 -5 )
						(45)	(17)
	NH7	11	70,750,209	72,476,622	3.7	57	1	1.0×10 -23
	NH8	24	52,606,336	54,542,635	1.9	40	7

Table 2 .2. Estimated effect of NH7 and NH8 on conception rate (CR) and non-return rate at 56 days (NRR56) in at-risk mating 1

 2 At risk mating are defined by mating between a carrier AI bull and a daughter of a carrier sire. Non-risk mating include all other combinations (carrier x non-carrier, non-carrier x carrier, and non-carrier x non-carrier).

	Haplotype	Trait	Category	Number of mating Non-risk At-risk	Total	Mean (µ %)	Effect on fertility (%) Expected 2 Observed
	NH7	CR NRR56	Heifer Cow Heifer Cow	2,981,552 7,008,668 1,870,062 3,999,038	14,474 40,618 7,513 22,494	2,996,026 7,049,286 1,877,575 4,021,532	57.37 47.20 76.62 67.13	-7.30 -6.01 -9.75 -8.55	-5.53 -5.56 -1.76 -1.52
	NH8	CR NRR56	Heifer Cow Heifer Cow	3,014,200 7,095,799 1,893,593 4,047,881	2,497 3,962 1,520 2,261	3,016,697 7,099,761 1,892,073 4,050,142	57.44 47.21 76.60 67.16	-7.25 -5.95 -9.66 -8.47	0.64 0.38 0.47 0.55
	1

Table 2 .3. Integrative annotation of whole-genome sequence (WGS) variants strongly correlated with the NH7 haplotype

 2

	CHROM:POS 1	REF/ALT	Sequence Ontology	Pearson's correlation (r)	Gene-level scores dN/dS 2 RVIS (%) 3	Z-score 4	Position-level scores 5 SIFT GERP++	PhyloP	PhastCons
	Chr11:71144269	C/T	Intron variant of PLB1 (ENSBTAG00000018669)	1	0.318	2.53 (98.5%)	-2.33	-6	1.14	0.18	0
	Chr11:71151125	C/T	Intron variant of PLB1 (ENSBTAG00000018669)	1	0.318	2.53 (98.5%)	-2.33	-	-3.37	-2.32	0
	Chr11:71282689	C/T	Intergenic variant	0.82	-	-	-	-	-	-	-
	Chr11:71695719	T/A	Intron variant of BRE (ENSBTAG00000031335)	1	0.095	-0.38 (32.33%)	1.13	-	-0.72	0.90	0
	Chr11:71766806	A/C	Intron variant of BRE (ENSBTAG00000031335)	1	0.095	-0.38 (32.33%)	1.13	-	3.45	1.32	0.016
	Chr11:71772162	G/A	Intron variant of BRE (ENSBTAG00000031335)	0.82	0.095	-0.38 (32.33%)	1.13	-	-4.32	-0.48	0
	Chr11:72384541	C/T	Intron variant of CAD (ENSBTAG00000017894)	0.82	0.051	-4.84 (0.18%)	4.85	-	1.69	0.70	0
	Chr11:72399397	T/C	Missense variant of CAD (ENSBTAG00000017894)	1	0.051	-4.84 (0.18%)	4.85	0.01	5.71	6.85	1

Table 2 .4. Genotypes at CAD g.72399397A>G mutation in 20 French cattle breeds 1 1

 2 In total, 381,916 animals were genotyped for the candidate-SNP with EuroG10K custom SNP chip(Boichard et al., 2018).

2 Genotype: AA= "CAD p.452Tyr/Tyr"; AG= "CAD p. 452Tyr/Cys"; GG= "CAD p.452Cys/Cys".

Table 2 .5. Concordance between NH7 status and genotypes for CAD g.72399397A>G in 33,323 Normande

 2

	cattle
	Genotype	+/+ 1	NH7/+	NH7/NH7	Total
	g.72399397 A/A	31,368	8	0	31,376
	g.72399397 A/G	237	1,710	0	1,947
	g.72399397 G/G	0	0	0	0
	Total	31,605	1,718	0	33,323
	1 "+": Non-NH7 haplotype

  • C for 3 min, 2) 35 cycles of 30 sec denaturation (94 • C), 30 sec hybridization (65.2 • C), 30 sec elongation (72 • C) and a final 3 min elongation (72 • C). PCR products were separated on a 2% agarose gel, purified and directly sequenced using the BigDye Terminator Cycle Sequencing Kit (Applied Biosystems, CA, US). Electrophoresis of sequencing reactions was performed on

	3500xL Genetic Analyzers (Applied Biosystems, CA, US), and sequences were visualized with Sequencher
	5.4.6 (Gene Codes Corporation, MI, US). A 977 bp control amplification, with a primer pair within the
	deletion (forward primer: 5′-CCCAATGCAAAATCACAAAA -3′ and reverse primer: 5′-
	CCAGAAAAGCTACACTTGAACTGA -3′), was performed using the same reaction conditions as above
	except hybridization was performed at 59.8 • C.

  We calculated fold enrichment for a trait, such as for Health related QTL: (No. of Health QTL on Deletions ∕ Total QTL on deletions) ∕ (Total Health QTL ∕ Total QTL in the dataset), and statistical significance using Fisher's exact test (two sided).

	Chapter 3
	considered for QTL enrichment analysis.
	Dickinson et al. (Dickinson et al., 2016)) using Wilcoxson test. Reported causal
	genes for cattle were also retrieved from OMIA database (http://omia.angis.org.au/, last accessed on May
	10, 2016) for dN/dS comparison.
	50

We retrieved cattle QTL from QTLdb database

[START_REF] Hu | Developmental progress and current status of the Animal QTLdb[END_REF]

(release 31; accessed on January 6, 2017); autosomal QTL from Holstein, Jersey, Nordic Red Cattle and Ayrshire, associated to any of the six trait classes, e.g. "Reproduction", "Milk", "Production", "Exterior", "Meat and Carcass", and "Health", were

Table 3 .1. False discovery rate (FDR) estimates of Genome STRiP' deletion calls using 777K BovineHD BeadChip intensity data

 3

	Array-probe Overlap	A ¥	B §	FDR*
	One-probe	497	28	11.3%
	>1 array-probe	206	3	2.9%
	>2 array-probe	113	1	1.8%
	Overall	703	31	8.8%
	¥ A=No. of sites with p-value; § B=No. of sites with p-value ≥ 0.5; *FDR estimates were based on Wilcoxson rank sum test
	using BovineHD chip intensity data of 26 Holstein animals. FDR was calculated as (B/A ×2 ×100).

Table 3 .2. Enrichment of QTL on deletions

 3

	Trait Classes *	Fold Enrichment	P value (Fisher's test)
	Health	2	8.91×10 -10
	Reproduction	1.5	7.4×10 -11
	Milk	0.8	2.45×10 -7
	Exterior	0.5	1.85×10 -4
	Production	0.5	0.002
	Meat and Carcass	0.5	0.058

* Trait classes are from cattleQTLdb

  1 Mb for the WGS variants, and 2 Mb for the 50k and 777k variants, following the

	SHAPEIT	guidelines	(last	accessed	on	May	8,	2019,	from:

Table 4 .2. Accuracy of imputation from the 777k genotypes to the WGS SNPs, indels and deletions of the 29 bovine autosomes

 4

	Category	Minimac r 2 (SD) SNPs	INDELs	Deletions
	Overall	0.928 (0.207)	0.900 (0.217)	0.904 (0.220)
	Common variants (MAF >5%)	0.954 (0.163)	0.927 (0.176)	0.957 (0.148)
	Low frequency variants (1≤ MAF ≤5%)	0.889 (0.250)	0.874 (0.245)	0.927 (0.154)
	Rare variants (MAF <1%)	0.527 (0.392)	0.549 (0.370)	0.713 (0.343)
	0.1%	0.119 (0.183)	0.173 (0.222)	0.313 (0.365)
	0.2%	0.201 (0.231)	0.262 (0.265)	0.633 (0.364)
	0.3%	0.275 (0.269)	0.349 (0.271)	0.785 (0.280)
	0.4%	0.330 (0.306)	0.388 (0.294)	0.764 (0.294)
	0.5%	0.443 (0.348)	0.498 (0.328)	0.822 (0.243)
	0.6%	0.503 (0.371)	0.552 (0.351)	0.819 (0.217)
	0.7%	0.601 (0.370)	0.617 (0.353)	0.846 (0.232)
	0.8%	0.653 (0.377)	0.671 (0.356)	0.803 (0.275)
	0.9%	0.718 (0.355)	0.691 (0.350)	0.880 (0.180)
	1.0%	0.731 (0.357)	0.722 (0.341)	0.871 (0.207)
	Total number of imputed WGS SNPs = 14,070,960, Indels = 122,054, and deletions = 5,730

Table 5 .1. Summary statistics of female fertility traits (de-regressed breeding values) analyzed

 5 AISc: number of inseminations per conception in cows, AISh: number of inseminations per conception in heifers, ICF: Interval (number of days) from calving to first insemination (in cows), IFLc: interval (number of days) from first to last insemination in cows, IFLh: interval (number of days) from first to last insemination in heifers, FI: fertility index, NRRc: non-return rate in cows, NRRh: non-return rate in heifers.

	Trait 1	Number	Holstein Mean (SD)	Range	Number	Jersey Mean (SD)	Range	Number	Nordic Red Dairy Cattle Mean (SD)	Range
	AISc	5,342	100.30 (11.38)	65-142	1,179	99.54 (10.11)	69-126	4,291	102.57 (10.56)	73-136
	AISh	5,487	100.06 (10.00)	72-131	1,204	100.22 (12.70)	60-133	4,345	95.27 (12.62)	59-132
	ICF	5,342	98.38 (11.19)	65-131	1,186	103.96 (10.9)	75-135	4,299	99.21 (10.20)	69-130
	IFLc	5,361	98.1 (12.09)	60-146	1,192	99.74 (11.53)	65-136	4,309	101.63 (10.17)	73-134
	IFLh	5,478	98.75 (10.97)	68-133	1,199	101.66 (11.37)	68-134	4,341	95.18 (11.64)	62-129
	FI	5,577	98.86 (10.77)	68-142	1,211	101.14 (11.42)	69-131	4,471	100.44 (9.67)	74-129
	NRRc	5,353	100.77 (10.21)	72-135	1,189	97.10 (9.67)	67-126	4,309	102.34 (10.56)	73-136
	NRRh	5,487	101.65 (11.38)	70-137	1,207	99.95 (10.86)	65-130	4,351	96.76 (10.98)	65-129
	1

Table 5 .2. Genomic inflation rate (lambda) before and after correcting for population stratification using first ten principal components (PCs)

 5 AISc: number of inseminations per conception in cows, AISh: number of inseminations per conception in heifers, ICF: Interval (number of days) from calving to first insemination (in cows), IFLc: interval (number of days) from first to last insemination in cows, IFLh: interval (number of days) from first to last insemination in heifers, FI: fertility index, NRRc: non-return rate in cows, NRRh: non-return rate in heifers.

	Trait 1	Without PCs	Holstein	With PCs	Without PCs	Jersey	With PCs	Nordic Red Dairy Cattle Without PCs With PCs
	AISh	1.76		1.65	1.39		1.29	1.41	1.37
	AISc	2.04		1.83	1.31		1.24	1.56	1.50
	ICF	2.38		2.19	1.36		1.29	1.63	1.57
	IFLh	1.93		1.81	1.37		1.29	1.42	1.37
	IFLc	2.24		2.02	1.32		1.25	1.59	1.53
	FI	2.26		2.05	1.33		1.27	1.59	1.52
	NRRh	1.99		1.84	1.37		1.29	1.41	1.38
	NRRc	2.20		1.98	1.39		1.28	1.53	1.48
	1

Table 5 .3. Summary of single-marker association results (markers with P-value <5 × 10 -8 )

 5 AISc: number of inseminations per conception in cows, AISh: number of inseminations per conception in heifers, ICF: Interval (number of days) from calving to first insemination (in cows), IFLc: interval (number of days) from first to last insemination in cows, IFLh: interval (number of days) from first to last insemination in heifers, FI: fertility index, NRRc: non-return rate in cows, NRRh: non-return rate in heifers.

	Trait	Holstein	Jersey	Nordic Red Dairy Cattle
	AISh	583	0	10,004
	AISc	2,012	0	16,089
	ICF	21,539	8	280
	IFLh	932	0	7,302
	IFLc	7,645	1	16,289
	FI	11,091	0	14,916
	NRRh	1,057	0	8,036
	NRRc	2,008	8	9,182
	1

Table 5 .4. Fertility associated large deletions (P-value <5×10 -8 )

 5 Marker ID: deletion IDs are from Database of Genomic Variants archive (DGVa), were extracted from Ensembl release 94; Position base pair positions in UMD3.1; A1/A2: Effect allele/other allele; Freq: frequency of A1 allele; b: allele substitution effect, se: standard error of b estimates; P: p-value from single marker GWAS.

	Chr Marker ID	Position	Size of Deletion	Trait	A1/A2	Freq	b	se	p
	Holstein cattle
	13	esv4017813	13:42,945,993-42,946,953 961 bp	ICF	DEL/G 0.271 -1.050 0.186 1.7×10 -8
	13	esv4018227	13:40,037,936-40,039,382 2.45 kb	FI	DEL/C 0.482 -0.894 0.155 8.1×10 -9
	26	esv4015626	26:16,220,854-16,223,815 2.96 kb	NRRc	DEL/T	0.097 1.418	0.248 1.0×10 -8
	Nordic Red Dairy cattle
	6	esv4013561	6:31,705,882-31,706,925	1.04 kb	NRRh	DEL/C 0.287 -1.341 0.242 3.0×10 -8
	6	esv4018710	6:37,999,532-37,999,979	448 bp	NRRh, NRRc	DEL/C 0.352 -1.525 0.236 9.4×10 -11
	6	esv4019275	6:39,473,114-39,473,611	498 bp	NRRc	DEL/A 0.389 1.260	0.215 4.9×10 -9
	6	esv4011295	6:39,711,454-39,712,517	1.06 kb	AISc	DEL/A 0.423 -1.208 0.215 1.9×10 -8
	6	esv4014663	6:39,856,457-39,856,658	202 bp	NRRc, AISc	DEL/A 0.466 1.502	0.211 1.2×10 -12
	12	esv4019299	12:17163364-17163728	365 bp	FI, AISh, AISc, IFLh, IFLc, NRRh	DEL/G 0.250 -1.655 0.226 2.2×10 -13
	12	esv4015629	12:20,100,648-20,763,119 662 kb	FI, AISh, AISc, IFLh, IFLc, NRRh, NRRc	DEL/G 0.137 -3.985 0.293 3.5×10 -42
	12	esv4012561	12:20,570,478-20,571,118 641 bp	FI, IFLc, IFLh, AISc	DEL/A 0.092 -2.010 0.323 5.0×10 -10
	12	esv4017726	12:22,961,548-22,961,822 275 bp	IFLc, AISc, FI	DEL/A 0.435 1.394	0.207 1.8×10 -11
	12	esv4013857	12:35,299,066-35,299,339 274 bp	IFLc	DEL/T	0.127 -1.789 0.315 1.4×10 -8
	12	esv4015687	12:43,999,233-44,009,032 9.8 kb	IFLc, AISc, NRRc, FI	DEL/T	0.102 -2.255 0.331 9.8×10 -12
	25	esv4014328	25:12,757,355-12,758,442 1.09 kb	FI, AISh, IFLh, AISc, IFLc	DEL/G 0.059 -2.432 0.401 1.3×10 -9
	25	esv4014654	25:13,705,829-13,706,486 658 bp	AISh, AISc, IFLc, IFLh, FI, NRRh	DEL/A 0.050 -3.833 0.522 2.1×10 -13
	Chr: Chromosome;

Table 5 .5. Annotation of fertility associated deletion loci

 5 Annotations are based on UMD3.1 and Ensembl release 94. Information in the parenthesis represents the distance between the top marker and the start position of the corresponding gene; MGI: Mouse Genome Informatics (http://www.informatics.jax.org/); RGD: Rat Genome Database (https://rgd.mcw.edu/wg/); MP: Mouse Phenotype (https://www.mousephenotype.org/); OMIM: Online Mendelian Inheritance in Man (https://www.omim.org/)

	Marker ID	Gene

Table 5 .6. Fertility index associated loci in Holstein cattle identified from GWAS with imputed sequence variants

 5 intron variant; 2 intergenic variant. Annotations are based on UMD3.1 and Ensembl release 94. Information in the parenthesis represents the distance between the top marker and the start position of the corresponding gene. MGI: Mouse Genome Informatics (http://www.informatics.jax.org/); RGD: Rat Genome Database (https://rgd.mcw.edu/wg/); MP: Mouse Phenotype (https://www.mousephenotype.org/); OMIM: Online Mendelian Inheritance in Man (https://www.omim.org/); OMIA: Online Mendelian Inheritance in Animals; ZFIN: The Zebrafish Information Network (http://zfin.org/).

	Chr	Top Marker	Position	A1/A2	Freq	Single-marker analysis b se	p	Conditional and joint analysis bJ bJ_se pJ
		rs137209107	68,642,776	A/G	0.143	-1.377	0.228	1.5×10 -9	-1.377	0.229	1.7×10 -9
		rs384339926	113,461,025	G/A	0.129	1.285	0.224	1.0×10 -8	1.285	0.225	1.1×10 -8
		rs207792133	23,885,602	C/T	0.131	-1.284	0.226	1.4×10 -8	-1.284	0.227	1.5×10 -8
		rs43306949	50,233,158	G/A	0.159	-1.147	0.209	3.9×10 -8	-1.147	0.209	4.2×10 -8
		rs137588365	29,789,003	A/G	0.376	-0.997	0.159	3.6×10 -10	-0.997	0.160	4.1×10 -10
		rs42253052	68,281,044	A/G	0.313	1.071	0.170	3.3×10 -10	1.071	0.171	3.8×10 -10
		rs109192619	79,679,521	T/C	0.479	-0.902	0.159	1.4×10 -8	-0.902	0.159	1.6×10 -8
		rs135727662	47,510,332	A/G	0.489	0.954	0.153	4.8×10 -10	0.954	0.154	5.4×10 -10
		rs385277008	65,156,839	G/A	0.172	1.270	0.207	7.9×10 -10	1.270	0.207	9.0×10 -10
		rs210968504	86,998,676	G/T	0.182	-1.127	0.205	3.7×10 -8	-1.127	0.205	4.0×10 -8
		rs110866176	93,466,329	C/G	0.449	-1.096	0.157	3.0×10 -12	-1.096	0.158	3.7×10 -12
		rs473589185	93,099,110	T/G	0.060	-1.827	0.327	2.3×10 -8	-1.827	0.328	2.5×10 -8
		rs382943070	27,286,359	T/C	0.112	-1.491	0.245	1.2×10 -9	-1.491	0.246	1.4×10 -9
		rs208577191	49,271,809	G/A	0.323	0.953	0.170	2.2×10 -8	0.953	0.171	2.4×10 -8
		rs380495923	68,425,703	T/A	0.320	-1.186	0.168	1.6×10 -12	-1.186	0.169	2.0×10 -12
		rs134897493	81,159,267	T/C	0.354	-1.123	0.168	2.4×10 -11	-1.123	0.169	2.9×10 -11
		rs134324803	23,631,737	G/C	0.242	-1.087	0.181	2.0×10 -9	-1.087	0.182	2.3×10 -9
		rs209793386	42,987,167	T/G	0.168	1.307	0.209	3.9×10 -10	1.307	0.210	4.4×10 -10
		rs384444008	25,178,930	G/T	0.214	1.040	0.190	4.6×10 -8	1.301	0.196	3.0×10 -11
		rs385034802	29,587,693	T/C	0.445	0.345	0.159	2.9×10 -2	0.990	0.168	4.2×10 -9
		rs434006863	36,616,682	T/C	0.148	1.849	0.215	8.8×10 -18	2.192	0.224	1.5×10 -22
		rs136213414	51,160,223	A/G	0.157	-1.437	0.215	2.6×10 -11	-1.437	0.216	3.1×10 -11
		rs42685926	34,562,085	C/G	0.132	1.361	0.232	4.4×10 -9	1.361	0.233	4.9×10 -9
		rs109077219	8,194,073	T/G	0.395	1.054	0.160	4.4×10 -11	1.054	0.161	5.2×10 -11
		rs41853464	72,010,870	A/G	0.371	-0.999	0.162	7.0×10 -10	-0.999	0.162	7.9×10 -10
		rs382112789	31,104,391	T/G	0.175	-1.166	0.208	2.0×10 -8	-1.166	0.208	2.2×10 -8
		rs450967176	44,576,200	T/G	0.117	-1.717	0.244	1.9×10 -12	-1.717	0.245	2.3×10 -12
		rs382522431	61,080,660	A/T	0.135	1.380	0.232	2.8×10 -9	1.380	0.233	3.1×10 -9
		rs383235276	4,888,887	T/C	0.253	1.053	0.186	1.4×10 -8	1.053	0.186	1.6×10 -8
		rs41911242	26,618,825	C/T	0.316	-0.978	0.166	3.8×10 -9	-0.978	0.166	4.2×10 -9
		rs715807267	36,483,718	A/G	0.100	1.630	0.259	2.9×10 -10	1.630	0.259	3.3×10 -10
		rs380439408	29,556,826	A/C	0.233	-1.413	0.187	4.0×10 -14	-1.413	0.188	5.4×10 -14
		rs134963365	46,786,442	T/G	0.485	0.988	0.160	6.2×10 -10	0.988	0.160	7.1×10 -10
		rs379672347	62,207,595	G/T	0.248	-1.010	0.181	2.6×10 -8	-1.010	0.182	2.9×10 -8

1

Table 5 .8. Fertility index associated loci in Nordic Red Dairy cattle identified from GWAS with imputed sequence variants

 5

	Chr.	Top Marker	Position	A1/A2	Freq.	Single-marker analysis b se p	Conditional and joint analysis bJ bJ_se pJ
		rs135868307	131,241,313	G/T	0.413	-1.198	0.202	2.9×10 -9	-1.198	0.203	3.3×10 -9
		rs450253115	93,470,575	G/C	0.047	-2.559	0.450	1.3×10 -8	-2.787	0.453	7.6×10 -10
		rs42490886	97,240,584	T/C	0.138	-1.498	0.279	8.1×10 -8	-1.646	0.281	4.7×10 -9
		rs210258782	36,974,279	T/C	0.194	-1.618	0.238	1.0×10 -11	-1.618	0.239	1.3×10 -11
		rs110759619	39,633,323	G/A	0.481	1.373	0.189	3.3×10 -13	1.373	0.190	4.4×10 -13
		rs379594438	95,384,366	G/A	0.055	2.381	0.421	1.6×10 -8	2.381	0.423	1.8×10 -8
		rs470078356	87,251,892	T/C	0.013	-4.800	0.861	2.4×10 -8	-4.800	0.864	2.7×10 -8
		rs132823348	103,372,929	A/G	0.190	-1.343	0.242	2.7×10 -8	-1.343	0.243	3.0×10 -8
		rs386105326	11,016,415	G/A	0.113	-2.336	0.308	3.2×10 -14	-2.336	0.310	4.6×10 -14
		rs385475503 (rs797826108)	21,149,563	G/T	0.141	-4.087	0.290	4.9×10 -45	-11.123	0.445	8.0×10 -138
		rs384604017	29,738,986	T/C	0.285	-0.482	0.221	2.9×10 -2	9.110	0.464	7.0×10 -86
		rs133494863	29,790,646	C/T	0.274	0.734	0.224	1.0×10 -3	-5.561	0.359	3.6×10 -54
		rs432066150	30,453,718	G/A	0.291	-0.962	0.226	2.1×10 -5	3.602	0.316	4.7×10 -30
		Chr12:31178612	31,178,612	C/T	0.128	-2.820	0.298	3.1×10 -21	-14.778	0.636	2.5×10 -119
		Chr12:43320610	43,320,610	A/G	0.087	-2.424	0.338	7.1×10 -13	-2.424	0.340	9.4×10 -13
		rs377982349	55,272,444	C/T	0.086	-1.944	0.348	2.3×10 -8	-1.944	0.349	2.6×10 -8
		rs133220791	69,485,528	T/C	0.234	-1.374	0.230	2.2×10 -9	-1.374	0.231	2.5×10 -9
		rs207566536	46,143,358	G/A	0.136	-1.794	0.277	1.0×10 -10	-1.794	0.279	1.2×10 -10
		rs378749244	56,756,424	T/C	0.475	1.246	0.206	1.4×10 -9	1.246	0.207	1.6×10 -9
		rs136268337	59,753,524	T/C	0.353	1.200	0.210	1.1×10 -8	1.200	0.211	1.2×10 -8
		rs135743550	39,835,929	A/G	0.302	-1.275	0.216	3.5×10 -9	-1.275	0.217	4.1×10 -9
		rs382027875	13,426,672	T/C	0.055	-3.091	0.420	1.8×10 -13	-3.091	0.422	2.4×10 -13
		rs481408336	35,317,193	A/G	0.061	2.147	0.363	3.3×10 -9

Table 6 .1. Summary of selected markers in the training dataset, used in constructing genomic-relationship matrices for genomic prediction

 6

	Type of variants 1		HOL	JER	RDC	GRM
	50k		40,917	36,400	41,394	G50k
	DEL		3,683	3,101	4,225	GDEL
		FI	12,990	3,465	16,541	GFI
		AISc	2,186	3,521	3,402	GAISc
	GWAS	AISh IFLc	673 929	3,679 1,507	1,901 1,408	GAISh GIFLc
		IFLh	461	2,710	128	GIFLh
		ICF	16,694	3,517	398	GICF
	1

Table 6 .2. Variance explained by different classes of variants

 6 Here, the order of the variance components are same as in corresponding values. V(G): Genetic variance at marker set G; V(G)/Vp: proportion of phenotypic variance explained by marker-set G; G corresponds to G50k, GFI and GDEL.

	Variance		HOL		JER		RDC
	Components 1	V(G)	V(G)/Vp	V(G)	V(G)/Vp	V(G)	V(G)/Vp
	G50k	65.721	72.77%	82.491	67.38%	69.016	71.49%
	GDEL	50.615	55.14%	60.951	50.99%	43.708	47.34%
	GFI	313.260	84.69%	329.756	83.92%	129.174	65.68%
	G50k + GDEL	63.096 + 2.752	69.90% + 3.05%	88.532 + 0.000	72.56% + 0.00%	66.336 + 3.061	68.64% + 3.17%
	G50k + GFI	45.233 + 38.300	41.49% + 35.13%	38.747 + 249.261	12.14% + 78.10%	54.704 + 25.971	50.41% + 23.93%
	G50k + GFI + GDEL	43.749 + 38.057 + 1.613	40.21% + 34.98% + 1.48%	37.555 + 249.350 + 1.209	11.76% + 78.11% + 0.38%	54.345 + 25.957 + 0.413	50.08% + 23.92% + 0.38%
	1

Table 6 .3. Slope of the regression of phenotype on the genomic predicted breeding value for fertility index using the 50k, selected GWAS loci, large deletions (>50 bp) and combinations of these markers Breed Regression slope of true vs. predicted breeding values

 6

		50k DEL	FI	50k & DEL	50k & FI	B50k, FI & DEL	FI & 5-subtraits	B50k, FI & 5-subtraits	B50k, FI, 5-subtraits & DEL
	HOL	0.94	0.88	0.81	0.94	0.78	0.77	0.36	0.67	0.67
	JER	1.18	1.22	0.99	1.12	0.18	0.18	0.42	0.33	0.32
	RDC	0.91	0.90	0.75	0.91	0.88	0.87	0.49	0.72	0.73

Table 6 .4. Reliability of genomic prediction (without scaling 1 ) for fertility index

 6

	Breed	Reliability of prediction 50k DEL FI	50k & DEL	50k & FI	B50k, FI & DEL	FI & 5-subtraits	B50k, FI & 5-subtraits	B50k, FI, 5-subtraits & DEL
	HOL	0.408 0.286 0.296	0.408		0.379	0.371	0.188	0.346	0.344
	JER	0.354 0.265 0.441	0.353		0.163	0.161	0.244	0.205	0.198
	RDC	0.320 0.203 0.151	0.318		0.334	0.327	0.153	0.293	0.291
	1 Reliability of prediction could be scaled using:	𝑎𝑎𝑠𝑠𝑠𝑠𝑎𝑎𝑖𝑖𝑔𝑔𝑠𝑠 𝑠𝑠𝑜𝑜𝑖𝑖𝑖𝑖𝑔𝑔𝑎𝑎𝑎𝑎𝑤𝑤𝑖𝑖𝑜𝑜𝑎𝑎 𝑏𝑏𝑔𝑔𝑤𝑤𝑤𝑤𝑔𝑔𝑔𝑔𝑎𝑎 𝐸𝐸𝐸𝐸𝐸𝐸𝑎𝑎 𝑎𝑎𝑎𝑎𝑠𝑠 𝐺𝐺𝐸𝐸𝐸𝐸𝐸𝐸𝑎𝑎 𝑖𝑖𝑔𝑔𝑎𝑎𝑖𝑖𝑎𝑎𝑏𝑏𝑖𝑖𝑎𝑎𝑖𝑖𝑤𝑤𝑟𝑟 𝑜𝑜𝑓𝑓 𝐸𝐸𝐸𝐸𝐸𝐸𝑎𝑎 𝑖𝑖𝑎𝑎 𝑤𝑤𝑔𝑔𝑎𝑎𝑤𝑤𝑖𝑖𝑎𝑎𝑔𝑔 𝑠𝑠𝑎𝑎𝑤𝑤𝑎𝑎𝑎𝑎𝑔𝑔𝑤𝑤

Acknowledgements: This study is part of the BOVANO project (ANR-14-CE19-0011) funded by the French Agence Nationale de la Recherche (Paris, France) and APIS-GENE (Paris, France). Md Mesbah-Uddin (MMU) benefited from a joint grant from the European Commission within the framework of the Erasmus-Mundus joint doctorate "EGS-ABG". MMU's PhD project is also supported by the Center for Genomic Selection in Animals and Plants (GenSAP) funded by Innovation Fund Denmark (grant 0603-00519B). The authors are grateful to the partners of the 1000 bull genomes consortium for their excellent collaboration. We are grateful to the Genotoul bioinformatics platform Toulouse Midi-Pyrenees (Bioinfo Genotoul) for providing computing and storage resources.

Acknowledgements

Md Mesbah-Uddin benefited from a joint grant from the European Commission within the framework of the Erasmus-Mundus joint doctorate "EGS-ABG". This research was supported by the Center for Genomic Selection in Animals and Plants (GenSAP) funded by Innovation Fund Denmark (grant 0603-00519B).

Acknowledgement

This research was supported by the Center for Genomic Selection in Animals and Plants (GenSAP) funded by Innovation Fund Denmark (grant 0603-00519B). Md Mesbah-Uddin benefited from a joint grant from the European Commission within the framework of the Erasmus-Mundus joint doctorate 'EGS-ABG'.

Authors also acknowledge the 1,000 Bull Genomes Project for sharing the VCF files. We are also grateful to the Genotoul bioinformatics platform Toulouse Midi-Pyrenees (Bioinfo Genotoul) for providing computing and storage resources.

Acknowledgement

We are grateful to the Nordic Cattle Genetic Evaluation (NAV), Aarhus, Denmark for providing the phenotypic data used in this study and Viking Genetics, Randers, Denmark for proving blood/semen samples for genotyping. This research was supported by the Center for Genomic Selection in Animals and Plants (GenSAP) funded by Innovation Fund Denmark (grant 0603-00519B). Md Mesbah-Uddin benefited from a joint grant from the European Commission within the framework of the Erasmus-Mundus joint doctorate 'EGS-ABG'. Authors also acknowledge the 1000 Bull Genomes Project for sharing the VCF files. Authors are also grateful to the Genotoul bioinformatics platform Toulouse Midi-Pyrenees (Bioinfo Genotoul) for providing computing and storage resources.

Acknowledgement

We are grateful to the Nordic Cattle Genetic Evaluation (NAV), Aarhus, Denmark for providing the phenotypic data used in this study and Viking Genetics, Randers, Denmark for proving blood/semen samples for genotyping. This research was supported by the Center for Genomic Selection in Animals and Plants (GenSAP) funded by Innovation Fund Denmark (grant 0603-00519B). Md Mesbah-Uddin benefited from a joint grant from the European Commission within the framework of the Erasmus-Mundus joint doctorate 'EGS-ABG'. Authors also acknowledge the 1000 Bull Genomes Project for sharing the VCF files. Authors are also grateful to the Genotoul bioinformatics platform Toulouse Midi-Pyrenees (Bioinfo Genotoul) for providing computing and storage resources.

selektion. det overordnede formål gennem fem artikler / manuskripter, der praesenteres i kapitel 2 til kapitel 6 i denne afhandling. for at prioritere indsamling af store maengder genotypeinformation om kandidatmutationen for at validere foreløbige resultater. deletioner ved brug af nordiske malkekvaegsracer. Denne undersøgelse . Dernaest beskrives funktionelle betydninger af de identificerede er opståede. Denne undersøgelse danner basis for de efterfølgende genetiske analyser, som er praesenteret i denne ph.d.-afhandling. deletioner. for prøverne; Dette var dog kendt for egenskaber med staerke imputerede) store deletioner når fertilitetsegenskaber kortlaegges. kromosomale deletioner. Ved baseret på denne ph.d.-afhandling.

vii

Data availability

Availability of data and materials

All relevant data are included in the manuscript and its additional files. WGS SNPs and INDEL are available at http://www.1000bullgenomes.com/doco/1000bulls_v6_annotated_snps.tab.gz, http://www.1000bullgenomes.com/doco/1000bulls_v6_annotated_indels.tab.gz, and large deletions are at ftp://ftp.ebi.ac.uk/pub/databases/dgva/estd234_Mesbah-Uddin_et_al_2017/. The imputation pipelines are available at https://github.com/MMesbahU/ImputeDelPipeline.git.

Abbreviations

AISc number of inseminations per conception in cows

AISh number of inseminations per conception in heifers Supplementary Material

Abstract

Large genomic deletions are potential candidate for loss-of-function, which could be lethal as homozygote.

Analyzing whole genome data of 175 cattle, we report 8,480 large deletions (199bp to 773KB) with an overall false discovery rate of 8.8%; 82% of which are novel compared to deletions in the dbVar database.

Breakpoint sequence analyses revealed that majority (24 of 29 tested) of the deletions contain microhomology/homology at breakpoint, and therefore, most likely generated by microhomology-mediated end joining. We observed higher differentiation among breeds for deletions in some genic-regions, such as ABCA12, TTC1, VWA3B, TSHR, DST/BPAG1, and CD1D. The genes overlapping deletions are on average evolutionarily less conserved compared to known mouse lethal genes (p-value=2.3×10 -6 ). We report 167 natural gene knockouts in cattle that are apparently nonessential as live homozygote individuals are observed. These genes are functionally enriched for immunoglobulin domains, olfactory receptors, and MHC classes (FDR=2.06×10 -22 , 2.06×10 -22 , 7.01×10 -6 , respectively). We also demonstrate that deletions are enriched for health and fertility related QTL (2 and 1.5 fold enrichment, Fisher's p-value=8.91×10 -10 and 7.4×10 -11 , respectively). Finally, we identified and confirmed the breakpoint of a ~525 KB deletion on Chr23:12,291,817,087 (overlapping BTBD9,GLO1 and DNAH8), causing stillbirth in Nordic Red Cattle.

Key words: dairy cattle, structural variants, whole genome sequence, population genetics, loss-of-function

Introduction

Embryonic lethality has become a challenge to cattle breeders, especially for dairy cattle where a limited number of bulls were extensively used in breeding for fast genetic progress in economic traits like milk and protein yield (Charlier et al., 2016). An estimated yearly loss of ~$10.74 million is attributed to known recessive lethals in four dairy cattle breeds from USA only, where Holstein accounts for ~70% of the total losses, followed by Jersey, Brown Swiss, and Ayrshire (Cole et al., 2016). Hence, understanding the genomic architecture of cattle populations is important, now more than ever, for optimizing genetic gain while constraining negative impact of deleterious mutations responsible for genetic defects and inbreeding depression.

Unlike single nucleotide polymorphism (SNP) and small insertion or deletion (indel), structural variants (SVs), i.e. DNA alterations larger than 50 base pairs (bp) that include insertions, deletions, duplications, inversions, and translocations [START_REF] Weischenfeldt | Phenotypic impact of genomic structural variation: insights from and for human disease[END_REF], are the least explored polymorphisms in cattle.

SVs contribute substantially to phenotypic variations and have a wide-spectrum of impact ranging from beneficial to lethal in both humans [START_REF] Weischenfeldt | Phenotypic impact of genomic structural variation: insights from and for human disease[END_REF][START_REF] Zarrei | A copy number variation map of the human genome[END_REF] and animals (Bickhart and Liu, 2014). The phenotypic impact of SVs in cattle is well evident from numerous studies. For example, Xu et al. [START_REF] Xu | Genome wide CNV analysis reveals additional variants associated with milk production traits in Holsteins[END_REF] showed that a combination of SNPs with SVs could explain additional genetic variance underlying milk production traits, while Charlier et al. (Charlier et al., 2012), [START_REF] Schutz | The Holstein Friesian Lethal Haplotype 5 (HH5) Results from a Complete Deletion of TBF1M and Cholesterol Deficiency (CDH) from an ERV-(LTR) Insertion into the Coding Region of APOB[END_REF], and Kadri et al. (Kadri et al., 2014) showed the lethal effect of large deletions in dairy cattle.

Sequence alignment to reference genome and SNPs/indels calling

Raw sequencing reads were filtered and FASTQ files were aligned to bovine reference genome assembly UMD3.1 using BWA software (Li and Durbin, 2009) to produce BAM files for subsequent variant calling.

In 1KBGP, SNPs and indels were called using SAMtools 0.1.18 mpileup software (Li et al., 2009), while GATK v1.6 software (McKenna et al., 2010) was used for Nordic WGS data (detailed method in (Daetwyler et al., 2014) and (Brondum et al., 2014), respectively). For all the analysis, bovine genome assembly UMD3.1 was used as the reference genome.

Discovery and genotyping of deletions

Structural variants can be detected from NGS data based on sequence signatures such as (discordant) readpair (RP), split-read (SP) and read-depth (RD), as well as de novo assembly of reads (Alkan et al., 2011).

However, approaches based on only one sequence signature could be constrained by high false discovery rate (FDR) (Handsaker et al., 2011), hence we employed a population scale SV detection method called "Genome STRucture in Populations (Genome STRiP)" (Handsaker et al., 2011) -which leverages technical (e.g. RP and RD signals) and population-level sequence features (e.g. coherence around shared alleles, and heterogeneity of evidentiary sequences in different genomes) for accurate discovery of deletions, and determines genotype (allelic state) of each locus from read-depth using a Gaussian mixture model.

Genome STRiP.

For deletion discovery and genotyping Genome STRiP software version 2.00.1678 (Handsaker et al., 2011) was used. Following the documentation, we built a custom reference metadata bundle for cattle samples that includes alignability mask, copy-number mask (CN2 mask), ploidy map, gender map. Alignability mask represents sites on the reference genome that are uniquely alignable by sequence read of a certain length (readLength). Our WGS data was a mixture of different Illumina pairedend reads ranging from 90 to 101 bp (Q1 = 90, median = 100, and Q3 = 100), hence genome alignability mask was prepared with readLength value of 90 using 'ComputeGenomeMask' utility from Genome STRiP.

Copy-number mask (CN2 mask), i.e. regions on the reference genome unlikely to be copy-number variable in most individuals, was produced for the bovine assembly UMD3.1 excluding sex chromosome X, unplaced contigs, and repeat sequences (retrieved from RepeatMasker track of UCSC Table Browser (Karolchik et al., 2004), accessed on July 4, 2016).

Preprocessing, deletion discovery and genotyping. We ran the preprocessing Queue script (dry run) to emit all the commands, prepared bash scripts to run in Portable Batch System job scheduler, and executed these commands proving 175 BAM files (one for each sample) as input.

Large deletions (100 bp ≤ size ≤ 1 MB) were discovered and filtered using SVDiscovery Queue script.

Discovered sites were filtered (default filters) if (i) the site contained too high or too low read pileup, (ii) read-pairs spacing was inconsistent with a single segregating deletion, (iii) read-depth and read-pair evidences were inconsistent across samples, (iv) read depth differences were not significant, and (v) read

Validation of deletions

Validation using 777K BovineHD BeadChip intensity data. We validated deletion calls using 777K

BovineHD BeadChip (Illumina, San Diego, CA, USA) intensity data on 26 Holstein samples that were both whole-genome sequenced and 777K chip typed. We calculated false discovery rate (FDR) for the deletion call-set using IntensityRankSum (IRS) test implemented in Genome STRiP. Intensity file was prepared from raw chip intensity data following the guideline for IRS test. Overall FDR for the call-set was calculated as two times the fraction of sites with IRS p-value ≥0.5 (i.e. sites with IRS p-value ≥0.5 to the sites with valid p-value). Details of IRS test could be found in (Handsaker et al., 2011[START_REF] Sudmant | An integrated map of structural variation in 2,504 human genomes[END_REF].

Validation by targeted assembly of breakpoint. Targeted iterative graph routing assembler (TIGRA-0.4.3) software [START_REF] Chen | TIGRA: a targeted iterative graph routing assembler for breakpoint assembly[END_REF] was used, with default parameters, for assembling deletion breakpoint sequences from a set of randomly selected deletions along with three previously known deletions segregating in the study populations. TIGRA extracted all reads mapped to 500 bp upstream and 50 bp downstream of start coordinate, and 50 bp upstream and 500 bp downstream of end coordinate of a given deletion; and reads were then assembled iteratively using de Bruijn graph assembler with multiple k-mers (e.g. 15 bp followed by 25 bp). We aligned the assembled contigs to UMD3.1 using Cow BLAT Search [START_REF] Kent | BLAT--the BLAST-like alignment tool[END_REF] from UCSC Genome Browser (https://genome.ucsc.edu/cgi-bin/hgBlat) to visualize and infer breakpoints from the alignments.

Validation by PCR and amplicon sequencing. We validated a previously reported ~525 KB deletion segregating in Nordic Red Cattle (Sahana et al., 2016) using PCR and amplicon sequencing. Genomic DNA was extracted as described previously by [START_REF] Miller | A simple salting out procedure for extracting DNA from human nucleated cells[END_REF] from semen sample of two bulls carrying the deletion and two non-carriers. The PCR reaction was done with the DyNAzyme II DNA Polymerase (Thermo Fisher, MA, US) in a 30µl volume of 1x PCR buffer, 0.2mM dNTPs, 10pmol primer mix (forward primer: 5′-AAGCCACCACAATGAGAAGC -3′ and reverse primer: 5′-TTTGGGGTAGGAGAAGTAGGG -3′) and 50 ng of genomic DNA. The cycling conditions were the breeds. Similar population structure (and substructure within Nordic Red Cattle) has been reported using genome-wide SNPs [START_REF] Mao | Genome-wide association studies of growth traits in three dairy cattle breeds using whole-genome sequence data[END_REF]. Nordic Red Cattle from Denmark showed closer relationship with the Holstein in our WGS samples (Figure S5). This is consistent with the known history of Holstein interrogation in Danish Red cattle, as previously reported based on imputed WGS SNP analysis on a larger sample [START_REF] Mao | Genome-wide association studies of growth traits in three dairy cattle breeds using whole-genome sequence data[END_REF]. It is also known from admixture analysis that genomes of Nordic Red Cattle are a mosaic of multiple ancestral populations, i.e. more ancestral components in Nordic Red Cattle than in Holstein and Jersey [START_REF] Bovine Hapmap | Genomewide survey of SNP variation uncovers the genetic structure of cattle breeds[END_REF][START_REF] Brondum | Reliabilities of genomic prediction using combined reference data of the Nordic Red dairy cattle populations[END_REF]; our deletion-based PCA largely corroborate that. VST analysis. We analyzed population stratification in terms of VST [START_REF] Redon | Global variation in copy number in the human genome[END_REF], a measure highly correlated with Wright's fixation index (FST) [START_REF] Wright | Evolution in Mendelian Populations[END_REF], to identify population differentiation. We calculated VST for each deletion pairwise amongst the breeds (e.g. Holstein vs Jersey+Nordic Red Cattle) from VAFs. We identified 158 highly stratified deletions (pairwise VST ≥ mean + 4 standard deviations)

among the breeds (Figure 3.6). Around 27% of these deletions overlap genic elements, i.e. exons, introns, or (upstream/downstream) untranslated regions (UTRs), and remaining 73% are intergenic variants (Tables S3-5). There were eleven sites shared between Holstein and Nordic Red Cattle, two sites between Holstein and Jersey, and one site between Nordic Red Cattle and Jersey.

Among these sites were gene variants, such as ABCA12 (Chr2:103,682,772-103,684,297 in Holstein & Jersey) associated with growth and development [START_REF] Cole | A genome-wide association study of calf birth weight in Holstein cattle using single nucleotide polymorphisms and phenotypes predicted from auxiliary traits[END_REF][START_REF] Xu | Genomic signatures reveal new evidences for selection of important traits in domestic cattle[END_REF], TTC1 (Chr7:73,725,513-Chapter 3 56 73,725,918 in Holstein & Nordic Red Cattle) with cold tolerance (Howard et al., 2014), VWA3B (Chr11:3,521,329-3,522,551 in Holstein & Nordic Red Cattle) with milk glycosylated kappa-casein percentage (Buitenhuis et al., 2016), and were intergenic variants, such as Chr15:41,393,393-41,393,780 (in Jersey & Nordic Red Cattle) and Chr20:26,812,159-26,812,834 (in Holstein & Jersey) overlap QTL associated with calving traits [START_REF] Mcclure | A genome scan for quantitative trait loci influencing carcass, post-natal growth and reproductive traits in commercial Angus cattle[END_REF], Sahana et al., 2011), Chr20:45,816,245-45,820,519 (in Holstein & Nordic Red Cattle) with meat and carcass trait (McClure et al., 2012), and Chr23:49,778,653-49,782,567 (in Holstein & Nordic Red Cattle) with body weight [START_REF] Snelling | Genome-wide association study of growth in crossbred beef cattle[END_REF]. We also identified a highly differentiated fertility associated gene DST/BPAG1 [START_REF] Lobago | Clinical features and hormonal profiles of cloprostenol-induced early abortions in heifers monitored by ultrasonography[END_REF], Cole et al., 2011) (Chr23:3,486,232-3,486,603) [START_REF] Araujo | Use of a candidate gene array to delineate gene expression patterns in cattle selected for resistance or susceptibility to intestinal nematodes[END_REF][START_REF] Sandri | Transcriptome profiles of whole blood in Italian Holstein and Italian Simmental lactating cows diverging for genetic merit for milk protein[END_REF], and also reported differentially expressed post intra-mammary infection [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]. Majority of these stratified deletion regions are novel compared to previous CNV studies in cattle, and therefore are interesting targets to investigate large deletions undergoing genetic drift or artificial selection.

Functional impact of deletions

We annotated all the deletions using Variant Effect Predictor (Ensembl 87). Around 71% (6,019 SVs) variants were intergenic and remaining 29% (2,461 SVs) overlapped genic elements, such as exons, introns, and untranslated regions (UTR). On average, high frequency gene disrupting deletions were somewhat depleted compared to intergenic variants (VAFintergenic > VAFgenic, p-value=0.04; one-sided Wilcoxon test).

Furthermore, we observed many common genic deletions. These genes are relatively less conserved, and majority has multiple paralogs (discussed later). However, deletions on known essential genes were only observed as heterozygote with relatively low VAF (<3%), and generally were private to a specific breed. For example, FANCI deletions (cause brachyspina (Charlier et al., 2012)) were only observed in Holstein, and RNASEH2B deletions (cause embryonic lethality (Kadri et al., 2014)) in Nordic Red Cattle.

Selective constraints on genes overlapping deletions. The relative abundance of high frequency genic and intergenic variants indicate that majority of these intersected genes are non-essential, and thus did not affect the viability or fecundity of the carriers. To test this hypothesis, we analyzed the selective constraints between deleted genes (overlap of any genic element) and known mouse lethal genes (from Dickinson et al. [START_REF] Dickinson | High-throughput discovery of novel developmental phenotypes[END_REF]) in terms of dN/dS ratio of cow-mouse 1-to-1 orthologs (Figure 3.7). Here, high dN/dS values indicate low selective constraints on genes, and low value indicates high constraints. We found that genes in deletions have significantly higher dN/dS ratios than lethal genes (p-value 2.3×10 -6 ; one-sided Wilcoxon test), and thus are evolutionarily less conserved. This is consistent with the rate of evolution seen in essential and non-essential genes -where mutations in essential genes were under strong purifying selection and thus evolved slowly (low dN/dS ratio), while non-essential genes were under relaxed selection, and hence, evolved faster (high dN/dS ratio) [START_REF] Hurst | Do essential genes evolve slowly?[END_REF]. Nonetheless, robustness of these processes is also evident in the evolution of human essential genes. Interestingly, ~77% human essential genes could even be traced back to pre-metazoans (Blomen et al., 2015).

Nonessential genes in cattle. In total, we found 5,000 deletions for which at least one individual was homozygous. In the set, we analyzed homozygous deletions in genes to find natural gene knockouts. We found 167 deleted genes (transcript-ablation or complete deletion) corresponding to 115 independent deletions that are apparently nonessential based on the occurrences of live homozygote individuals. This is ~45% more than the previous report (Boussaha et al., 2015). Nonetheless, we found ~44% fewer genes compared to in humans (240 nonessential genes) [START_REF] Sudmant | An integrated map of structural variation in 2,504 human genomes[END_REF], which could be due to the differences in sample size (175 vs 2,504 individuals) and study populations (3 vs 26 populations in human).

Among these genes, ~83% (139 genes) are protein-coding, 12% pseudogenes, and the rest are different types of small RNAs (Table S6). Most of these genes belong to multigenic families and are not highly conserved (median cow-mouse dN/dS of 0.17 vs OMIA genes dN/dS of 0.11; Figure S6), as expected for homozygous deletion [START_REF] Sudmant | An integrated map of structural variation in 2,504 human genomes[END_REF]. Moreover, this set of genes are functionally enriched in immunoglobulin domains, olfactory receptors, and MHC classes (FDR = 2.06×10 -22 , 2.06×10 -22 , 7.01×10 -6 , respectively), along with other related domains (Tables S7-9). Similar functional enrichment of nonessential genes was also seen in humans [START_REF] Sudmant | An integrated map of structural variation in 2,504 human genomes[END_REF]

Abstract

Genotype imputation, often focused on single-nucleotide polymorphisms (SNPs) and small insertions and deletions (indels; size <100 bp), is a crucial step for association mapping and estimation of genomic breeding values. Here, we present strategies to impute genotypes for large chromosomal deletions (size >100 bp), along with SNPs and indels in cattle. The pipelines include a strategy for extending the whole-genome sequence reference panel for large deletions, a two-step genotype refinement approach using Beagle4 and SHAPEIT2 software, and finally, joint imputation of SNPs, indels and large deletions to the existing SNP array-typed population using Minimac3 software. Using these pipelines we achieved an imputation accuracy of r 2 >0.6 at minor allele frequencies as low as 0.7% for SNPs and indels, and 0.2% for large deletions. This highlights the potential of our approach to build haplotype reference panel and impute different classes of sequence variants across a wide allele frequency-spectrum with high accuracy.

Key words: Gaussian mixture model, imputation, structural variant, copy-number variation, chromosomal deletion

Introduction

The population-scale whole-genome resequencing in different cattle breeds has generated a comprehensive catalog of DNA polymorphisms, including single-nucleotide polymorphisms (SNPs), small insertions and deletions (indels) (Daetwyler et al., 2014), copy-number variations (CNVs) (Letaief et al., 2017), and structural variants (SVs) (Boussaha et al., 2015, Chen et al., 2017). In livestock, there is a growing interest in understanding the phenotypic impacts of SVs, e.g., large deletions and duplications, because a substantial proportion of SVs overlap genes and may have a biological effect (Bickhart and Liu, 2014). Such genomic studies require SV genotypes in large cohorts of animals, which are rarely available. However, one way to generate such SV datasets is to design a customized SNP array using the precise SV breakpoints and genotype the population directly. The custom part of the EuroG10K BeadChip (Boichard et al., 2018) is an example of this. Alternatively, one can perform imputation for common bi-allelic SVs in an approach similar to that of bi-allelic SNPs and indels-since common SVs are often tagged well by the surrounding SNP haplotypes (as shown by studies in humans, e.g., Conrad et al. (2010) and Handsaker et al. (2015)). The later approach will allow us to impute SVs into the large cohorts of animals genotyped using the existing SNP arrays.

However, the majority of whole-genome sequence (WGS) based studies in cattle were either constrained by (small) sample size [START_REF] Shin | Deleted copy number variation of Hanwoo and Holstein using next generation sequencing at the population level[END_REF][START_REF] Gao | CNV discovery for milk composition traits in dairy cattle using whole genome resequencing[END_REF][START_REF] Mielczarek | Analysis of copy number variations in Holstein-Friesian cow genomes based on whole-genome sequence data[END_REF] or (unavailability of) SV genotypes (Chen et al., 2017). Therefore, building a suitable reference panel for imputing SVs into large cohorts of SNP array-typed cattle population remains a challenge. Earlier, we reported ~8,500 large deletions 68 Format (VCF) files of the 1,000 Bull Genomes Project (1KBGP). In this study, we present a strategy for computing deletion genotype likelihoods using the SNP read depth information from a VCF file when access to BAM files are not available, building a haplotype reference panel and jointly imputing WGS SNPs, indels, and deletions.

Methods

SNP Genotype Dataset

A total of 13,307 JER, HOL and RDC were genotyped using either Illumina BovineSNP50 (50k, hereafter)

or BovineHD BeadChip (777k, hereafter). The 50k and the 777k datasets comprised 9,932 (838 JER, 5,371 HOL, and 3,723 RDC) and 3,375 (835 JER, 1,215 HOL, and 1,325 RDC) animals, respectively.

From both datasets, we excluded SNPs that were monomorphic, had a GenCall <0.60, or deviated from Hardy-Weinberg proportions (P value ≤1 × 10 -5 ) in any of the three breeds. We also removed SNPs, which overlapped a deletion locus from our deletion dataset (see "WGS Dataset" subsection "Deletions"). After filtering, the 50k and 777k datasets had 47,237 and 697,296 SNPs, respectively.

WGS Dataset

In total, 772 animals were available in our WGS dataset (92 JER, 132 RDC and 548 HOL). Among these animals, 175 and 597 were from the Nordic sequence project [START_REF] Brøndum | Strategies for imputation to whole genome sequence using a single or multi-breed reference population in cattle[END_REF] and the Run-6 of the 1KBGP (Daetwyler et al., 2014), respectively. GATK v1.6 (McKenna et al., 2010) and SAMtools 0.1.18 mpileup (Li et al., 2009) software were used to call SNPs and indels from the Nordic and the 1KBGP WGS dataset, respectively (for detail, see [START_REF] Brøndum | Strategies for imputation to whole genome sequence using a single or multi-breed reference population in cattle[END_REF] and Daetwyler et al. (2014)).

SNPs and indels.

In this study, we only considered Phred-scaled genotype likelihoods (PL) for bi-allelic SNPs and indels, which occurred in both the Nordic and the 1KBGP dataset. We removed a marker from the dataset when any of the following conditions was met, namely: (i) minor allele count (MAC) less than 5

(vcftools' --mac 5), (ii) deviation from Hardy-Weinberg proportion in any one breed (P-value ≤1 × 10 -5 ),

(iii) missing genotypes in more than 10% of the animals (--max-missing 0.1), (iv) marker within 5 bp from one another (--thin 5), (v) marker quality less than 50 (--minQ 50), (vi) markers within 100 bp of gaps in the UMD3.1 bovine genome assembly [START_REF] Zimin | A whole-genome assembly of the domestic cow, Bos taurus[END_REF], or (vii) markers within 1 Kb of a deletion locus from our deletion list (discussed later). After filtering, 14,800,299 bi-allelic markers remained in our WGS dataset:

14,678,220 SNPs and 122,079 indels.

Deletions.

In addition, Mesbah-Uddin et al. (2018a) estimated PL for ~8,500 deletions on Nordic samples from the sequences aligned to UMD3.1 (BAM files) using GenomeSTRiP-2 software (Handsaker et al., 2011). Among them, we selected 5,798 deletions after a thinning of 100 Kb using VCFtools' "--thin 100000" option, i.e., excluding deletions within 100 Kb distance from one other. From this dataset, we extracted deletion PL for the Nordic animals, while for the remaining 597 animals, we estimated PL using a modified version of the read depth (RD) genotyping method used by Mesbah-Uddin et al. (2018b), as explained below.

Imputation

The imputation was performed in two steps using Minimac3 v.2.0.1 software (Das et al., 2016) (hereafter referred to as Minimac) with default parameters. In the first step, we imputed 9,932 animals from 50k to 777k using the 777k haplotype reference panel consisted of 3,375 animals. There were few animals common in both 50 k and 777 k dataset or 777 k and WGS dataset. Here, we kept the sample with highest genotype density (i.e. WGS >777k >50k) in the reference panel and remove the sample from the corresponding imputation target dataset. We then combined both imputed and genotyped 777k dataset to prepare a nonredundant target of 12,966 animals (1,640 JER,6,371 HOL and 4,955 RDC), after excluding the overlap between the 777k and WGS dataset. Next, using the haplotype reference panel of 772 WGS animals we imputed our target population from 777k to WGS SNPs, indels and deletions.

Evaluation of the Reference Panel

First, we prepared a phased reference panel using Beagle with the following inputs: (Phred-scaled) genotype likelihoods (PLs) for WGS SNPs and indels in 772 animals, deletion PLs in 175 animals from the study by Mesbah-Uddin et al. (2018a), deletion PLs estimated for 597 animals using auxiliary read depth data from VCF files with missing deletion PLs coded as "." for EM non-convergence. Second, to get the final reference panel, we re-phased and refined genotypes at loci with Beagle GP less than 0.99 using SHAPEIT. We refer to the first and second reference panels as RefBeagle and RefBeagle_SHAPEIT, respectively. We assessed the performance of this two-step approach on chromosome 29. The performance differences between the two reference panels, RefBeagle_SHAPEIT vs. RefBeagle, were tested using two-tailed paired t-test.

Cross validation.

We prepared two reference panels of chromosome 29: RefBeagle and RefBeagle_SHAPEIT.

The leave-out-trials were performed on 112 deletions with MAC ≥5 common in both dataset. For each deletion, we extracted WGS markers within 1 Mb on either side of the deletion, after excluding markers within and ±500 bp surrounding regions. Next, we removed the phasing for all the markers, and performed 77 leave-out-trials by randomly masking deletion genotypes for a set of non-overlapping (i.e. sampled without replacement) 10 animals in each trial (12 animals in the last trial). We used the same random seed to mask the genotypes of samples in both dataset. These analyses were performed using Beagle with the same parameters as in the "Phasing" step.

Here, we considered the deletion genotypes in the corresponding reference panel as the true genotype for that scenario. To evaluate the imputation accuracy at each deletion locus, we extracted the imputed and true genotypes for the corresponding animals from the 77 trials, and performed Pearson's correlation (r) of allele Chapter 4 72 dosages (coded as 0, 1, or 2 for homozygous reference, hemizygous or homozygous deletion, respectively)

between the reference and the imputed dataset. Finally, we reported the squared version of the correlation coefficient (r 2 ) in the Results section.

Downstream imputation performance.

For chromosome 29, we also assessed downstream performances of the two reference panels (RefBeagle and RefBeagle_SHAPEIT) to impute jointly the target 777k dataset to WGS SNPs, indels, and deletions using Minimac. Here, we evaluated the performance based on the imputation accuracy estimates (𝑣𝑣̂2) provided in the Minimac "info" output file. Under the assumption that genotype counts for a poorly imputed marker will shrink towards its

𝑔𝑔(1 -𝑔𝑔)

Here, 𝑔𝑔̂ = the estimated alternate allele frequency, 𝐶𝐶 𝑖𝑖 = the imputed alternate allele probability at the i th haplotype, and 𝑚𝑚 = the number of samples in the imputation target panel.

Tools used for data manipulation and analyses

We used following tools for data preparation and analyses: R software version 3.1.2 (R Core

Principal Component Analysis

To distinguish the three breeds present in our dataset, we performed a principal component (PC) analysis using the imputed WGS variants with MAF >1%. Within this MAF threshold, we randomly selected 1% SNPs, 10% indels, and 100% deletions from each of the 29 autosomes. The final dataset contained 137,064

SNPs, 11,710 indels, and 4,670 deletions. The analysis was performed using PLINK (v2.00a2 AVX2) software [START_REF] Chang | Second-generation PLINK: rising to the challenge of larger and richer datasets[END_REF]. For each of the three variant classes, we presented the first two PCs and along with the proportion of variance explained by the corresponding PC.

Results and discussion

Deletion Genotyping using RD data from the 1KBGP VCF file

The WGS reference panel was extended from 175 animals to 772 animals by data from the 1KBGP. In the extended reference panel, genotype likelihoods (PL) for 14,800,299 bi-allelic markers (SNPs + indels) on all the animals were available. However, the deletion PL for 5,798 selected markers were only available on RefBeagle: reference panel phased using Beagle; RefBeagleSHAPEIT: reference panel phased using Beagle followed by re-phasing using SHAPEIT 2 Performance gains between the two reference panels were calculated using R software function t.test(RefBeagle_SHAPEIT, RefBeagle, alternative = "two.sided", paired = TRUE, conf.int = TRUE)

3 WGS markers of chromosome 29: 107 deletions, 2,723 indels and 326,838 SNPs 4 The chosen Beagle GP threshold was 0.99 5 On each dataset, three different tests were performed; therefore, a P-value threshold of 0.017 (0.05/3) was selected at 5% level of significance, after multiple-testing correction 6 LoT: Leave-out-trials 7 Downstram imputation performance of the two reference panels to impute from the 777k to the WGS SNPs, indels and deletions of chromosome 29 using Minimac 8 The P-value <1.0 × 10 -16

Authors' contributions MMU, BG, MSL and GS conceived and designed the study. MMU performed computational analyses.

MMU drafted the manuscript. BG, MSL, DB and GS jointly supervised the study. BG, DB and GS critically revised the manuscript. All authors read and approved the final manuscript.

Chapter 5. we performed single-marker analyses followed by stepwise selection of independent markers using conditional and joint association (COJO) analyses.

Genome

From single-marker analyses for fertility traits, we reported genome-wide significant association (P-value <5×10 -8 ) of 30,384 SNPs, 178 indels and 3 deletions in HOL, 17 SNPs in JER, and 23,481 SNPs, 189 indels and 13 deletions in RDC. After COJO analyses of FI associated markers, we identified 37 and 23 independent associations in HOL and RDC, respectively. For these QTL region (top marker ±500 kb), we proposed several candidate genes that had functional annotations such as embryonic lethality, male and female infertility, oocyte degeneration, abnormal estrous cycle, decreased ovulation rate, and other fertility related phenotypes in mouse, cattle, and zebrafish.

Inclusion of these QTL information in genomic prediction is expected to provide gains in accuracy.

However, to maximize this gain, direct genotyping of these top markers could be considered by including them in the SNP array used for routine genotyping. Finally, validation of these QTL in an independent population should also be considered to ascertain the true nature of these associations.

Key words: GWAS, fertility, deletion, Holstein, Jersey, Nordic Red, dairy cattle

Introduction

Over the years, genome-wide association studies (GWASs) have provided a better understanding of underlying genetic architectures for many economical traits in dairy cattle. These include milk yield (Iso-Touru et al., 2016), milk compositions [START_REF] Gebreyesus | Multi-population GWAS and enrichment analyses reveal novel genomic regions and promising candidate genes underlying bovine milk fatty acid composition[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 Montbeliarde cows[END_REF], milking speed [START_REF] Jardim | Association analysis for udder index and milking speed with imputed whole-genome sequence variants in Nordic Holstein cattle[END_REF], Marete et al., 2018), clinical mastitis [START_REF] Kadri | Genetic dissection of milk yield traits and mastitis resistance quantitative trait loci on chromosome 20 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], body conformation (Abo-Ismail et al., 2017), stature (Bouwman et al., 2018), and so on. Fertility is also an important economic trait that was difficult to select before genomic selection due to its very low heritability. GWAS for female fertility provided an opportunity to detect quantitative trait loci (QTL) [START_REF] Hoglund | Validation of associations for female fertility traits in Nordic Holstein, Nordic Red and Jersey dairy cattle[END_REF], Hoglund et al., 2015a, Hoglund et al., 2015b) and to include these QTL in genomic prediction (Brøndum et al., 2015).

However, the accuracy of prediction for female fertility is lower compared with production traits. There is a potential to identify QTL or markers in strong linkage disequilibrium (LD) with QTL using GWAS for female fertility to utilize the QTL information in breeding value prediction. Two factors could be considered

Abstract

We performed within-breed genomic prediction analyses using genomic best linear unbiased prediction (GBLUP) method. Besides Bovine50k BeadChip (50k) markers, in GBLUP, we used one or several (genetic) variance components to investigate the effect of imputed whole-genome sequence (WGS) variants and large chromosomal deletions (deletions) on prediction accuracy of female fertility in Nordic Holstein (HOL), Danish Jersey (JER) and Nordic Red Dairy cattle (RDC) populations of artificial insemination bulls. In HOL, we observed same prediction accuracy (0.639) using either 50k (one-component model) or 50k and deletion markers (two-component model). In JER, we achieved a prediction accuracy of 0.664 with fertility index (FI) associated markers alone, which was ~10% higher than the prediction accuracy with 50k markers only (0.595). In RDC, we observed gains in prediction accuracy when imputed WGS variants as well as deletions were included in addition to 50k markers; we achieved a 2.1% gains in prediction accuracy in RDC with 50k and FI associated markers (in two-component model) vs. that of only 50k markers (one-component model).

Interestingly, we also achieved higher gains in prediction accuracy for FI in all three breeds compared with the previous reports. Finally, our results showed that imputed WGS variants and deletions could improve prediction accuracy. However, direct genotyping is necessary to explore the full potential of these markers as predictor in routine genomic evaluation program.

Key words: genomic prediction, large deletion, fertility, GBLUP, dairy cattle

Introduction

Estimation of genetic merits of animals using genomic markers (known as 'genomic selection') has revolutionized livestock breeding (Georges et al., 2019). Since the introduction of genomic selection in 2008 in US dairy cattle, the rates of annual genetic gain has increased by 50-100% on the breeding objectives (Garcia-Ruiz et al., 2016). This increase is attributed to a strong reduction in generation interval, reduced from 7 y to 2.5 y on the sire-daughter and sire-son pathways. For fitness traits, such as productive life and fertility, this relative increase is much larger, 3 to 4-fold in Holstein cattle, despite low heritability (0.04-0.08) of these traits (Garcia-Ruiz et al., 2016). This performance is due to the higher accuracy of genomic estimated breeding values for these traits. Nevertheless, prediction accuracy for breeding values is markedly lower than unity and can still be increased by various means, such as, larger reference population, improvement of statistical models, use of trait-associated markers from whole genome sequence variants and putative causal variants, use of other class of markers, besides SNPs and indels, like structural variants etc. (as discussed in Chapter 1 and Chapter 7 in this thesis).

Several studies investigated the effect of using the Bovine50k BeadChip (50k) markers along with few known quantitative trait loci (QTL), identified from imputed whole-genome sequence (WGS) variants, on prediction accuracy for female fertility (Brøndum et al., 2015, Su et al., 2015[START_REF] Ma | Use of a Bayesian model including QTL markers increases prediction reliability when test animals are distant from the reference population[END_REF]; however, very little gain in accuracy was reported for this trait. Structural variants (SVs) such as large chromosomal deletions (>50 bp, referred to as deletions, hereafter), on the other hand, has the potential to be causal, as Chapter 6 112 shown previously (Charlier et al., 2012, Kadri et al., 2014, Sahana et al., 2016), and thus when included as predictors could improve prediction accuracy, especially for fitness traits. To our knowledge, only one study investigated the effect of deletions and duplications along with other markers for predicting several growth traits in Nellore cattle (Hay et al., 2018). [START_REF] Liu | Using additional single nucleotide polymorphisms selected from whole genome sequence data for genomic prediction in Danish Jersey[END_REF] also reported SVs along with other WGS variants in their dataset; however, the primary focus of that study was to investigate the effect of training population size on prediction reliability. The objective of this study was to investigate the effect of imputed WGS selected markers, such as markers with suggestive association, and deletions, with or without 50k markers, on accuracy of genomic prediction for female fertility in HOL, JER and RDC.

Methods

Phenotype

We considered deregressed proofs of fertility index (FI) as the trait of interest for investing the improvement in accuracy of genomic prediction when various classes of DNA markers were used. This study was conducted in three Nordic dairy breeds, namely HOL (5,577), JER (1,211), and RDC (4,471) bulls with FI in our dataset.

Marker selection

We used imputed WGS variants, namely SNPs, indels, and deletions from our joint imputation study presented in Chapter 4 in this thesis. From this dataset, we selected 50k markers, deletions, and a set of markers associated with FI or its component traits that we identified in a within-breed GWAS (presented in Chapter 5). We only included markers with minor allele frequency (MAF) >0.01 and Minimac3 software (Das et al., 2016) imputation accuracy R 2 >0.1 in the training dataset. The GWAS markers were included when the P-value of association with FI or any of the five-component traits were <5×10 -5 in JER, while <10 - 7 in HOL/RDC. Five FI component traits were number of inseminations per conception in heifers (AISh) or in cows (AISc), interval (number of days) from calving to first insemination in cows (ICF), and interval (number of days) from first to last insemination in heifers (IFLh) or in cows (IFLc). Markers were selected within-breed; each selected marker was allocated to a class corresponding to a trait. When a single marker was associated with multiple traits within a breed, it was assigned to a trait according to the preference of FI >AISc >AISh >IFLc >IFLh >ICF. This trait preference was based on relative importance of a component trait on fertility index calculation as described in (NAV, 2013). A summary of the selected variants is presented in Table 6.1.

Prediction model

We performed within-breed prediction analysis using genomic best linear unbiased prediction (GBLUP) model, implemented in GCTA software (Yang et al., 2011), with a minimum of one-component variance model to a maximum of eight-component variance model. For each breed, we built eight genomic relationship matrices (GRM, or simply, G), each corresponding to a class of variants, G50k, GDeletions and one class for each trait (Table 6.1). We used the following GBLUP model:

Here, 𝒚𝒚 is the vector of phenotypes (de-regressed breeding values), 𝜇𝜇 is the phenotypic mean, n is the number of variant classes included in the model (𝑚𝑚 = 1, 2, … , 8); 𝒈𝒈 𝑖𝑖 is a vector of random genetic effects corresponding to i th variant class with 𝒈𝒈 𝒊𝒊 ~𝑁𝑁(𝟎𝟎, 𝑮𝑮 𝑖𝑖 𝜎𝜎 𝑖𝑖 2 ), Gi was constructed following VanRaden (2008) method-1 using GCTA software, 𝜎𝜎 𝑖𝑖 2 was the additive genetic variance explained by the i th variant class; and

𝜺𝜺 is a vector of residual effect with 𝜀𝜀~𝑁𝑁(𝟎𝟎, 𝑰𝑰𝜎𝜎 𝜀𝜀 2 ) where 𝑰𝑰 is an identity matrix and 𝜎𝜎 𝜀𝜀 2 is the residual variance.

Differences in accuracy of deregressed proofs were neglected. Variance components were estimated in GTCA software by restricted maximum likelihood (REML) approach with a combination of expectationmaximization (EM-REML) and average-information (AI-REML) algorithms (Yang et al., 2011).

For each breed, the GBLUP models were trained on animals born before the year 2005 (4,275 HOL, 1,031 JER and 3,626 RDC), and tested on animals born in the year 2005 or later (1,302 HOL, 180 JER and 845 RDC).

Results and Discussion

In a within-breed genomic prediction using the GBLUP method, we investigated the effect of a group of selected markers, each with a separate variance component (Table 6.1), on the accuracy of prediction for female fertility index in populations of HOL, JER, and RDC bulls. The average reliability of FI in our dataset was (train vs. test) 0.729 vs. 0.645 in HOL, 0.635 vs. 0.587 in JER, and 0.811 vs. 0.654 in RDC. Additive genetic variance and the proportion of phenotypic variance explained by different markers set is presented in Table 6.2. We present the prediction accuracy in Figure 6.1 and prediction bias in Table 6.3. In Table 6.4, we report reliability of predicted breeding values (without scaling) for the ease of comparisons with previous results on the same population.

In HOL, the prediction accuracy was highest (Figure 6.1a) and prediction bias (Table 6.3) was lowest when 50k markers in a one-component or 50k and deletion markers in a two-component GBLUP model were used.

However, when 50k and FI associated markers were used in a two-component GBLUP model prediction accuracy decreased by 0.023 compared with that of 50k, or 50k and deletion markers. Furthermore, using FI and its component traits associated markers only (without 50k) resulted in a ~31% relative reduction in prediction accuracy (Figure 6.1a) and ~61% increase in prediction bias where bulls' breeding values were overestimated (Table 6.3). Higher reliability of prediction with 50k markers was also reported previously for this breed (Brøndum et al., 2015). However, using 50k (or 50k with deletion) markers we had a gain of 36% ( 0.408 0.645 in Table 6.4 vs. 0.402 in Brøndum et al. (2015) Table A1) in reliability of prediction compared to that of Brøndum et al. (2015), with similar prediction bias (0.939 in this study vs. 0.943 in Brøndum et al.

(2015) Table A4). Interestingly, contrary to Brøndum et al. (2015), prediction reliability of 50k with GWAS selected markers was relatively lower in our study compared with that of 50k markers only. In JER, using only FI associated markers we observed a substantial gain in prediction accuracy (and reliability), with a regression slope (true vs. predicted) close to unity (Figure 6.1b and Table 6.3 and 6.4).

Using only 50k markers or 50k and deletion markers reduced prediction accuracy by 0.069 and 0.070, respectively, compared with that of FI associated markers only (Figure 6.1b). Although predictions with both one-component model of 50k markers and two-component model of 50k and deletion markers were biased upwards (breeding value were underestimated), the latter model was relatively less biased (Table 6.3); although deletion component had nearly zero variance in this model (Table 6.2). Importantly, using 50k or FI associated markers only we had a relative gain of ~12% (0.354 in this study vs. 0.311 previously reported) and ~29% (0.441 in this study vs. 0.311 previously reported) in prediction reliability of FI compared with previous reports on the same population [START_REF] Thomasen | Reliabilities of genomic estimated breeding values in Danish Jersey[END_REF], Su et al., 2015, Liu et al., 2018).

Chapter 7. General Discussion

variants, besides SNPs and small indels, like SVs for large-scale genome-wide screening.

Consequently, to our knowledge, these procedures are not yet included in any routine evaluation system because SVs lack flexibility and fully automatic pipelines.

In this PhD, we devoted two studies in mapping and imputing large deletions using WGS:

 In Chapter 3, we showed a systematic approach for high-resolution mapping of large deletions from WGS data of 175 animals, and subsequent genotyping and validation of the discoveries. Using deletion genotypes, we illustrated the population-genetic properties such as population diversity, stratification and differentiation. For essential gene deletion, we could only observe heterozygotes in our dataset. While for nonessential genes, homozygotes were observed among live animals. Using QTL enrichment analysis, we showed that deletion could be better predictor for health and fertility related traits; this was our hypothesis for performing subsequent imputation, GWAS and prediction studies. This study has made a catalogue of more than 8,000 large deletions segregating in three Nordic dairy cattle breeds. This catalogue of deletions will facilitate further research for mapping gene, improving accuracy in genomic selection and population genetics study.

 In Chapter 4, we described a Gaussian mixture model-based approach to estimate deletion-genotype likelihood using SNP read-depth data from VCF file. This was needed for enlarging the haplotype reference panel for imputing deletions into our bull population. There was a possibility of including ~600 additional animals (from HOL, JER and Norwegian RDC) from Run-6 of the 1000 Bull Genome Project, only if we could have deletion genotypes for these animals. However, access to the raw sequence (BAM) file is limited in this case. Nevertheless, read-depth data supporting each SNP genotype call present in VCF file is an indirect representation of number of sequencing reads mapped to the SNPs neighbouring region. Furthermore, number of reads aligned to a position is proportional to the copy-number, and there is a linear relationship between number of reads and copy-number of a region. We used sequencing depth from VCF as proxy and applied a similar genotyping approach that is used for SNP and SVs with raw sequencing data. We validated our approached using the known carrier status for Brachyspina deletion among the sequenced bulls. This technique has several limitations, to name few, it can only genotype copy-number loss or gain, will perform purely in regions with shallow coverage, fully depends on the presence of SNPs.

Chapter 7 123 Nevertheless, using this deletion genotyping approach, we were able to include additional WGS animals. Next, we built haplotype reference panel and subsequently imputed our array-typed population for SNPs, indels and deletions. Imputation accuracy for SNPs and indels were higher than in previous reports. And imputation accuracy for deletions were also high. Our results showed that common deletions could be imputed jointly with SNPs and indel with comparable accuracy. We used this imputed genotypes to perform GWAS and genomic prediction analyses.

Genomic Prediction

In dairy cattle, genomic selection provided substantial genetic gains, primarily by reducing the generation interval, for traits with low, moderate to high heritability (Garcia-Ruiz et al., 2016, Georges et al., 2019). Indeed, in contrast with phenotypic selection, high accuracy of genomic breeding values can be achieved with an appropriate reference population: many individuals with both phenotypic and genomic information can compensate for a low heritability. In addition, in its early stage in dairy cattle, genomic selection has been able to take advantage of a tremendous already existing resource, i.e. the large number of progeny tested bulls with very accurate breeding values even for low heritability traits. In this particular situation, the parameter of interest is not the heritability but the reliability of the evaluation, and the estimated breeding value based on progeny information is equivalent to an own performance of the bull with a heritability equal to the reliability.

However, the accuracy of genomic prediction has not reached to the label that was usually obtained through progeny testing, especially for populations with limited reference population size. Further increase in the accuracy of genomic prediction is necessary to achieve higher genetic gain. Several options have been studied in the literature in order to increase accuracy of genomic evaluation:

-Increase size of the reference population. This may be achieved (1) by genotyping more individuals with phenotype. However, this option is valid only if genotyping is cheap (or not directly supported by the breeding program) and if the population is large enough. Because the most informative individuals were genotyped first (i.e. the progeny tested bulls), the additional individuals are less informative and many more need to be genotyped to observe a significant gain. In addition, if the heritability of the trait is low such as for fertility, this option requires many individuals and may be unrealistic. Increasing size of the reference population may be also achieved by (2) building large consortia within breed such as EuroGenomics in Holstein (Lund et al., 2011, Boichard et al., 2018) or Intergenomics in Brown Swiss (https://www.brown-swiss.org/genetics), or by exchanging information across countries in a bilateral way, such as between Denmark and USA in Jersey [START_REF] Wiggans | Short communication: Improving accuracy of Jersey genomic evaluations in the United States and Denmark by sharing reference population bulls[END_REF]. This is an easy and cost-effective way to proceed, but it is possible only for international breeds and when 𝐺𝐺 × 𝐸𝐸 interactions are limited. (3) Many efforts have been dedicated to the use of a multi-breed reference population (Calus et al., 2018, van den Berg et al., 2019). Theoretically, this approach is very appealing because it reduces the fixed cost of the reference population in each breed.

Its success relies on different assumptions, e.g. a strong common basis of the genetic determinism (i.e. many common QTL across breeds) and a strong conserved linkage disequilibrium phase across breeds. In spite of all these efforts, results are not fully convincing and significant gains in accuracy have been obtained only for closely related breeds and for traits with high heritability. Due to the major associated stakes, these efforts are continued based on sequence data and methods able to extract the appropriate information.

-Feature Selection. Ideally, inclusion of causal variants as predictor will provide maximum prediction accuracy, as seen in simulation study (Meuwissen and Goddard, 2010). Inclusion of markers in high LD with the causal variants will also improve prediction accuracy (van den Berg et al., 2016a, van den Berg et al., 2016b). One extreme example of this approach is the prediction of stature in cattle using only 163 highly predictive markers: multi-breed metaanalysis to map QTL and subsequent selection of markers (Bouwman et al., 2018). Here, homogeneity of genetic architecture for the trait across breeds is vital, i.e. if different sets of QTL affect the trait of interest in different breeds, this will not work as expected. This point remains an open question mark. Indeed, many QTL, especially those with moderate effects, are not found to be shared across populations. But most mapping designs have limited detection power for these small QTLs. Therefore, it is difficult to conclude whether the genetic determinism is different or if it is the same but we are not able to show that.

Nevertheless, the reported gains in prediction accuracy for fertility traits were negligible at best, although the gains were substantial for production traits [START_REF] Brøndum | Reliabilities of genomic prediction using combined reference data of the Nordic Red dairy cattle populations[END_REF], Brøndum et al., 2015). This could be a circular-problem: it is hard to map fertility associated QTL with the available mapping population, and thus inclusion of less informative markers rather adds noise to prediction models and often resulted in reduction in accuracy. Noteworthy, prediction Chapter 7 125 performance with the same mapping population is relatively higher for traits with high heritability.

However, in this PhD thesis, we have demonstrated that improving imputation accuracy of WGS variants (Chapter 4) and subsequently, selecting fewer predictive markers from GWAS (Chapter 5) can enhance accuracy in genomic prediction for female fertility (Chapter 6):

 In Chapter 5, we described single-marker association analyses for eight fertility traits using imputed SNPs, indels and deletions from Chapter 4. After the initial GWAS scan, we performed conditional and joint analyses. We detected several QTL regions, some were previously reported, and some new discoveries. Besides SNPs and indels, several deletions passed the GWAS significant threshold of <5×10 -8 . However, only BTA12 deletion, which was previously reported by Kadri et al. (2014), was among the top signals. Nevertheless, association p-values for majority of the significant loci in our results were stronger (small) compared with previous report, which could be due to the higher accuracy of our imputed dataset. Overall, our results highlights the importance of including SVs for association mapping in livestock species.

 In Chapter 6, we presented the results of genomic prediction analysis. Here, we investigated the impact of imputed WGS variants on accuracy of genomic prediction.

We investigated several scenarios, such as, use of 50k data only, imputed deletion, fertility index associated markers, and a combination of these. We used GBLUP method with one or several genomic-relationship matrices derived from these marker sets. This method allows considering several classes of markers with different variance priors, while keeping a good simplicity of implementation. We observed substantial gain in prediction accuracy for fertility compared with previous reports. Our results demonstrated that deletions in our imputed dataset were also good predictors for fertility.

Evolutionary Conservation as a Tool for Identifying Lethal Genes

An interesting observation from three of the mapping studies (Chapter 2, 3 and 5) is that higher organisms share a common set of essential genes. These genes are indispensable for the organism's survival. On the other hand, higher organisms also share a variable set of genes that are dispensable due to redundancies in the genome. In fact, this is the basis of evolutionary The overall aim of this PhD thesis is to identify causal variants for recessive lethal mutations and select a set of predictive markers that are in high linkage-disequilibrium with the causal variants for female fertility in dairy cattle. We addressed this broad aim under five articles: (i) describes a systematic approach of mapping recessive lethals in French Normande cattle using homozygous haplotype deficiency (HHD). This study shows the influence of sample size, quality of genotypes, quality of (genotype) phasing and imputation, age of haplotype (of interest), and last but not the least, multiple testing corrections, on discovery and replicability of HHD results. It also illustrates the importance of fine-mapping with pedigree and whole-genome sequence (WGS) data, (cross-species) integrative annotation to prioritize candidate mutation, and finally, large-scale genotyping of the candidate mutation, to validate or invalidate initial results. (ii) describes a high-resolution populationscale mapping of large chromosomal deletions from whole-genome sequences of 175 animals from three Nordic dairy breeds. This study employs three different approaches to validate identified deletions. Next, it describes population genetic properties and functional importance of these deletions. (iii) deals with three main issues related to imputation of structural variants, in this case, large chromosomal deletions, e.g. availability of deletion genotypes, size of haplotype reference panel, and finally, imputation itself. To address the first two issues, this study describes a Gaussian mixture model-based approach where read-depth data from the variant call format (VCF) file is used to genotype a known deletion locus, without the need for raw sequence (BAM) file. Finally, it presents a pipeline for joint imputation of WGS variants along with large chromosomal deletions. (iv) describes genome-wide association studies for female fertility in three Nordic dairy cattle breeds using imputed WGS variants including large chromosomal deletions. This study is based on the analyses of eight fertility related traits using single-marker association, conditional and joint analyses. This study illustrates that inflation in association test-statistics could be seen even after correcting for population stratification using (genomic) principal components, and relatedness among the samples using genomic relationship matrices; however, this was known for traits with strong polygenic effects, among other factors. Finally, mapping of several new quantitative trait loci (QTL), along with the previously known ones, are reported in this study. This study also highlights the importance of including (imputed) large deletions for association mapping of fertility traits. (v) describes prediction of genomic breeding values for fertility using SNP array-chip genotypes, selected QTL and large chromosomal deletion. Using genomic best linear unbiased prediction (GBLUP) method with one or several genomic-relationship matrices derived from a set of selected markers, this study reports higher prediction accuracy compared with previous report. This study also highlights the influence of selecting markers with best predictability, especially for a breed with small training population, in accuracy of genomic prediction. The results demonstrate that large deletions in general have a high predictive performance. L'objectif général de cette thèse est d'identifier les variants causaux ou, à défaut, un ensemble de marqueurs prédictifs -qui présentent un déséquilibre de liaison élevé avec les variants causaux -pour la fertilité des vaches laitières. Nous avons abordé cet objectif général dans cinq articles: (i) décrit une approche systématique de cartographie des variants létaux récessifs chez les bovins Normands français basée sur la recherche de déficit en haplotypes homozygotes (HHD). Cette étude montre l'influence de la taille de l'échantillon, de la qualité des génotypes, de la qualité du phasage des génotypes en haplotypes et de l'imputation, de l'âge de l'haplotype et enfin, de la définition des seuils de signification prenant en compte les tests multiples, sur la découverte et la reproductibilité des résultats de HHD. Elle illustre également l'importance de la cartographie fine avec les données de généalogie et de séquence de génome entier (WGS), l'annotation intégrative (entre espèces) pour hiérarchiser les mutations candidates et, enfin, le génotypage à grande échelle de la mutation candidate, pour valider ou invalider les mutations initiales. (ii) décrit une cartographie à haute résolution de grandes délétions chromosomiques de séquences du génome dans une population de 175 animaux appartenant à trois races laitières nordiques. Cette étude utilise trois approches différentes pour valider les résultats de la cartographie.

Individual Training Plan

Le chapitre décrit les propriétés génétiques des populations et l'importance fonctionnelle des délétions identifiées. (iii) traite de trois questions liées à l'imputation de variants structuraux, ici de délétions chromosomiques importantes: la disponibilité des génotypes de délétion, la taille du panel de référence d'haplotypes et, enfin, l'imputation elle-même. Pour aborder les deux premières questions, cette étude décrit une approche basée sur un modèle de mélange gaussien dans laquelle les données de profondeur de lecture provenant de fichiers au format VCF (variant call format) sont utilisées pour génotyper un locus de délétion connu, en l'absence d'information sur la séquence brute. Enfin, il présente un pipeline pour l'imputation conjointe de variants WGS et de grandes délétions chromosomiques. (iv) décrit des études d'association pangénomiques de la fertilité femelle dans trois races de bovins laitiers nordiques à l'aide de variants WGS imputés et de grandes délétions chromosomiques. Cette étude concerne huit caractères de fertilité et utilise des analyses d'association mono-marqueur, conditionnelles et conjointes. Cette étude montre qu'une surestimation, ou « inflation », des statistiques de test peut être observée même après correction pour la stratification de la population à l'aide de composantes principales génomiques et pour les structures familiales à l'aide de matrices de relations génomiques. Ce biais était connu pour les caractères très polygéniques. Enfin, cette étude présente plusieurs locus de traits quantitatifs (QTL) nouveaux et confirme plusieurs autres déjà connus. Elle souligne également l'importance d'inclure les grandes délétions (imputées) pour la cartographie par association des caractères de fertilité. (v) décrit la prédiction des valeurs génomiques de fertilité (ou indice de fertilité) à l'aide de génotypes à puces SNP, de QTL sélectionnés et de délétions chromosomiques importantes. En utilisant la méthode de meilleure prédiction linéaire sans biais génomique (GBLUP) avec une ou plusieurs matrices de relations génomiques dérivées d'un ensemble de marqueurs sélectionnés, cette étude rapporte une précision de prédiction améliorée. Cette étude met également en évidence l'influence de la sélection des marqueurs les plus prédictifs, en particulier pour une race ayant une population d'apprentissage réduite, sur la précision des prédictions génomiques. Enfin, les résultats démontrent que les grandes délétions ont en général un pouvoir prédictif élevé.