# Uncovering Adaptation from Sequence Data: Lessons from Genome Resequencing of Four Cattle Breeds

Simon Boitard, Mekki Boussaha, Aurelien Capitan, Dominique Rocha, Bertrand Servin

<figure><img src="image_1.png" /></figure>

## To cite this version:

Simon Boitard, Mekki Boussaha, Aurelien Capitan, Dominique Rocha, Bertrand Servin. Uncovering Adaptation from Sequence Data: Lessons from Genome Resequencing of Four Cattle Breeds. Genetics, 2016, 203 (1), pp.433-450. 10.1534/genetics.115.181594. hal-01533825

HAL Id: hal-01533825

https://hal.science/hal-01533825v1

Submitted on 6 Sep 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
Uncovering adaptation from sequence data: lessons from genome
resequencing of four cattle breeds

Simon Boitard¹,², Mekki Boussaha¹, Aurélien Capitan¹,³, Dominique Rocha¹, and Bertrand Servin⁴

¹GABI, INRA, AgroParisTech, Université Paris-Saclay, Jouy-en-Josas, 78350, France

²Ecole Pratique des Hautes Etudes, UMR 7205 ISYEB (Muséum National d'Histoire Naturelle & CNRS & EPHE & Université Pierre et Marie Curie), Paris, 75005, France

³Union Nationale des Coopératives d'Elevage et d'Insémination Animale, Paris, 75595, France

⁴INRA, UMR1388 GenPhySE (INRA & Institut National Polytechnique de Toulouse), Castanet-Tolosan, 31326, France
11 **Running title:** Detecting adaptation in cattle using resequencing data

12 **Keywords:** FST, domestication, linkage disequilibrium, next generation sequencing, selective sweeps.

13 **Corresponding author:** Simon Boitard

14 Mailing adress: INRA GABI, Domaine de Vilvert, 78352 Jouy-en-Josas, France.

15 Phone: 0033 1 40 79 81 66

16 E-mail: sboitard@mnhn.fr
Detecting the molecular basis of adaptation is one of the major questions in population genetics. With the advance in sequencing technologies, nearly complete interrogation of genome-wide polymorphisms in multiple populations is becoming feasible in some species, with the expectation that it will extend quickly to new ones. Here, we investigate the advantages of sequencing for the detection of adaptive loci in multiple populations, exploiting a recently published dataset in cattle (*Bos taurus*). We used two different approaches to detect statistically significant signals of positive selection: a within-population approach aimed at identifying hard selective sweeps and a population-differentiation approach that can capture other selection events such as soft or incomplete sweeps. We show that the two methods are complimentary in that they indeed capture different kind of selection signatures. Our study confirmed some of the well known adaptive loci in cattle (e.g. MC1R, KIT, GHR, PLAG1, NCAPG/LCORL) and detected some new ones (e.g. ARL15, PRLR, CYP19A1, PPM1L). Compared to genome scans based on medium or high-density SNP data, we found that sequencing offered an increased detection power and a higher resolution in the localization of selection signatures. In several cases, we could even pinpoint the underlying causal adaptive mutation, or at least a very small number of possible candidates (e.g. MC1R, PLAG1). Our results on these candidates suggest that a vast majority of adaptive mutations are likely to be regulatory rather than protein coding variants.

Introduction

Detecting the molecular basis of adaptation in natural species is one of the major questions in popula-
tion genetics. With the spectacular progress of genotyping and sequencing technologies, genome-wide
scans for positive selection have been performed in multiple species and populations within the last
decade. Livestock species provide a considerable ressource for these selection scans, because they have
been subjected to strong artificial selection since their initial domestication, leading to a large variety
of breeds with disctinct morphology, coat color or specialized production. Besides, the economic value
of these species and the need to improve them has motivated the developement of standardized single
nucleotide polymorphism (SNP) chips and the genotyping of millions of animals using these chips,
providing considerable data for population genetics analyzes. For instance, in taurine cattle, at least
21 genomic scans for selection have already been published and were reviewed in Gutierrez-Gil et al.
(2015). Numerous genomic scans for selection have also been published in other livestock species, see
Gouveia et al. (2014) for a review.

The regions detected by these studies are generally convincing, because they contain interesting
positional and functionnal candidate genes (e.g. Fariello et al. (2014)) and / or are statistically
enriched with genes from regulation pathways related to production traits (e.g. Flori et al. (2009)).
Nevertheless, these regions often span several Megabases and typically include up to tens of genes, so
determining the exact gene(s) under selection in each region, and even more the causal mutation(s),
remains difficult from these studies. This might change with the recent advent of next generation
sequencing technologies, which allows to characterize a very large proportion of the variants in the
genome of a species. However, although several examples of genomic scans for selection based on large
sequencing samples are already available in livestock species such as pig (Rubin et al., 2012; Li et al.,
2013), cattle (Qanbari et al., 2014) or chicken (Rubin et al., 2010; Roux et al., 2015), the advantage of using sequencing data for detecting selection signatures has still not been widely discussed.

Another interesting question is the small overlap observed between studies, even when these studies focus on similar populations (Biswas and Akey, 2006; Qanbari et al., 2011). This can largely be explained by the fact that many different detection methods have been applied, whose sensitivity depends on the type of selection event: recent or old, complete or ongoing, from a new variant (i.e. “hard”) or standing variation (i.e. “soft”) etc. These differences between methods have been described by several studies (Biswas and Akey, 2006; Sabeti et al., 2006; Gouveia et al., 2014), but the practical implications of these differences when comparing the regions detected by different studies, or by different approaches within a same study, are rarely discussed.

Here we detected selection signatures in 4 European taurine cattle breeds (Angus, Fleckvieh, Holstein and Jersey), using large samples of sequencing data which have been recently published by the 1000 bull genomes project (Daetwyler et al., 2014). We used two different statistical approaches: a within-breed approach detecting genomic regions with low genetic diversity (Boitard et al., 2009), and a between population approach detecting genomic regions with large allele (Bonhomme et al., 2010) or haplotype (Fariello et al., 2013) frequency differences between breeds. These two analyzes provided regions whose genomic features are significantly different from what would be expected under neutrality, even when accounting for the effects of past population size changes, population structure and gene flow on the four breeds of our study. We show that applying the above methods to sequencing data improves the detection power of selection signatures and reduces considerably the length of detected regions. In some particular situations, it even leads to identify the exact mutation under selection. We also provide a detailed characterization of the regions that are detected only by the within-population approach (in one or several populations), only by the between population approach, or jointly by the two approaches.

# Results

## Genomic regions with low within-population diversity: hard sweeps signatures

We looked for hard sweep signatures within each breed using the method of Boitard et al. (2009), as described in the Methods. This method detects regions showing an excess of low and high frequency derived alleles compared to the rest of the genome. Although this pattern is typically expected under a hard sweep sscenario, i.e. when a new mutation appears in the population and goes to fixation due to positive selection, it may also arise from purely demographic events, in particular bottlenecks. In order to test if our analysis could be influenced by such false positive signals, we first simulated genomic samples under two different neutral demographic models allowing to reproduce the genetic diversity of the breeds under study, and applied the method of Boitard et al. (2009) to these samples.

### Analysis of neutral samples
In the first demographic model, we considered the joint history of the four breeds and accounted for several important features of this history: (i) the shared ancestry
of the four breeds, which diverged recently from an ancestral pool of european domestic animals, (ii) the population size differences between breeds since their divergence, and (iii) the possible existence of gene flow between breeds. Moreover, because recent studies suggested that effective population size in taurine cattle strongly declined since domestication (MacLeod et al., 2013; Boitard et al., 2016), we allowed one population size change in the ancestral population. We estimated the parameters of this model from the joint allele frequency spectra observed in our data for all breed pairs, using the approach of Excoffier et al. (2013) (see the methods for more details), and obtained the demography shown in Figure 1.

In this estimated demography, the ancestral population size change was not a decline related to domestication, but an older expansion occuring in the wild population, approximately 120,000 years before present. However, a very strong population decline was found at the time where the four breeds diverged, from an order of 100,000 individuals to an order of 100 individuals. Interestingly, the estimated divergence time (500 years before present) was consistent with a geographic isolation process starting a few hundred years before the strict separation of these populations, induced by the creation of modern breeds (Felius et al., 2011). Besides, the order of magnitude of estimated recent effective sizes (100) and the ranking of breeds according to these sizes were consistent with previous studies (Leroy et al., 2013; The Bovine HapMap Consortium, 2009; MacLeod et al., 2013; Boitard et al., 2016). Thus, this simple model seemed to provide a reasonable approximation of the demography of the four breeds under study. When genomic samples were simulated from this model, the average number of sweeps detected per genome was equal to 0 in Fleckvieh and Holstein, 3.75 (± 1.93) in Angus and 17.25 (± 4.15) in Jersey.

We performed the same test using a second demographic model, which was estimated in another study (Boitard et al., 2016) based on the same dataset as that considered here. This model focuses on each breed independently of the others, so it does not account for shared ancestry or gene flow. However, it accounts for the variations of population size over time more accurately than our first model, because population size in each breed is modelled as a stepwise process with 21 time windows (Figure 6 in (Boitard et al., 2016)). The population size within each time window is estimated from the allele frequency and linkage disequilibrium patterns observed in the breed, using an approximate Bayesian computation approach. When genomic samples were simulated from this second model, the average number of sweeps detected per genome was even lower than with the first model: 0 in Fleckvieh, Holstein and Jersey, and 0.75 (± 0.87) in Jersey.

Overall, these results indicate that the number of false hard sweep signals detected by the method of Boitard et al. (2009) should be negligible in Angus, Fleckvieh and Holstein, and relatively small in Jersey, even when accouting for the demography of these breeds.

**Overview of the detected signals** When analyzing the cattle data with the same approach, we detected 1,057 hard sweep signals: 226, 384, 316 and 131 in Holstein, Angus, Jersey and Fleckvieh respectively. According to the simulation results presented above, the false discovery rate associated with this analysis should be below 7% in Jersey (21.4/316), and close to 0 in the other breeds. The size of detected regions ranged from 8.2kb to 948kb, with a median of 78.7kb. Some signals were overlapping between breeds so that after merging them we obtained 798 sweeps that were unique to
one of the breeds (159, 297, 249 and 93 respectively) and 118 that were shared between at least two
breeds. Overall this provided 916 regions covering about 4.3% of the autosomal genome. Among the
916 regions covered by hard sweep signatures, 450 included no (protein coding) gene, 268 included a
single gene, 154 included between 2 and 5 genes, and 44 included more than 5 genes, with a maximum
of 19 genes. Overall, 1088 genes were included in sweeps windows, which represents about 5.7% of
all annotated genes in the bovine genome, so there was a slight enrichment of protein coding genes
within sweep regions. The list of all detected regions and of genes included in these regions is given
in File S2.

In a recent genome scan for selection focusing on the Fleckvieh breed (Qanbari et al., 2014), the 43 Fleckvieh sequences considered in our study were analyzed using the CLR method (Nielsen et al., 2005) and the iHS method (Voight et al., 2006). Seventy-three hard sweep signals were found with the former approach and 67 with the latter. Since the HMM approach used in this study aims at capturing the same allele frequency patterns than the CLR method, we checked if our results in Fleckvieh were consistent with those in Qanbari et al. (2014). To this end, we compared the distributions of CLR p-values within regions associated with selective sweeps to their distribution on the rest of the genome. Figure 2 (left) plots the ratio of the two densities (on a log2 scale) for increasing levels of significance of the CLR test. It shows a very strong enrichment of low CLR p-values in the sweep windows we detected in Fleckvieh, compared to the rest of the genome. We performed the same analysis with iHS p-values and found a similarly strong enrichment (Figure S6), which can be explained by the fact that iHS also tries to detect hard sweep patterns, even if the information used (the length of haplotypes) is different.

More surprisingly, we also observed an enrichment, albeit of lower intensity, of CLR and iHS low p-values in the sweep windows detected in other breeds than the Fleckvieh. For example, Figure 2 (right) shows the enrichment in low Fleckvieh CLR p-values in selective sweep regions detected only in the Holstein breed. A similar trend was observed in selective sweep regions specific to the Jersey and Angus breeds (not shown). Hence some of the hard sweeps detected in one breed also have probably taken place in the other breeds, but to a slightly lower extent that did not lead to a significant signal. These signatures must be related to favorable alleles that either started to increase in frequency before the divergence of the breeds or were selected in parallel in different breeds. However, selection signatures that are specific to one breed are interesting because they illustrate the importance of this breed for cattle functional diversity (Gutierrez-Gil et al., 2015). We therefore derived a way of finding clear breed-specific sweeps as follows.

**Hard sweep regions specific to one population** The HMM approach “tags” a region as selected by reconstructing a hidden state at each position of the genome. Based on the observation above, we suspected that some regions were not tagged as selected, but might still have a non negligible probability of being adaptive under the HMM model, explaining the enrichment patterns observed in Figure 2. To investigate this possibility, we derived a statistic ($T_{OR}$) measuring the strength of evidence for selection in a breed, measured as the log odds ratio of selection over neutrality in a region (see methods for details). As expected, this statistic showed clearly different distributions in regions tagged as selected in one breed and in regions where no selection was detected in any breeds. In
regions where selection was detected in another breed, the distribution of $T_{OR}$ was skewed toward higher values compared to clearly neutral regions (Figure S9). We exploited this to call breed-specific sweeps regions where $T_{OR}$ was unambiguously consistent with the neutral density in all other breeds (see methods for details). 55 breed-specific regions were detected, and we could check that this time the sweeps specific to Holstein did not show any enrichment in low Fleckvieh CLR p-values (Figure S10). The 12 sweeps exhibiting the most contrasted patterns for $T_{OR}$ are listed in Table S2.

### Hard sweep regions shared by all populations
We also looked for sweep signals shared by all breeds, as they might correspond to older selection events, possibly related to initial cattle domestication and at least anterior to the divergence of the 4 breeds considered here. We found only one region with a sweep detected in all 4 breeds, but we also considered regions where a sweep was detected in 3 breeds and where allele frequencies in the fourth one were almost consistent with a sweep. This provided 8 candidate regions, 4 of which include a single gene (Table 1).

Several of these genes represent natural detection targets in cattle, as they are related to husbandry, metabolism or fertility. On chromosome 1, we found evidence for selection in a region 10Kb upstream the OLIG1 gene, encompassing a lincRNA, orthologous to the Human gene LINC00945, which expression has been shown to be associated with polledness in Holstein and Fleckvieh (Allais-Bonnet et al., 2013). PPM1L is a protein phosphatase that has been shown to be involved in the response to exercise in humans (Tonevitsky et al., 2013); SLC25A33, located in the middle of one of the shared sweeps windows, encodes for a mitochondrial pyrimidine nucleotide transporters and is essential for mitochondrial DNA and RNA metabolism in humans (Di Noia et al., 2014). CYP19A1 encodes the key enzyme for estrogen biosynthesis. Many studies have documented its role during the development of bovine follicles, and it has been found more abundant in bovine cells of twinters versus controls (Echternkamp et al., 2012). To illustrate the allele frequency patterns observed in such regions, allele frequencies in the polled locus region are provided in Figure 3.

### Genomic regions with large genetic differentiation between breeds.

We applied two approaches to detect genome regions that exhibit outlying divergence in single sites (Bonhomme et al., 2010) or haplotype (Fariello et al., 2013) frequencies between populations. The population tree of the 4 breeds (Figure 4) was approximately star shaped, (i.e. breeds essentially evolved in parallel from an ancestral population), although we estimate a small shared history between the Fleckvieh and Jersey populations. Fixation indices, represented by the branch length from the root to the tips of the tree, had similar values in all populations, the Jersey's one being slightly higher. The genome-wide level of differentiation between populations was rather large, with between population $F_{ST}$ values ranging from 0.22 to 0.4.

When genetic drift is large, single marker tests are expected to have low power because even large allele frequency differences can be explained by drift alone. This is indeed the case here for the single marker-based differentiation analysis (FLK): the smallest observed p-value was $6.10^{-7}$, which corresponded to an FDR of about 10% when applying the approach of Storey and Tibshirani (2003). Although it is not clear how to correct p-values for correlation between markers in such a setting (genome-wide differentiation based tests), even this smallest p-value did not provide clear evidence of
being related to a selection signature. This does not mean that there is no selection in these data, only
that genuine selection signatures cannot be discriminated from background noise provoked by drift
when looking at sites independently. However, as illustrated below, given a region where a selection
signature has been found, FLK can help in identifying the mutations that have likely been under
selection.

To account for LD patterns when looking for selection signatures, we used a method, hapFLK
(Fariello et al., 2013), similar to FLK but that incorporates LD information through the exploitation
of a multilocus locus LD model (Scheet and Stephens, 2006). Because hapFLK combines information
across multiple sites, it has been shown to have better detection power than single site statistics
(Fariello et al., 2013, 2014). This was confirmed when applied to this dataset, as we could confidently
find 67 significant regions using a FDR threshold of 15%. Ten of these regions likely resulted from
assembly errors and were thus not considered in the rest of our analysis (Table S1). The cumulated
length of the remaining 57 regions, listed and annotated in table S7, was 9.1 Megabases, 0.36 % of
the total autosome length. Detected regions spanned from a few hundred basepairs to more than one
megabase with a median length of about 20 kilobases (Figure S7). The median size of detected regions
was thus considerably smaller than that obtained in previous studies where hapFLK was applied to
60K data in sheep (Fariello et al., 2014; Kijas, 2014). Nineteen of the regions encompassed at least one
gene while 38 contained no gene. In total, 82 genes were included in hapFLK regions, which represents
about 0.4% of the total number of genes in the genome, corresponding to a small enrichment in protein
coding genes in hapFLK signatures.

To investigate whether sequencing improves detection power with hapFLK, we thinned the dataset
by considering only SNPs present on the Illumina BovineHD BeadChip. We found, as shown in Figure
5, that the excess of small p-values was much larger when applying hapFLK to all sites identified in
the 1000 bull genomes project, than when applying it only to the SNPs that are included in the SNP
chip. Note that this was not the case with FLK, where detection power was low with both sequencing
or SNP chip data due to the amount of drift, as already discussed above (Figure S8).

Hard sweep regions showing a strong differentiation signal

Eight selection signatures were found with both the HMM and the hapFLK analyzes, pointing out
regions that combine a low diversity within at least one breed and a strong haplotypic differentiation
beween breeds (Table 2). Although this is significantly more than expected if the two kind of signals
were independent (p = 1.4 10⁻⁴), it is somewhat surprising to observe that many hard sweeps were
not detected with hapFLK.

A prevalent reason for this is the large genome-wide level of differentiation between the four breeds,
which reduces power of differentiation based tests (Fariello et al., 2013). Indeed, hard sweep regions
detected in one or two populations showed a clear enrichment in low hapFLK p-values (Figure 6).
This indicates that many hard sweeps exhibit a mild differentiation signal, although the power to
detect them with hapFLK is not sufficient. In addition, some hard sweep regions did not show any
differentiation signal, because the same haplotype was fixed or at least increased in frequency in all
populations. This is typically the case of hard sweep regions detected in three or four populations,
for which there was a depletion of low hapFLK p-values (Figure 6). This may also concern regions where an hard sweep was detected in one or two populations, but where the swept haplotype was also at quite high frequency in other populations, as already discussed above.

Among the regions with evidence for both a hard-sweep and an extreme differentiation signature, the top three corresponded to genes and mutations of known phenotypic effects that recapitulate the most obvious phenotypic divergence of the four breeds in this dataset. The most differentiated region corresponded to the KIT gene, which has been shown to be associated with white spotting patterns in the Holstein (Hayes et al., 2010) and the Fleckvieh (Qanbari et al., 2014), while the Jersey and the Angus are non spotted breeds. The next most differentiated region harbored the MC1R gene, for which previous studies have identified two causal polymorphisms for coat color (Klungland et al., 1995). Finally, the third signature was a small genomic region comprising the PLAG1 gene (Figure 7, top). Hard sweeps identified in this region indicate a past selection event affecting the Angus and the Holstein breeds, whereas no selection was evidenced in the Fleckvieh and the Jersey breeds (Table 2), as can also be seen from heterozygosity patterns in the region (Figure 7, bottom). The region surrounding PLAG1 was previously demonstrated to harbour a QTL for calving ease in the Fleckvieh (Pausch et al., 2011) and one for stature in a Holstein × Jersey cross (Karim et al., 2011). Our results are consistent with these studies and suggest that the allele favoring high stature was selected in Holstein and Angus, but not in Jersey and Fleckvieh.

The other common regions include the ASIP gene, which plays a crucial role in adipocyte development and seems to be expressed in a wide set of tissues in different cattle breeds (Albrecht et al., 2012) and the ANKRD55 region that is strongly associated with autoimmune disorders in humans (in particular multiple sclerosis (Stahl et al., 2010)) and type 2 diabetes (Morris et al., 2012).

## hapFLK signatures of soft or incomplete sweeps

Apart from the 8 signatures above, none of the other hapFLK signatures match hard sweeps detected by the HMM approach. These are most likely incomplete sweeps for which the selected allele did not reach fixation, or soft sweeps where selection targetted an allele that was already at intermediate frequency in the population. Indeed, such signals do not lead to the skewed allele frequency patterns that are looked for by the approach of Boitard et al. (2009). This hypothesis was confirmed by examining haplotype diversity patterns in hapFLK signatures, which typically show haplotypes of large but not fixed frequency in at least one population (e.g. Figure S11 for region 1 in table S7).

Two signatures are located close to the homologous region of a human promoter region near ROBO1 on one hand and the prolactin receptor gene PRLR on the other, hinting that the causal mutation is likely regulatory in nature. ROBO1 has been shown to be involved in early follicular development in sheep (Dickinson and Duncan, 2010; Dickinson et al., 2010), while the prolactin and its receptor are involved in a large range of biological functions (development, metabolism, immunology, reproduction ...) (Bole-Feysot et al., 1998).

The PRLR gene lies ≈ 6 Mb to the growth hormone receptor (GHR) gene, itself close (≈ 140Kb) to another hapFLK signature. It has been evidenced that two QTLs affecting milk, fat and protein yield segregate near these two genes in a Finnish dairy cattle population (Viitala et al., 2006). Our results
suggest that these QTLs might be in regulatory regions of these genes and that they have responded to selection in populations from the 1000 bull genomes. We found another selection signature, in Jersey, within a gene potentially involved in milk production (Figure S12): ARL15, seven highly differentiated variants (FLK p-value < 10⁻⁵) lying in an ARL15 intron. ARL15 is a protein of unknown function that has been shown to be strongly associated with adiponectin levels in humans (Richards et al., 2009). Adiponectin is an hormone involved in glucose metabolism, low concentration of adiponectin being associated with insulin resistance. In dairy cows, insulin resistance is maintained in early lactation, favoring mammary glucose uptake. Giesy et al. (2012) showed that the reduction of plasma adiponectin concentration in early lactating cows was not associated to changes in the adiponectin expression itself. Our results could suggest a potential role of ARL15 in this process, with a particular adaptation of the Jersey dairy cattle at this gene.

Apart from the PLAG1 signature, several other include genes involved in morphology and growth: RUNX3 (Yoshida et al., 2004; Soung do et al., 2007), STARD3NL (Rivadeneira et al., 2009) and RASSF2 (Song et al., 2012) are involved in bone development, NCAPG and/or LCORL match a large QTL for many growth traits in cattle, horse and sheep (Eberlein et al., 2009; Weikard et al., 2010; Lindholm-Perry et al., 2011; Signer-Hasler et al., 2012; Makvandi-Nejad et al., 2012; Bongiorni et al., 2012; Metzger et al., 2013; Tetens et al., 2013; Lindholm-Perry et al., 2013; Kijas, 2014; Xu et al., 2015; Randhawa et al., 2015; Sahana et al., 2015) and CTNNBL1 (Liu et al., 2008) is associated with obesity traits in humans.

## Identifying causal adaptive polymorphisms

As demonstrated above, genomic scans for selection based on sequencing data have a higher detection power than those based on genotyping chip data, and locate the selection signatures with higher precision. This is an expected outcome of the higher marker density, and the same could be said when comparing high density to medium density chip data. But a more fundamental difference between dense SNP chip data and sequencing data is that, with the latter, one can reasonably expect the causal polymorphism under selection to be included in the observed data. Clearly, not all selection signatures can be related to a single polymorphism, and even in this case this polymorphism might be absent in our data due to insufficient coverage or remaining calling issues. Still, the favorable situation where a single polymorphism leads to a selective advantage and is present in the data should also occur. One natural question is thus to determine if such variants can be identified only from genetic diversity patterns. We show below that this is indeed possible.

## HapFLK signals

Haplotype frequency differences detected by hapFLK typically result from the increase in frequency of one particular allele in a population due to positive selection, which implied the increase in frequency of one or several haplotypes carrying this allele by genetic hitchhiking. Thus, in regions detected by hapFLK, the causal polymorphism under selection should be the one with the largest allele frequency differentiation, and a natural strategy to detect this polymorphism should be to look at the variants with the largest FLK value. Two of the regions identified by the HMM and the hapFLK approach
listed in Table 2 validate this strategy.

Within the MC1R selection signature, we found three polymorphisms with clear outlying FLK values (Figure 8), and two of these correspond to the known causal mutations mentioned previously: a single base mutation at position 14757910 (rs109688013), responsible for the black pigmentation in Holstein and Angus breeds and a single base deletion at position 14757924 (rs110710422), responsible for the red pigmentation in the Fleckvieh breed (Klungland et al., 1995). The third outlying polymorphism (rs110494166) in the region was located at position 14678403, within an intron of the nearby FANCA gene, and exhibited the same allele frequencies as rs110710422.

In the PLAG1 region, the causal mutation was shown to be one of 8 candidate QTNs (Karim et al., 2011). We found 7 polymorphisms exhibiting a high level of differentiation between breeds based on the FLK statistic (p-value below $10^{-4}$) (Figure 7 middle, Table S4). Out of these 7 candidates, 6 were common with the candidate QTNs listed in the Table 2 of Karim et al. (2011). This indicates that the causal mutation for stature identified by Karim et al. (2011) has been under selection, and that this selection applied not only in Holstein but also in Angus. rs134215421 and rs109815800 showed particularly extreme allele frequency differences between Holstein and Angus on the one hand and Jersey and Simmental on the other hand. These two mutations are located at the 3' end of the PLAG1 gene, rs134215421 being 1Kb downstream and rs109815800 within and intron of the PLAG1 gene. One SNP (rs134029466) was not identified as a potential QTN in Karim et al. but showed a strong signal for differentiation. Two of the mutations in Karim et al., rs209821678 and rs210030313, which are considered by the authors as the most serious candidates because they affect a highly conserved element, were not available in the 1000 bull genomes project. One is a VNTR, a kind of polymorphism that is hardly callable from short-read sequence data and the other one is a SNP that lies 44 base pair from the VNTR, exhibited very low quality scores and was therefore not called in the 1000 bull genomes data. While these polymorphisms cannot be considered as disqualified based on our study, their positions, highlighted in Figure 7, lie in a region where the differentiation signal was not significantly elevated.

## Hard sweep signals

Hard sweep signals are expected when an allele goes from very low frequency to almost fixation in a population due to positive selection. In these regions, detecting the causal selected variant only from genetic data of the swept population is impossible, because all physical positions show either extreme allele frequencies or no polymorphism at all. However, if we assume that other sampled populations evolved neutrally in the region, alleles that were initially at low frequency in the swept population have likely remained at relatively low frequency in these other populations. This should result in high genetic differentiation at and around the selected polymorphism, which again can be detected using the FLK statistic.

For all hard sweep regions except the ones shared by all populations, we thus tried to identify the selected site by looking for polymorphisms with a high FLK value ($p \le 10^{-4}$). To ensure that the detected polymorphisms could indeed be the causal ones, we further required a high allele frequency ($\ge 75\%$) in the swept population(s), and only in those ones. We found only 12 sweep regions exhibiting
such causal candidates (Tables 3 and S5). Again, this small proportion is related to the fact that
in most cases the selected allele must actually be at relatively high frequency even in non swept
populations, due to undetected ongoing sweeps or just random drift.

The regions detected by this approach include those of KIT, PLAG1 and MC1R. In all these
regions, a very limited number of potentially causal polymorphisms was detected, and in the case of
MC1R these candidates included the true causal variants. This provides a validation of the detection
strategy considered here and suggests that other regions of Table 3 should also contain interesting
putative causal polymorphisms.

Of particular interest, a sweep region on chromosome 20 included a single causal candidate at
position 39872347 (Figures 9 and S13). This polymorphism is located downstream SLC45A2, a gene
that has been associated to fertility (Killeen et al., 2014) and residual feed intake (Karisa et al., 2013)
in cattle, and to pigmentation in several other species including human (Sturm, 2009; Stefanaki et al.,
2013; Morice-Picard et al., 2014), dog (Wijesena and Schmutz, 2015) or tiger (Xu et al., 2013). It is
also located downstream RXFP3, a gene known to be involved in food intake regulation and body
weight in mice (Ganella et al., 2012; Smith et al., 2014)

Another interesting region was found on chromosome 22, where two strong candidate polymor-
phisms were located 30kb upstream MAGI1 (Figures S14 and S15) Although not reported in Table 3
due to p-values slightly greater than 10<sup>-4</sup>, 5 other suggestive polymorphisms are located within MAGI1
introns. One of these, at position 36011838, is located within a highly conserved region according to
a multiple alignment of 36 eutherian mammals (www.ensembl.org). At this position, all considered
species carry allele G, which is also the most frequent allele in Angus, Fleckvieh and Holstein, while
allele A swept in Jersey. MAGI1 is a scaffolding protein present in tight junction of epithelial cells
and might be implied in nervous system functions. Adaptive selection around this gene was already
reported in some West African cattle breeds (Gautier et al., 2009).

A more complex situation was observed on chromosome 7, where 8 and 16 candidate polymorphisms were found in two sweep windows distant from 300kb. Among all these candidates, the highest FLK value was reached at position 26459812, in an intergenic region between SLC27A6 and FBN2. Interestingly, this polymorphism was located in a very conserved region, and it was the only one in this case among all candidates in the sweep window. Allele T was found in almost all species at this position and was almost fixed in Angus, Fleckvieh and Holstein, while allele A swept in Jersey. SLC27A6 has been shown related to fatty acid metabolism in cattle (Bionaz and Loor, 2008; Nafikov et al., 2013), while FBN2 has been related to several development processes, including for instance bone formation, in mice (Nistala et al., 2010).

Strong candidate mutations were also found in two regions without protein-coding genes on the
current bovine annotation, and in general we note that all candidate polymorphisms were located in
intergenic or regulatory regions. This implies that validating the effect of these polymorphisms will be
difficult, but this also outlines the potential of selective sweep studies to improve genome annotation.
Hard sweep signals shared by all breeds

In regions where hard sweep signals were shared by all breeds (Table 1), identifying causal poly-
morphisms only from our dataset is impossible. Indeed, the majority of polymorphisms have similar
diversity patterns, with a low minor allele frequency in all 4 breeds. Besides, positions that were
monomorphic in our data set are also good candidates, since they might result from the complete
fixation of a favorable allele in the 4 breeds.

In order to identify potential causal variants in these regions, one possible approach can be to use
data from related species, looking for highly conserved positions for which the major (or the only)
allele of our bovine data set is absent in other mammals. To illustrate this approach, we implemented
the following procedure. Based on a multiple alignment of the bovine reference sequence with 10 other
mammal reference sequences (Rocha et al., 2014), we selected the positions where (i) the bovine allele
(or the major bovine allele in the case of polymorphic positions) was distinct from the yak allele, (ii)
the bovine allele was unobserved among the 10 other species, and (iii) the yak allele was observed
in more than 7 (out of 9) other species, including at least buffalo or sheep (the two closest species).
Such positions represent convincing candidates for two reasons. First, alleles that are observed in
other mammal species must have been present among bovines for a very long time and are thus very
unlikely to produce a hard sweep pattern, even if they became positively selected at some point in
the bovine history. Second, highly conserved positions are more likely functional, and thus subject to
selection, than less conserved ones.

Among the 1,382,681 positions included in the regions of Table 1, only 91 satisfied the 3 conditions
above. Further removing 7 positions for which the minor bovine allele was at quite high frequency, we
finally obtained 84 causal candidates (Table S6). All sweep regions except one exhibited convincing
causal variants located in coding or regulatory regions, but more work will be needed to determine
the exact causal variants in these regions.

Discussion

We performed in this study a genomic scan for selection in European taurine cattle, based on large
samples of sequencing data from four different breeds. We used two detection approaches, based
respectively on the genetic diversity within breeds and on the genetic differentiation between breeds,
and compared the signals detected by these two approaches.

One important conclusion from our analysis is that sequencing data represents a great opportunity
in the context of genomic scans for selection. Indeed, the detection power of hapFLK was higher when
applied to sequencing data than when applied to high-density chip data (Figure 5). This was consistent
with previous studies looking for selection signatures in cattle (Ma et al., 2015) and humans (Liu et al.,
2014), which also found, based on computer simulations, that higher power could be expected from
sequencing data compared to SNP chip data. The localization of selection signatures was also found
more accurate with sequencing data, both for hapFLK and the within population approach, as the
median size of detected regions (a few tens of kilobases) was considerably lower than with SNP chip
data. Note that the reduced size of detected regions does not mean that we detected older selection
442

443

444

445

446

447

448

449

450

451

452

453

454

455

456

457

458

459

460

461

462

463

464

465

466

467

468

469

470

471

472

473

474

475

476

477

478

479

480

481

482

events. It just comes from the fact that, in many regions, only the most significant part of the
sweep was identified. For instance, the selection events detected by hapFLK are quite recent, because
they must have occured more recently than breed divergence (approximately 500 years before present
according to our demographic analysis). Still, the region with excess differentiation, which is captured
by hapFLK, is generally smaller than the sweep itself, because part of the swept haplotype can be
shared between breeds. An example of such a signature can be seen in the MC1R signature (Figure
S16). JE NE SUIS PAS SUR QUE CETTE FIGURE SOIT TRES PARLANTE. EN PLUS DANS
LA TABLE 2 ON VOIT QUE LA TAILLE DU SWEEP DONNEE PAR HAPFLK EST EN FAIT
PLUS GRANDE QUE CELLE DONNE PAR LE HMM. PEUX TU PRECISER OU CHANGER D
EXEMPLE? Importantly, the higher precision provided by the use of NGS data reduces the number
of genes included in each window, allowing potentially an easier exploration of the molecular functions
driving selection.

Since sequencing data capture a large proportion of the SNPs and small indels segregating in a
sample, genomic scans for selection based on such data can be expected to identify even the causal
polymorphism under selection in a given sweep region. We demonstrated that this is indeed the case in
some regions, where the few variants with the highest allele frequency differentiation between breeds,
as measured by the FLK statistic, included the causal variant. For instance, in the MC1R region, the
two mutations that have been shown to affect coat color in the breeds of our study, were included in the
top three FLK values. Similarly, in the PLAG1 region, a good overlap was observed between variants
with top FLK values and the QTNs found by Karim et al. (2011) in a GWAS analysis on stature. In
several other regions, variants with top FLK values also lead to promising candidate polymorphisms
(Tables 3 and S5), which were often located in highly conserved regions or/and in the vicinity of genes
with interesting metabolic functions. However, it is important to point out that identifying the selected
variant from FLK values will essentially be possible in population-specific hard sweep scenarios, where
the favorable allele has almost fixed in the selected population(s), and has remained absent or at
very low frequency in the other sampled populations. In all other situations, and for instance those
where the favorable allele is at intermediate frequency also in neutral populations, the FLK value of
the causal variant will be reached just by chance by many neutral variants, so identifying the causal
variant will require additional data, for instance phenotypes related to the selection constraint.

Our study also illustrated the small overlap between the selection signatures detected by the within
and the between population approach, and allowed to better characterize the regions detected by each
of these approaches. Among the 916 hard sweep regions detected by the within-population approach,
only 8 were detected by hapFLK. One reason is that, in many of the regions where an hard sweep
was detected, the swept allele was likely at quite high frequency also in the other breeds, reducing
the differentiation signal. Indeed, in regions where a sweep was detected in one of the breeds, the
TOR distribution in the other breeds was shifted toward that of swept regions, even if this signal
was not strong enough to be detected by the HMM (Figure S9). This can be explained by the long
shared history between the 4 breeds considered in this study, which have a common geographic origin
and have been strictly isolated only in the last 200 years. Convergent selection of similar alleles in
different breeds may also have occurred, as the same traits were selected in these breeds, but this
is not the most parcimonious hypothesis. Finally, the small proportion of hard sweeps detected
by hapFLK was also likely related to a power issue, because none of the most breed-specific hard sweeps listed in Table S2 was detected by this approach. Although such regions clearly showed some signal of differentiation, this signal was not strong enough to be detected by hapFLK, because the very high levels of genetic drift in cattle breeds implies that only extremely differentiated regions can be considered as significantly under selection. Actually, we can see from Table 2 that in most hard sweep regions detected by hapFLK, there was not one hard sweep in a single population, but at least two sweeps (either complete or incomplete) implying distinct haplotypes in the four breeds, which represents an even stronger differentiation signal (Figure S16). We expect the 1000 bull genomes dataset to grow by inclusion of new populations, that will most likely increase power of differentiation approaches such as FLK and hapFLK by adding information on breeds more closely related to each other than on the initial release.

On the other hand, among the 57 selection signatures detected by hapFLK, 49 were not detected by the within-population approach. This comes from the fact that this approach is specifically designed to detect hard sweep signals, while hapFLK also detects incomplete or soft sweeps (Fariello et al., 2013, 2014). For the 49 regions detected by hapFLK and not by the within-breed approach, the evidence of a hard sweep, as measured by statistic $T_{OR}$, was always very low in all breeds, indicating that the signal detected by hapFLK had nothing to do with a hard sweep in one of the breeds. For example, in the ARL15 signature, the swept haplotype clearly segregates at moderate frequency (~15%) in the Jersey population (Figure S12). This ability of hapFLK to detect incomplete and soft sweeps is extremely interesting for the study of livestock species. Indeed, recent intensive selection in modern livestock breeds has mainly targetted polygenic traits such as milk or meat production, which is much more likely to produce soft and incomplete sweep patterns than hard sweep patterns (Pritchard et al., 2010; Hernandez et al., 2011).

## Conclusion

Our study illustrates how sequencing data offer great advantages for the detection of adaptive loci: higher power and better precision in localizing adaptive genes, and even in some rare cases causal mutations. However, we found that many of our strongest signals and mutations lie in non-coding regions of the genome, hinting that a majority of adaptive mutations are regulatory in nature. A better annotation of genomes, such as obtained by high-throughput post-genomic approaches, in combination with phenotyping in large population samples will be key in uncovering the biological basis of adaptation.

## Materials and Methods

### Samples and sequencing

A total of 234 genome sequences were obtained from the 1000 bull genomes project, Run II (Daetwyler et al., 2014). These included 129 Holsteins (125 Black and 4 Red), 43 Fleckviehs, 47 Angus, and 15 Jersey. We considered all these sequences except the 4 Red Holsteins. We based our analyzes on the
phased and corrected autosomal data produced by the 1000 bull genomes project (Daetwyler et al., 2014). These data included 27,535,425 bi-allelic SNPs and 1,507,728 bi-allelic indels.

## Choice of unrelated animals

To remove potential biases arising from sample size heterogeneity between breeds and inbreeding within breeds, we selected a subset of 25 unrelated animals in Holstein, Fleckvieh and Angus, while keeping all 15 Jersey animals. Within each breed, we computed the genetic relationship matrix (GRM) of all available animals based on SNPs from Chromosome 1 with minor allele frequency (MAF) above 10%, using the GCTA 1.04 software (Yang et al., 2011). We then selected unrelated animals as follows. First, we removed all animals with inbreeding value (the diagonal term of the GRM) above 1.5 (this threshold was chosen because inbreeding values above 1.5 clearly appeared as outliers compared to the rest of the distribution, Figure S1). Second, we considered all animal pairs with genetic relationship above 0.3 in absolute value and removed the most inbred animal of each pair. Third, we performed a hierarchical clustering of the remaining animals based on the distance $d_{ij} = M - GRM_{i,j}$, where $GRM_{i,j}$ is the genetic relationship between animals i and j and M is the maximum value of the GRM, and sampled the 25 most distant animals.

## Estimation of a demographic model for the joint history of the four breeds

In order to model the demographic history of the four breeds under study, we assumed that these breeds diverged simultaneously from a common ancestral population, $T_{DIV}$ generations ago. Although the population tree estimated by FLK (Figure 4) would suggest a slightly more recent divergence between Fleckvieh and Jersey, we considered that this difference was negligible and preferred reducing the number of parameters to be estimated. Based on previous studies suggesting that effective population size in taurine cattle strongly declined since domestication (MacLeod et al., 2013; Boitard et al., 2016), we allowed one population size change from $N_{ANC}$ (the ancestral population size) to $N_{DOM}$ (the “domestic” population size) in the ancestral population, $T_{ANC}$ generations ago. The divergence between breeds implied a second population size change: after this event, each breed i was assumed to have a specific population size $N_i$, which possibly differed from $N_{DOM}$. Finally, we assumed that each breed received, every generation, a proportion m of migrants from any of the other breeds.

We estimated the parameters of this model using the composite likelihood approach implemented in fastsimcoal2 (Excoffier et al., 2013), which is based on the joint Site Frequency Spectrum (SFS). In our dataset, this joint SFS had a very high dimension (51×51×51×31). Consequently, we rather considered the collection of joint SFS obtained for all 6 population pairs (Angus × Fleckvieh, Angus × Holstein ...), following the recommendations of the authors. We computed these observed joint SFS from our data using home-made scripts and provided them as input to fastsimcoal2, version 5.2.8 (may 2015). We performed 50 independent EM estimations and selected that with the highest composite likelihood. For each estimation, we used the folded SFS (option -m) and the default settings -n100000 -N100000 -M0.001 -l10 -L40. In fastsimcoal2, time is scaled in generations. Time in years were obtained by assuming a generation time of 5 years.
Detection of genomic regions with low within-breed diversity

**Model** We looked for hard sweep signatures within each breed using the Hidden Markov Model (HMM) of Boitard et al. (2009). In this model, only bi-allelic variants are considered and the derived allele frequency at variant *i*, denoted $Y_i$, is taken as the observed state of the HMMM at this position. Each variant *i* is also assumed to have a hidden state $X_i$, which can take 3 different values: “Selection”, for variants that are very close to a swept site, “Neutral”, for variants that are far away from any swept site, and “Intermediate” for variants in between. These three values are associated with different allele frequency distributions. The “Neutral” allele frequency distribution is estimated using all variants in the genome, assuming most of them have indeed evolved under neutrality. Allele frequency distributions in states “Intermediate” and “Selection” are deduced from this “Neutral” distribution using the derivations in Nielsen et al. (2005), and are typically more skewed towards extreme allele frequencies. The hidden states $X_i$ form a Markov chain along the genome with a per bp probability *p* of switching state, so that close variants tend to be in the same hidden state. Under this HMM, the most likely sequence of hidden states can be predicted from the sequence of observed states using the Viterbi algorithm. Each set of consecutive variants with predicted state “Selection” is called a sweep window.

**Implementation** Ancestral *Bovinae* alleles at 448,289 SNPs included in the Illumina BovineHD BeadChip were obtained from Utsunomiya et al. (2013). To check this information we also aligned the bovine sequence against the sequence of three *Bovidae* species Rocha et al. (2014). For 365,146 SNPs, the ancestral allele in our alignment was consistent with that reported by Utsunomiya et al. (2013) so we used this information for sweep detection. For all other SNPs, we used a folded allele frequency distribution, i.e allele frequencies $Y_i$ and $1 - Y_i$ were considered as the same observed state. Indels were not included at this stage of the analysis (they were only considered when looking at candidate polymorphisms within the region).

To reduce computation time, we estimated the neutral SFS in each breed using only 5% of the SNPs from each chromosome, which were selected at random. These SFS are shown in Figure S2.

The type I error of the above method, i.e. the probability that it detects a sweep window in a population that has evolved under neutrality, depends on parameter *p* (see Boitard et al. (2009) for more details). To control the genome wide number of false positives, we simulated 5,000 samples of length 500kb under neutral evolution using *ms* (Hudson, 2002) and adjusted *p* so that sweeps were detected in only 0.1% of these samples. We performed this calibration for each breed, using the same sample size and proportion of unfolded sites as in the data used for sweep detection. Parameter $\theta$ was also estimated from this data using Waterson's estimator (Waterson, 1975) and $\rho$ was taken equal to $\theta/10$, which is rather low compared to current estimates in cattle (Sandor et al., 2012; Ma et al., 2015). As the detection sensitivity of the HMM method increases when $\rho/\theta$ decreases (Boitard et al., 2009), our adjusted value of *p* should be conservative. Assuming a 2.5Gb genome as that of bovine (focusing on autosomes) is equivalent to 5,000 windows of 500kb, so we expected no more than 5 false positive signals over the genome with this value of *p*.
**Influence of demography on sweep detection** We evaluated the robustness of the above detection approach under two neutral demographic scenarios: the multi-population model estimated by fastsimcoal2 (Figure 1) and the single-population model proposed in (Boitard et al., 2016). For each model, we simulated 20,000 samples of length 500kb using *ms*, assuming a mutation rate and a recombination rate of 1e-8 per bp and generation, and we applied the HMM procedure to these samples. For the multi-population model, each simulated sample included genomes from the 4 breeds, which were splitted in order to apply the HMM within each population. For the single-population scenario, breed-specific samples were directly simulated independently of each others, each breed being associated to a different population size history. For each scenario and breed, the total length of simulated samples was equivalent to that of 4 cattle genomes. Consequently, the estimated proportion *m* of false positive signals per genome was given by the total number of sweeps detected in the simulations, divided by four, and the variance of this estimation was equal to *m*/4. The confidence interval of this estimation was approximated by

$$[m - 2\sigma(m); m + 2\sigma(m)] = [m - \sqrt{m}; m + \sqrt{m}]$$

**Comparing Hard-sweep signatures detected in different populations** For a given variant *i*, the evidence for a hard sweep window occurring around this variant in population *j* was measured by the statistic $T_{OR}(i, j) = \log_{10}(q_i^j/(1 - q_i^j))$, where $q_i^j$ was the posterior probability of hidden state “Selection” returned by the backward-forward algorithm applied to the HMM in population *j*. When considering a given region of the genome, the evidence for a hard sweep in population *j* was quantified by the median of $T_{OR}(i, j)$ over the variants of the region. In order to detect selection signatures that are really breed-specific, we computed for each breed the distribution of $T_{OR}$ in 3 different classes of regions: (i) those where a hard sweep was detected in this breed, (ii) those where a hard sweep was detected in another breed, and (iii) those where no hard sweep was detected (figure S9). Obviously, class (i) and class (iii) regions lead to very different distributions, with lower values $T_{OR}$ in the latter (i.e. in the completely neutral regions). In addition, the distribution of $T_{OR}$ in class (ii) regions was not similar to that found in class (iii) regions, but was shifted towards that of class (i) regions. Based on these observations, we therefore considered that a hard sweep was breed-specific when, in all other breeds, the value of $T_{OR}$ of the region was below a given quantile *q* of the class (iii) distribution. For *q* = 0.5, 55 breed-specific sweeps were detected (listed in File S3) and for *q* = 0.25, 12 breed-specific regions were detected.

**Detection of genomic regions with large differentiation between breeds**

We applied two methods for the detection of genomic regions exhibiting large genetic differentiation between populations: FLK (Bonhomme et al., 2010), a single marker approach based on allele frequency differences, and its haplotypic extension hapFLK (Fariello et al., 2013), which exploits the Linkage Disequilibrium information to capture differences of haplotype frequencies.

**Kinship matrix** In contrast with the $F_{ST}$ statistic, FLK and hapFLK account for the population history through a kinship matrix, which captures (i) differences in effective population sizes between
populations and (ii) possible shared ancestry between populations. The kinship matrix is inferred
from a population tree, with branch length expressed in units of drift, i.e. measured in fixation indices
$f \approx t/2N$ where $t$ is the number of generations from the root and $N$ the effective population size.
We estimated the population tree using neighbour-joining on the Reynold's genetic distances between
populations (see Bonhomme et al. (2010) for details). We used the ancestral allele reconstruction of
Utsunomiya et al. (2013) to root the population tree and estimate the population kinship matrix.

**FLK** For the single marker analysis, we performed the FLK test on all variants and computed p-values using the theoretical $\chi^2(3)$ distribution, which was a good fit to the observed distribution (Figure S3).

**hapFLK** For the hapFLK test, to save computation time, we removed variants that had low minor allele frequency (< 10%) when merging animals from all breeds. Note that as the analysis looks for signal of differentiation, removing these variants does not preclude detection. In subsequent reanalysis of small genomic regions, we kept all variants in the analysis. hapFLK makes use of the local clustering approach in Scheet and Stephens (2006) to model haplotype diversity. This model requires specifying a number of haplotype clusters as input parameter. Using the cross validation procedure implemented in the fastPHASE software, we found that 15 clusters provided the smallest imputation error rate. hapFLK can be computed on unphased or phased genotype data. Genotype calling made use of imputation approaches, and data was therefore already phased (Daetwyler et al., 2014). We computed hapFLK both on the haplotype data and on the genotype (but imputed) data and found the two analyzes provided similar results (not shown).

The distribution of hapFLK is not known, but, from theoretical arguments hapFLK is a deviance statistic. However, the variance parameter of this deviance is not known. If this parameter was known then hapFLK should follow a $\chi^2((N-1)(K-1))$ distribution, where N is the number of populations and K the number of haplotype clusters. Building on this fact, we compared a set of quantiles (from 0.05 to 0.95 every 0.05) of the $\chi^2(42)$ distribution $t_q$ to the observed quantiles of the hapFLK statistic $o_q$. We found that the relationship between $t_q$ and $o_q$ was very close to linear (Figure S4). Thus, we used the parameter of the linear model to scale the hapFLK statistic to a $\chi^2(42)$ distribution which was then used to compute p-values.

**Extracting significantly differentiated regions** In order to call significant regions, we applied the approach of Storey and Tibshirani (2003) aimed at controlling the False Discovery Rate (FDR), at the 15% level, i.e. we called significant variants with q-values < 0.15 and selection signatures regions where more than one variant was called significant. We note however that our level of control of the FDR for the regions themselves may not be well calibrated as the tests are correlated. This procedure still provides a significant threshold so that among marker discoveries there is a six fold enrichment in favor of the number of statistics under the alternative.
652

References

653

654

Albrecht E, Komolka K, Kuzinski J, Maak S, 2012. Agouti revisited: Transcript quantification of the asip gene in bovine tissues related to protein expression and localization. *PLoS ONE* 7:e35282.

655

656

657

Allais-Bonnet A, Grohs C, Medugorac I, Krebs S, Djari A, Graf A, Fritz S, Seichter D, Baur A, Russ I, et al., 2013. Novel insights into the bovine polled phenotype and horn ontogenesis in <italic>bovidae</italic>. *PLoS ONE* 8:e63512.

658

659

660

Bionaz M, Loor JJ, 2008. Acs11, agpat6, fabp3, lpin1, and slc27a6 are the most abundant isoforms in bovine mammary tissue and their expression is affected by stage of lactation. *The Journal of nutrition* 138:1019-1024.

661

Biswas S, Akey JM, 2006. Genomic insights into positive selection. *Trends in Genetics* 22:437-446.

662

663

664

Boitard S, Rodriguez W, Jay F, Mona S, Austerlitz F, 2016. Inferring population size history from large samples of genome wide molecular data - an approximate bayesian computation approach. *Bioarxiv*.

665

666

Boitard S, Schlotterer C, Futschik A, 2009. Detecting selective sweeps: a new approach based on hidden markov models. *Genetics* 181:1567-1578.

667

668

669

Bole-Feysot C, Goffin V, Edery M, Binart N, Kelly PA, 1998. Prolactin (prl) and its receptor: actions, signal transduction pathways and phenotypes observed in prl receptor knockout mice. *Endocrine Reviews* 19:225-268. PMID: 9626554.

670

671

Bongiorni S, Mancini G, Chillemi G, Pariset L, Valentini A, 2012. Identification of a short region on chromosome 6 affecting direct calving ease in Piedmontese cattle breed. *PLoS One* 7:e50137.

672

673

Bonhomme M, Chevalet C, Servin B, Boitard S, Abdallah J, Blott S, SanKristobal M, 2010. Detecting selection in population trees: the lewontin and krakauer test extended. *Genetics* 186:241-262.

674

675

676

Daetwyler HD, Capitan A, Pausch H, Stothard P, van Binsbergen R, Brondum RF, Liao X, Djari A, Rodriguez SC, Grohs C, et al., 2014. Whole-genome sequencing of 234 bulls facilitates mapping of monogenic and complex traits in cattle. *Nat Genet* 46:858-865.

677

678

679

Di Noia MA, Todisco S, Cirigliano A, Rinaldi T, Agrimi G, Iacobazzi V, Palmieri F, 2014. The human slc25a33 and slc25a36 genes of solute carrier family 25 encode two mitochondrial pyrimidine nucleotide transporters. *Journal of Biological Chemistry* 289:33137-33148.

680

681

Dickinson RE, Duncan WC, 2010. The SLIT-ROBO pathway: a regulator of cell function with implications for the reproductive system. *Reproduction* 139:697-704.

682

683

684

Dickinson RE, Hryhorskyj L, Tremewan H, Hogg K, Thomson AA, McNeilly AS, Duncan WC, 2010. Involvement of the SLIT/ROBO pathway in follicle development in the fetal ovary. *Reproduction* 139:395-407.
Eberlein A, Takasuga A, Setoguchi K, Pfuhl R, Flisikowski K, Fries R, Klopp N, Furbass R, Weikard R, Kuhn C, 2009. Dissection of genetic factors modulating fetal growth in cattle indicates a substantial role of the non-SMC condensin I complex, subunit G (NCAPG) gene. *Genetics* 183:951–64.

Echternkamp S, Aad P, Eborn D, Spicer L, 2012. Increased abundance of aromatase and follicle stimulating hormone receptor mrna and decreased insulin-like growth factor-2 receptor mrna in small ovarian follicles of cattle selected for twin births. *Journal of animal science* 90:2193–2200.

Excoffier L, Dupanloup I, Huerta-SÁnchez E, Sousa VC, Foll M, 2013. Robust demographic inference from genomic and snp data. *PLoS Genet* 9:e1003905.

Fariello MI, Boitard S, Naya H, SanCristobal M, Servin B, 2013. Detecting signatures of selection through haplotype differentiation among hierarchically structured populations. *Genetics* 193:929–941.

Fariello MI, Servin B, Tosser-Klopp G, Rupp R, Moreno C, Cristobal MS, Boitard S, Consortium ISG, 2014. Selection signatures in worldwide sheep populations. *PLoS ONE* 9:e103813.

Felius M, Koolmees PA, Theunissen B, Consortium ECGD, Lenstra JA, 2011. On the breeds of cattle—historic and current classifications. *Diversity* 3:660–692.

Flori L, Fritz S, Jaffrézic F, Boussaha M, Gut I, Heath S, Foulley JL, Gautier M, 2009. The genome response to artificial selection: A case study in dairy cattle. *PLoS ONE* 4:e6595.

Ganella DE, Ryan PJ, Bathgate RA, Gundlach AL, 2012. Increased feeding and body weight gain in rats after acute and chronic activation of rxfp3 by relaxin-3 and receptor-selective peptides: functional and therapeutic implications. *Behavioural pharmacology* 23:516–525.

Gautier M, Flori L, Riebler A, Jaffrezic F, Laloe D, Gut I, Moazami-Goudarzi K, Foulley JL, 2009. A whole genome bayesian scan for adaptive genetic divergence in west african cattle. *BMC Genomics* 10:550.

Giesy SL, Yoon B, Currie WB, Kim JW, Boisclair YR, 2012. Adiponectin deficit during the precarious glucose economy of early lactation in dairy cows. *Endocrinology* 153:5834–44.

Gouveia JJdS, Silva MVGBd, Paiva SR, Oliveira SAMPd, 2014. Identification of selection signatures in livestock species. *Genetics and Molecular Biology* 37:330–342.

Gutierrez-Gil B, Arranz JJ, Wiener P, 2015. An interpretive review of selective sweep studies in bos taurus cattle populations: identification of unique and shared selection signals across breeds. *Frontiers in Genetics* 6.

Hayes BJ, Pryce J, Chamberlain AJ, Bowman PJ, Goddard ME, 2010. Genetic architecture of complex traits and accuracy of genomic prediction: coat colour, milk-fat percentage, and type in holstein cattle as contrasting model traits. *PLoS Genetics* 6:e1001139.

Hernandez RD, Kelley JL, Elyashiv E, Melton SC, Auton A, McVean G, Sella G, Przeworski M, et al., 2011. Classic selective sweeps were rare in recent human evolution. *science* 331:920–924.
720

Hudson R, 2002. Generating samples under a Wright-Fisher neutral model of genetic variation. *Bioinformatics* 18:337-338.

721

722

Karim L, Takeda H, Lin L, Druet T, Arias JA, Baurain D, Cambisano N, Davis SR, Farnir F, Grisart B, et al., 2011. Variants modulating the expression of a chromosome domain encompassing plagl influence bovine stature. *Nature genetics* 43:405-413.

723

724

725

Karisa B, Thomson J, Wang Z, Stothard P, Moore S, Plastow G, 2013. Candidate genes and single nucleotide polymorphisms associated with variation in residual feed intake in beef cattle. *Journal of animal science* 91:3502-3513.

726

727

728

Kijas JW, 2014. Haplotype-based analysis of selective sweeps in sheep. *Genome* 57:433-7.

729

Killeen AP, Morris DG, Kenny DA, Mullen MP, Diskin MG, Waters SM, 2014. Global gene expression in endometrium of high and low fertility heifers during the mid-luteal phase of the estrous cycle. *BMC genomics* 15:234.

730

731

732

Klungland H, Vage D, Gomez-Raya L, Adalsteinsson S, Lien S, 1995. The role of melanocyte-stimulating hormone (msh) receptor in bovine coat color determination. *Mammalian genome* 6:636-639.

733

734

735

Leroy G, Mary-Huard T, Verrier E, Danvy S, Charvolin E, Danchin-Burge C, 2013. Methods to estimate effective population size using pedigree data: Examples in dog, sheep, cattle and horse. *Genetics Selection Evolution* 45:1.

736

737

738

Li M, Tian S, Jin L, Zhou G, Li Y, Zhang Y, Wang T, Yeung CK, Chen L, Ma J, et al., 2013. Genomic analyses identify distinct patterns of selection in domesticated pigs and tibetan wild boars. *Nature genetics* 45:1431-1438.

739

740

741

Lindholm-Perry AK, Kuehn LA, Oliver WT, Sexten AK, Miles JR, Rempel LA, Cushman RA, Freetly HC, 2013. Adipose and muscle tissue gene expression of two genes (NCAPG and LCORL) located in a chromosomal region associated with cattle feed intake and gain. *PLoS One* e80882.

742

743

744

Lindholm-Perry AK, Sexten AK, Kuehn LA, Smith TP, King DA, Shackelford SD, Wheeler TL, Ferrell CL, Jenkins TG, Snelling WM, et al., 2011. Association, effects and validation of polymorphisms within the NCAPG - LCORL locus located on BTA6 with feed intake, gain, meat and carcass traits in beef cattle. *BMC Genet* 12:103.

745

746

747

748

Liu X, Saw WY, Ali M, Ong RTH, Teo YY, 2014. Evaluating the possibility of detecting evidence of positive selection across asia with sparse genotype data from the hugo pan-asian snp consortium. *BMC genomics* 15:332.

749

750

751

Liu YJ, Liu XG, Wang L, Dina C, Yan H, Liu JF, Levy S, Papasian CJ, Drees BM, Hamilton JJ, et al., 2008. Genome-wide association scans identified ctnnbl1 as a novel gene for obesity. *Human Molecular Genetics* 17:1803-1813.

752

753
754

755

756

Ma L, O'Connell JR, VanRaden PM, Shen B, Padhi A, Sun C, Bickhart DM, Cole JB, Null DJ, Liu GE, et al., 2015. Cattle sex-specific recombination and genetic control from a large pedigree analysis. *PLoS Genet* 11:e1005387.

757

758

759

MacLeod IM, Larkin DM, Lewin HA, Hayes BJ, Goddard ME, 2013. Inferring demography from runs of homozygosity in whole-genome sequence, with correction for sequence errors. *Molecular Biology and Evolution* 30:2209–2223.

760

761

762

Makvandi-Nejad S, Hoffman GE, Allen JJ, Chu E, Gu E, Chandler AM, Loredo AI, Bellone RR, Mezey JG, Brooks SA, et al., 2012. Four loci explain 83% of size variation in the horse. *PLoS One* 7:e39929.

763

764

Metzger J, Schrimpf R, Philipp U, Distl O, 2013. Expression levels of LCORL are associated with body size in horses. *PLoS One* 8:e56497.

765

766

767

768

Morice-Picard F, Lasseaux E, Cailley D, Gros A, Toutain J, Plaisant C, Simon D, François S, Gilbert-Dussardier B, Kaplan J, et al., 2014. High-resolution array-cgh in patients with oculocutaneous albinism identifies new deletions of the tyr, oca2, and slc45a2 genes and a complex rearrangement of the oca2 gene. *Pigment cell & melanoma research* 27:59–71.

769

770

771

Morris AP, Voight BF, Teslovich TM, Ferreira T, Segre AV, Steinthorsdottir V, Strawbridge RJ, Khan H, Grallert H, Mahajan A, et al., 2012. Large-scale association analysis provides insights into the genetic architecture and pathophysiology of type 2 diabetes. *Nat Genet* 44:981–90.

772

773

774

775

Nafikov R, Schoonmaker J, Korn K, Noack K, Garrick D, Koehler K, Minick-Bormann J, Reecy J, Spurlock D, Beitz D, 2013. Association of polymorphisms in solute carrier family 27, isoform a6 (slc27a6) and fatty acid-binding protein-3 and fatty acid-binding protein-4 (fabp3 and fabp4) with fatty acid composition of bovine milk. *Journal of dairy science* 96:6007–6021.

776

777

Nielsen R, Williamson L, Kim Y, Hubisz M, Clark A, Bustamante C, 2005. Genomic scans for selective sweeps using SNP data. *Genome Research* 15:1566–1575.

778

779

780

Nistala H, Lee-Arteaga S, Smaldone S, Siciliano G, Carta L, Ono RN, Sengle G, Arteaga-Solis E, Levasseur R, Ducy P, et al., 2010. Fibrillin-1 and-2 differentially modulate endogenous tgf-β and bmp bioavailability during bone formation. *The Journal of cell biology* 190:1107–1121.

781

782

783

Pausch H, Flisikowski K, Jung S, Emmerling R, Edel C, Götz KU, Fries R, 2011. Genome-wide association study identifies two major loci affecting calving ease and growth-related traits in cattle. *Genetics* 187:289–297.

784

785

Pritchard JK, Pickrell JK, Coop G, 2010. The genetics of human adaptation: hard sweeps, soft sweeps, and polygenic adaptation. *Current Biology* 20:R208–R215.

786

787

788

Qanbari S, Gianola D, Hayes B, Schenkel F, Miller S, Moore S, Thaller G, Simianer H, 2011. Application of site and haplotype-frequency based approaches for detecting selection signatures in cattle. *BMC Genomics* 12:318.
789

Qanbari S, Pausch H, Jansen S, Somel M, Strom TM, Fries R, Nielsen R, Simianer H, 2014. Classic selective sweeps revealed by massive sequencing in cattle. *PLoS Genet* 10:e1004148.

790

791

Randhawa IA, Khatkar MS, Thomson PC, Raadsma HW, 2015. Composite Selection Signals for Complex Traits Exemplified Through Bovine Stature Using Multi-Breed Cohorts of European and African Bos taurus. *G3 (Bethesda)*.

792

793

794

Richards JB, Waterworth D, O'Rahilly S, Hivert MF, Loos RJ, Perry JR, Tanaka T, Timpson NJ, Semple RK, Soranzo N, et al., 2009. A genome-wide association study reveals variants in ARL15 that influence adiponectin levels. *PLoS Genet* 5:e1000768.

795

796

797

Rivadeneira F, Styrkarsdottir U, Estrada K, Halldorsson BV, Hsu YH, Richards JB, Zillikens MC, Kavvoura FK, Amin N, Aulchenko YS, et al., 2009. Twenty bone-mineral-density loci identified by large-scale meta-analysis of genome-wide association studies. *Nat Genet* 41:1199-206.

798

799

800

Rocha D, Billerey C, Samson F, Boichard D, Boussaha M, 2014. Identification of the putative ancestral allele of bovine single-nucleotide polymorphisms. *Journal of Animal Breeding and Genetics* 131:483-486.

801

802

803

Roux PF, Boitard S, Blum Y, Parks B, Montagner A, Mouisel E, Djari A, Esquerré D, Désert C, Boutin M, et al., 2015. Combined qtl and selective sweep mappings with coding snp annotation and cis-eqtl analysis revealed park2 and jag2 as new candidate genes for adiposity regulation. *G3: Genes/Genomes/Genetics* 5:517-529.

804

805

806

807

Rubin CJ, Megens HJ, Barrio AM, Maqbool K, Sayyab S, Schwochow D, Wang C, Carlborg, Jern P, Jørgensen CB, et al., 2012. Strong signatures of selection in the domestic pig genome. *Proceedings of the National Academy of Sciences* 109:19529-19536.

808

809

810

Rubin CJ, Zody MC, Eriksson J, Meadows JR, Sherwood E, Webster MT, Jiang L, Ingman M, Sharpe T, Ka S, et al., 2010. Whole-genome resequencing reveals loci under selection during chicken domestication. *Nature* 464:587-591.

811

812

813

Sabeti P, Schaffner S, Fry B, Lohmueller J, Varilly P, Shamovsky O, Palma A, Mikkelsen T, Altshuler D, Lander E, 2006. Positive natural selection in the human lineage. *science* 312:1614-1620.

814

815

Sahana G, Hoglund JK, Guldbrandtsen B, Lund MS, 2015. Loci associated with adult stature also affect calf birth survival in cattle. *BMC Genet* 16:47.

816

817

Sandor C, Li W, Coppieters W, Druet T, Charlier C, Georges M, 2012. Genetic variants in <italic>rec8</italic>, <italic>rnf212</italic>, and <italic>prdm9</italic> influence male recombination in cattle. *PLoS Genet* 8:e1002854.

818

819

820

Scheet P, Stephens M, 2006. A fast and flexible statistical model for large-scale population genotype data: applications to inferring missing genotypes and haplotypic phase. *The American Journal of Human Genetics* 78:629-644.

821

822
Signer-Hasler H, Flury C, Haase B, Burger D, Simianer H, Leeb T, Rieder S, 2012. A genome-wide association study reveals loci influencing height and other conformation traits in horses. *PLoS One* 7:e37282.

Smith CM, Chua BE, Zhang C, Walker AW, Haidar M, Hawkes D, Shabanpoor F, Hossain MA, Wade JD, Rosengren KJ, et al., 2014. Central injection of relaxin-3 receptor (rxfp3) antagonist peptides reduces motivated food seeking and consumption in c57bl/6j mice. *Behavioural brain research* 268:117-126.

Song H, Kim H, Lee K, Lee DH, Kim TS, Song JY, Lee D, Choi D, Ko CY, Kim HS, et al., 2012. Ablation of Rassf2 induces bone defects and subsequent haematopoietic anomalies in mice. *EMBO J* 31:1147-59.

Soung do Y, Dong Y, Wang Y, Zuscik MJ, Schwarz EM, O'Keefe RJ, Drissi H, 2007. Runx3/AML2/Cbfa3 regulates early and late chondrocyte differentiation. *J Bone Miner Res* 22:1260-70.

Stahl EA, Raychaudhuri S, Remmers EF, Xie G, Eyre S, Thomson BP, Li Y, Kurreeman FA, Zhernakova A, Hinks A, et al., 2010. Genome-wide association study meta-analysis identifies seven new rheumatoid arthritis risk loci. *Nat Genet* 42:508-14.

Stefanaki I, Panagiotou OA, Kodela E, Gogas H, Kypreou KP, Chatzinasiou F, Nikolaou V, Plaka M, Kalfa I, Antoniou C, et al., 2013. Replication and predictive value of snps associated with melanoma and pigmentation traits in a southern european case-control study. *PloS one* 8:e55712.

Storey JD, Tibshirani R, 2003. Statistical significance for genomewide studies. *Proceedings of the National Academy of Sciences* 100:9440-9445.

Sturm RA, 2009. Molecular genetics of human pigmentation diversity. *Human molecular genetics* 18:R9-R17.

Tetens J, Widmann P, Kuhn C, Thaller G, 2013. A genome-wide association study indicates LCORL/NCAPG as a candidate locus for withers height in German Warmblood horses. *Anim Genet* 44:467-71.

The Bovine HapMap Consortium, 2009. Genome-wide survey of snp variation uncovers the genetic structure of cattle breeds. *Science* 324:528-532.

Tonevitsky AG, Maltseva DV, Abbasi A, Samatov TR, Sakharov DA, Shkurnikov MU, Lebedev AE, Galatenko VV, Grigoriev AI, Northoff H, 2013. Dynamically regulated mirna-mrna networks revealed by exercise. *BMC physiology* 13:9.

Utsunomiya YT, Perez OBrien AM, Sonstegard TS, Van Tassell CP, do Carmo AS, MÃ©szÃ¡ros GA, SÃlkner J, Garcia JAF, 2013. Detecting loci under recent positive selection in dairy and beef cattle by combining different genome-wide scan methods. *PLoS ONE* 8:e64280.
857

Viitala S, Szyda J, Blott S, Schulman N, Lidauer M, MĀki-Tanila A, Georges M, Vilkki J, 2006. The role of the bovine growth hormone receptor and prolactin receptor genes in milk, fat and protein production in finnish ayrshire dairy cattle. *Genetics* 173:2151–2164.

858

859

860

Voight BF, Kudaravalli S, Wen X, Pritchard JK, 2006. A map of recent positive selection in the human genome. *PLoS Biol* 4:e72.

861

862

Watterson G, 1975. On the number of segregating sites in genetical models without recombination. *Theor Popul Biol* 7:256–276.

863

864

Weikard R, Altmaier E, Suhre K, Weinberger KM, Hammon HM, Albrecht E, Setoguchi K, Takasuga A, Kuhn C, 2010. Metabolomic profiles indicate distinct physiological pathways affected by two loci with major divergent effect on Bos taurus growth and lipid deposition. *Physiol Genomics* 42A:79–88.

865

866

867

868

Wijesena HR, Schmutz SM, 2015. A missense mutation in slc45a2 is associated with albinism in several small long-haired dog breeds. *Journal of Heredity* p. esv008.

869

870

Xu L, Bickhart DM, Cole JB, Schroeder SG, Song J, Tassell CP, Sonstegard TS, Liu GE, 2015. Genomic signatures reveal new evidences for selection of important traits in domestic cattle. *Mol Biol Evol* 32:711–25.

871

872

873

Xu X, Dong GX, Hu XS, Miao L, Zhang XL, Zhang DL, Yang HD, Zhang TY, Zou ZT, Zhang TT, et al., 2013. The genetic basis of white tigers. *Current Biology* 23:1031–1035.

874

875

Yang J, Lee SH, Goddard ME, Visscher PM, 2011. Gcta: A tool for genome-wide complex trait analysis. *The American Journal of Human Genetics* 88:76 – 82.

876

877

Yoshida CA, Yamamoto H, Fujita T, Furuichi T, Ito K, Inoue K, Yamana K, Zanma A, Takada K, Ito Y, et al., 2004. Runx2 and Runx3 are essential for chondrocyte maturation, and Runx2 regulates limb growth through induction of Indian hedgehog. *Genes Dev* 18:952–63.

878

879
<figure><img src="image_3.png" /><figcaption>Figure 1: **Joint demography of the four cattle breeds**, estimated from our data using the approach of Excoffier et al. (2013). Population sizes correspond to the number of haploid individuals. Parameter values correspond to the EM iteration with the highest likelihood, but similar values were obtained from the second and third best EM iterations.</figcaption></figure>
overlap with Fleckvieh sweeps

<figure><img src="image_4.png" /></figure>

overlap with Holstein sweeps

<figure><img src="image_5.png" /><figcaption>**Figure 2: Comparison between HMM and CLR results.** Proportion of Fleckvieh CLR p-values (obtained from Qanbari et al. (2014)) in sweep windows vs other windows in the genome. On the left panel, the sweep windows considered were those detected in Fleckvieh, independently of what happened in other breeds. On the right panel, the sweep windows considered were those detected only in Holstein.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Figure 3: Allele frequencies in the sweep region at the polled locus. For SNPs where the ancestral allele is known (in red), the frequency is that of the derived allele. For other SNPs (in black) the frequency is that of the minor allele (among all breeds). Vertical red bars delimit the union of detected regions among the four breeds.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Figure 4: Population tree estimated on the 1000 bull genomes data. Branch length is measured in units of drift (≈ t/2N, where t is the time in generations and N the effective population size).</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Figure 5: **Influence of NGS on hapFLK detection power.** PP plot of the hapFLK test applied to SNPs of the Illumina BovineHD SNP chip (left) or to all sites of the 1000 bull genomes project (right)</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Figure 6: **Distribution of hapFLK within hard sweep regions.** Comparison of the distribution of hapFLK p-values in hard sweep regions identified within populations and in the rest of the genome. The plot shows the ratio of the p-values distribution in hard sweeps detected in 1 to 4 populations to their distribution on the part of the genome where no hard sweep is detected. y-axis is on a log2 scale and x-axis on a log10 scale.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Figure 7: **Selection signature around PLAG1.** HapFLK (top) and FLK (middle) p-values (log10 scale) for the selection signature around PLAG1, and local heterozygosity in the four breeds (bottom) for the same region. On the top and middle panels, genes are indicated by purple full rectangles, and red full triangles correspond to the candidate QTNs of Karim et al. (2011).</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Figure 8: **Selection signature around MC1R.** HapFLK (gray) and FLK (black) p-values (log10 scale) for the selection signature around MC1R (purple vertical line). Red triangle highlight the two known causal mutations for red and black coat color in cattle.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>Figure 9: Allele frequencies in the sweep region including SLC45A2 and RXFP3. For SNPs where the ancestral allele is known (in red), the frequency is that of the derived allele. For other SNPs (in black) the frequency is that of the minor allele (among all breeds).</figcaption></figure>
<table><caption>Table 1: Sweep regions shared among all breeds.</caption><thead><tr><th>chr</th><th>start (Mb)</th><th>end (Mb)</th><th>Genes</th></tr></thead><tbody><tr><td>1</td><td>1.781</td><td>1.818</td><td>lincRNA, Polled locus (Allais-Bonnet et al., 2013)</td></tr><tr><td>1</td><td>107.452</td><td>107.557</td><td>PPM1L</td></tr><tr><td>1</td><td>107.571</td><td>107.749</td><td>ARL14</td></tr><tr><td>5</td><td>68.675</td><td>68.751</td><td>SLC41A2</td></tr><tr><td>7</td><td>4.574</td><td>4.745</td><td>FKBP8, ELL, ISYNA1, SSBP4, LRRC25, GDF15</td></tr><tr><td>10</td><td>59.148</td><td>59.338</td><td>CYP19A1</td></tr><tr><td>16</td><td>44.672</td><td>44.956</td><td>CLSTN1, PIK3CD, TMEM201, SLC25A33</td></tr><tr><td>16</td><td>45.644</td><td>45.903</td><td>RERE, SLC45A1</td></tr></tbody></table>
<table><caption>Table 2: **Hard sweep selection signatures associated to significant differentiation signals.** Signatures are ordered by decreasing hapFLK p-values. Region coordinates are expressed in Megabases on assembly UMD 3.1. Population abbreviation: **Angus**, **Fleckvieh**, **Holstein**, **Jersey**</caption><thead><tr><th rowspan="2">BTA</th><th colspan="3">Hard sweep region</th><th colspan="3">hapFLK</th><th rowspan="2">Candidate genes</th></tr><tr><th>start</th><th>end</th><th>population</th><th>start</th><th>end</th><th>p-value</th></tr></thead><tbody><tr><td>6</td><td>71.439</td><td>71.558</td><td>F</td><td>70.332</td><td>71.607</td><td>1.2 10<sup>-12</sup></td><td>KIT</td></tr><tr><td>18</td><td>14.755</td><td>14.963</td><td>F</td><td>14.305</td><td>14.872</td><td>2.1 10<sup>-8</sup></td><td>MC1R</td></tr><tr><td>14</td><td>24.805</td><td>25.076</td><td>H,A</td><td>24.937</td><td>25.070</td><td>2.2 10<sup>-7</sup></td><td>PLAG1</td></tr><tr><td>10</td><td>5.736</td><td>5.843</td><td>A</td><td>5.736</td><td>5.782</td><td>4.9 10<sup>-6</sup></td><td>intergenic</td></tr><tr><td>13</td><td>64.149</td><td>64.197</td><td>A</td><td>63.879</td><td>64.546</td><td>2.0 10<sup>-5</sup></td><td>ASIP</td></tr><tr><td>20</td><td>22.923</td><td>23.203</td><td>J</td><td>23.110</td><td>23.125</td><td>4.2 10<sup>-5</sup></td><td>ANKRD55</td></tr><tr><td>7</td><td>43.436</td><td>43.542</td><td>A,J</td><td>43.473</td><td>43.474</td><td>7.0 10<sup>-5</sup></td><td>OR cluster</td></tr><tr><td>24</td><td>14.006</td><td>14.040</td><td>F</td><td>14.021</td><td>14.023</td><td>8.2 10<sup>-5</sup></td><td>intergenic</td></tr></tbody></table>
Uncovering adaptation from sequence data: lessons
from genome resequencing of four cattle breeds

Simon Boitard¹,², Mekki Boussaha¹, Aurélien Capitan¹,³, Dominique Rocha¹, and Bertrand Servin⁴

¹INRA, UMR1313 GABI (INRA & AgroParisTech), Jouy-en-Josas, 78352, France

²Ecole Pratique des Hautes Etudes, UMR 7205 ISYEB (Muséum National d'Histoire Naturelle & CNRS & EPHE & Université Pierre et Marie Curie), Paris, 75005, France

³Union Nationale des Coopératives d'Elevage et d'Insémination Animale, Paris, 75595, France.

⁴INRA, UMR1388 GenPhySE (INRA & Institut National Polytechnique de Toulouse), Castanet-Tolosan, 31326, France
<table><caption>Table S1: hapFLK significant regions corresponding to likely assembly errors. These regions correspond exactly to contigs or scaffolds that show conserved synteny with human sex chromosomes, and usually with the human Y chromosome, and their borders correspond to gaps in the genome assembly. Sex chromosomes exhibit a different evolutionary history than autosomes, in particular they have a reduced effective population size. Thus, if sex chromosomes sequences are analyzed as autosomes, they will stand out as outliers. Here, we believe these regions correspond to sex chromosome sequences wrongly assigned to autosomes and therefore they were not considered in the rest of our analysis.</caption><thead><tr><th>Chromosome</th><th>start</th><th>end</th><th>hapFLK p-value</th></tr></thead><tbody><tr><td>1</td><td>83,218,059</td><td>83,237,473</td><td>6.610<sup>-7</sup></td></tr><tr><td>5</td><td>22,519,930</td><td>22,585,490</td><td>6.010<sup>-8</sup></td></tr><tr><td>7</td><td>44,194,255</td><td>44,428,089</td><td>1.310<sup>-12</sup></td></tr><tr><td>7</td><td>79,127,958</td><td>79,269,166</td><td>1.410<sup>-14</sup></td></tr><tr><td>13</td><td>5,373,476</td><td>5,575,795</td><td>3.810<sup>-16</sup></td></tr><tr><td>17</td><td>35,794,432</td><td>36,050,136</td><td>2.410<sup>-10</sup></td></tr><tr><td>17</td><td>50,710,246</td><td>50,780,102</td><td>1.710<sup>-6</sup></td></tr><tr><td>17</td><td>51,861,569</td><td>52,083,308</td><td>1.110<sup>-7</sup></td></tr><tr><td>17</td><td>72,519,569</td><td>72,522,168</td><td>3.410<sup>-5</sup></td></tr><tr><td>22</td><td>39,728,012</td><td>40,072,931</td><td>2.010<sup>-11</sup></td></tr></tbody></table>
<table><caption>Table S2: Hard sweep regions specific to one breed.</caption><thead><tr><th>Population</th><th>chr</th><th>start</th><th>end</th><th>genes</th><th>putative function or trait</th></tr></thead><tbody><tr><td>Holstein</td><td>2</td><td>100.51</td><td>100.6</td><td>ERBB4</td><td>neural crest development</td></tr><tr><td>Jersey</td><td>4</td><td>18.06</td><td>18.12</td><td>intergenic</td><td></td></tr><tr><td>Fleckvieh</td><td>10</td><td>83.12</td><td>83.2</td><td>RPL7</td><td>post-weaning gain in sheep [32]</td></tr><tr><td>Angus</td><td>16</td><td>27.64</td><td>27.71</td><td>CAPN8</td><td>gastric mucosal defense [8]</td></tr><tr><td>Angus</td><td>19</td><td>2.51</td><td>2.54</td><td>intergenic</td><td></td></tr><tr><td>Jersey</td><td>20</td><td>0.15</td><td>0.22</td><td>PANK3<br>SPZ1</td><td>bone formation [34]<br>spermatogenesis [11]</td></tr><tr><td>Jersey</td><td>20</td><td>29.67</td><td>29.73</td><td>intergenic</td><td></td></tr><tr><td>Angus</td><td>20</td><td>66.19</td><td>66.22</td><td>intergenic</td><td></td></tr><tr><td>Angus</td><td>21</td><td>49.09</td><td>49.14</td><td>intergenic</td><td></td></tr><tr><td>Angus</td><td>22</td><td>5.73</td><td>5.77</td><td>intergenic</td><td></td></tr><tr><td>Jersey</td><td>27</td><td>3.54</td><td>3.64</td><td>intergenic</td><td></td></tr></tbody></table>
<table><caption>Table S3: Population fixation indices calculated on real and simulated data</caption><thead><tr><th></th><th>Population</th><th>F data</th><th>F simulations</th></tr></thead><tbody><tr><td rowspan="4">centering</td><td>Angus</td><td>0.174</td><td>0.162</td></tr><tr><td>Fleckvieh</td><td>0.127</td><td>0.092</td></tr><tr><td>Holstein</td><td>0.106</td><td>0.106</td></tr><tr><td>Jersey</td><td>0.225</td><td>0.184</td></tr></tbody></table>
<table><thead><tr><th rowspan="2">position</th><th rowspan="2">P-value (FLK)</th><th colspan="4">Allele frequencies</th><th rowspan="2">QTN</th></tr><tr><th>Holstein</th><th>Angus</th><th>Fleckvieh</th><th>Jersey</th></tr></thead><tbody><tr><td>24805830</td><td>7.9 10<sup>-4</sup></td><td>0.87</td><td>0.96</td><td>0.12</td><td>0.03</td><td>No</td></tr><tr><td>24973953*</td><td>8.1 10<sup>-5</sup></td><td>1.0</td><td>0.96</td><td>0.10</td><td>0.03</td><td>Yes</td></tr><tr><td>24974221*</td><td>1.0 10<sup>-4</sup></td><td>1.0</td><td>0.96</td><td>0.12</td><td>0.03</td><td>Yes</td></tr><tr><td>24974811*</td><td>7.9 10<sup>-5</sup></td><td>1.0</td><td>0.96</td><td>0.12</td><td>0.0</td><td>Yes</td></tr><tr><td>24991209</td><td>6.3 10<sup>-5</sup></td><td>1.0</td><td>0.96</td><td>0.10</td><td>0.0</td><td>No</td></tr><tr><td>25003338*</td><td>6.3 10<sup>-5</sup></td><td>1.0</td><td>0.96</td><td>0.10</td><td>0.0</td><td>Yes</td></tr><tr><td>25006125*</td><td>3.9 10<sup>-5</sup></td><td>1.0</td><td>0.96</td><td>0.06</td><td>0.0</td><td>Yes</td></tr><tr><td>25015640*</td><td>3.1 10<sup>-5</sup></td><td>1.0</td><td>0.96</td><td>0.04</td><td>0.0</td><td>Yes</td></tr><tr><td>25019900</td><td>1.6 10<sup>-3</sup></td><td>1.0</td><td>0.96</td><td>0.56</td><td>0.07</td><td>No†</td></tr><tr><td>25022815</td><td>6.8 10<sup>-4</sup></td><td>0.99</td><td>0.96</td><td>0.60</td><td>0.0</td><td>No†</td></tr><tr><td>25031172</td><td>9.0 10<sup>-4</sup></td><td>1</td><td>0.97</td><td>0.54</td><td>0.03</td><td>No†</td></tr><tr><td>25052396</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>Yes</td></tr><tr><td>25052440</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>Yes</td></tr><tr><td>25059742</td><td>7.6 10<sup>-4</sup></td><td>1</td><td>0.96</td><td>0.54</td><td>0.0</td><td>No†</td></tr><tr><td>25061179</td><td>7.6 10<sup>-4</sup></td><td>1.00</td><td>0.96</td><td>0.54</td><td>0.0</td><td>No†</td></tr></tbody></table>

† : possible QTN in the Holstein × Jersey cross, but ruled out as QTN in an association study in the Fleckvieh population [12].

* : Mutations that are still candidates after combining selection tests and QTN results.

<caption>Table S4: Candidate selected mutations and/or QTN (from [12]) at the selection signature around the PLAG1 gene on BTA14</caption>
<table><caption>Table S5: Candidate mutations in private sweep regions. Horizontal lines are used to group closely related sweep windows, which might result from the same selection event.</caption><thead><tr><th>chr</th><th>position</th><th>pval</th><th>Jersey</th><th>Angus</th><th>Fleckvieh</th><th>Holstein</th><th>annotation</th></tr></thead><tbody><tr><td>6</td><td>71517158</td><td>1</td><td>0</td><td>0</td><td>0.92</td><td>0.02</td><td>intergenic</td></tr><tr><td>6</td><td>71552373</td><td>1</td><td>0</td><td>0</td><td>0.92</td><td>0</td><td>intergenic</td></tr><tr><td>7</td><td>25574759</td><td>5</td><td>0.77</td><td>0.01</td><td>0</td><td>0</td><td>intergenic</td></tr><tr><td>7</td><td>25580485</td><td>7</td><td>0.8</td><td>0.01</td><td>0</td><td>0</td><td>intergenic</td></tr><tr><td>7</td><td>25581841</td><td>7</td><td>0.8</td><td>0.01</td><td>0</td><td>0</td><td>intergenic</td></tr><tr><td>7</td><td>25582672</td><td>7</td><td>0.8</td><td>0.01</td><td>0</td><td>0</td><td>downstream KIAA1024L</td></tr><tr><td>7</td><td>25582815</td><td>7</td><td>0.8</td><td>0.01</td><td>0</td><td>0</td><td>downstream KIAA1024L</td></tr><tr><td>7</td><td>25589078</td><td>7</td><td>0.8</td><td>0.01</td><td>0</td><td>0</td><td>intron KIAA1024L</td></tr><tr><td>7</td><td>25673520</td><td>3</td><td>0.77</td><td>0.01</td><td>0</td><td>0</td><td>intron ADAMTS19</td></tr><tr><td>7</td><td>25756290</td><td>3</td><td>0.77</td><td>0</td><td>0</td><td>0</td><td>intron ADAMTS19</td></tr><tr><td>7</td><td>26335268</td><td>2</td><td>0.8</td><td>0</td><td>0</td><td>0</td><td>intergenic</td></tr><tr><td>7</td><td>26459812</td><td>1</td><td>0.83</td><td>0</td><td>0</td><td>0</td><td>intergenic</td></tr><tr><td>7</td><td>26472093</td><td>3</td><td>0.83</td><td>0.01</td><td>0.02</td><td>0</td><td>intergenic</td></tr><tr><td>7</td><td>26523973</td><td>3</td><td>0.83</td><td>0.01</td><td>0.01</td><td>0</td><td>intergenic</td></tr><tr><td>7</td><td>26740228</td><td>3</td><td>0.77</td><td>0</td><td>0</td><td>0</td><td>intron FBN2</td></tr><tr><td>7</td><td>26743727</td><td>7</td><td>0.8</td><td>0</td><td>0.06</td><td>0</td><td>intron FBN2</td></tr><tr><td>7</td><td>26744241</td><td>7</td><td>0.8</td><td>0</td><td>0.06</td><td>0</td><td>intron FBN2</td></tr><tr><td>7</td><td>26744717</td><td>3</td><td>0.77</td><td>0</td><td>0</td><td>0</td><td>intron FBN2</td></tr><tr><td>7</td><td>26753548</td><td>2</td><td>0.8</td><td>0</td><td>0</td><td>0</td><td>intron FBN2</td></tr><tr><td>7</td><td>26761008</td><td>2</td><td>0.8</td><td>0</td><td>0</td><td>0</td><td>intron FBN2</td></tr><tr><td>7</td><td>26762977</td><td>2</td><td>0.8</td><td>0</td><td>0</td><td>0</td><td>intron FBN2</td></tr><tr><td>7</td><td>26768920</td><td>3</td><td>0.77</td><td>0</td><td>0</td><td>0</td><td>intron FBN2</td></tr><tr><td>7</td><td>26772586</td><td>2</td><td>0.8</td><td>0</td><td>0</td><td>0</td><td>intron FBN2</td></tr><tr><td>7</td><td>26776782</td><td>3</td><td>0.77</td><td>0</td><td>0</td><td>0</td><td>intron FBN2</td></tr><tr><td>7</td><td>26789277</td><td>2</td><td>0.8</td><td>0</td><td>0</td><td>0</td><td>intron FBN2</td></tr><tr><td>7</td><td>26871380</td><td>7</td><td>0.77</td><td>0</td><td>0</td><td>0.01</td><td>intron FBN2</td></tr><tr><td>14</td><td>24973953</td><td>8</td><td>0.97</td><td>0.04</td><td>0.86</td><td>0</td><td>downstream MOS</td></tr><tr><td>14</td><td>24974811</td><td>8</td><td>1</td><td>0.04</td><td>0.85</td><td>0</td><td>downstream MOS</td></tr><tr><td>14</td><td>24991209</td><td>6</td><td>1</td><td>0.04</td><td>0.86</td><td>0</td><td>intergenic</td></tr><tr><td>14</td><td>24995794</td><td>37</td><td>0.83</td><td>0.04</td><td>0.79</td><td>0</td><td>intergenic</td></tr><tr><td>14</td><td>25003338</td><td>6</td><td>1</td><td>0.04</td><td>0.86</td><td>0.01</td><td>downstream PLAG1</td></tr><tr><td>14</td><td>25006125</td><td>4</td><td>1</td><td>0.04</td><td>0.91</td><td>0.01</td><td>downstream PLAG1</td></tr><tr><td>14</td><td>25015640</td><td>3</td><td>1</td><td>0.04</td><td>0.92</td><td>0.02</td><td>intergenic</td></tr><tr><td>18</td><td>14757910</td><td>4</td><td>0.93</td><td>0.02</td><td>0.99</td><td>0.07</td><td>missense MC1R</td></tr><tr><td>18</td><td>14757923</td><td>1</td><td>0.03</td><td>0.02</td><td>0.95</td><td>0.04</td><td>frameshift MC1R</td></tr><tr><td>18</td><td>14843827</td><td>8</td><td>0.07</td><td>0.02</td><td>0.86</td><td>0.08</td><td>intron GAS8</td></tr><tr><td>20</td><td>24307385</td><td>3</td><td>0.77</td><td>0</td><td>0</td><td>0</td><td>intergenic</td></tr><tr><td>20</td><td>25255113</td><td>7</td><td>0.8</td><td>0.04</td><td>0</td><td>0</td><td>intergenic</td></tr></tbody></table>

*Continued on next page*
<table><thead><tr><th>chr</th><th>position</th><th>pval</th><th>Jersey</th><th>Angus</th><th>Fleckvieh</th><th>Holstein</th><th>annotation</th></tr></thead><tbody><tr><td>20</td><td>25267029</td><td>2</td><td>0.8</td><td>0</td><td>0</td><td>0</td><td>intergenic</td></tr><tr><td>20</td><td>25267814</td><td>2</td><td>0.8</td><td>0</td><td>0</td><td>0</td><td>intergenic</td></tr><tr><td>20</td><td>25287735</td><td>3</td><td>0.77</td><td>0</td><td>0</td><td>0</td><td>intergenic</td></tr><tr><td>20</td><td>25287759</td><td>3</td><td>0.77</td><td>0</td><td>0</td><td>0</td><td>intergenic</td></tr><tr><td>20</td><td>25297390</td><td>3</td><td>0.77</td><td>0</td><td>0</td><td>0</td><td>intergenic</td></tr><tr><td>20</td><td>25306614</td><td>2</td><td>0.8</td><td>0</td><td>0</td><td>0</td><td>intergenic</td></tr><tr><td>20</td><td>25308282</td><td>2</td><td>0.8</td><td>0</td><td>0</td><td>0</td><td>intergenic</td></tr><tr><td>20</td><td>25308605</td><td>3</td><td>0.77</td><td>0</td><td>0</td><td>0</td><td>intergenic</td></tr><tr><td>20</td><td>25312532</td><td>2</td><td>0.8</td><td>0</td><td>0</td><td>0</td><td>intergenic</td></tr><tr><td>20</td><td>25313689</td><td>3</td><td>0.77</td><td>0</td><td>0</td><td>0</td><td>intergenic</td></tr><tr><td>20</td><td>25314539</td><td>3</td><td>0.77</td><td>0</td><td>0</td><td>0</td><td>intergenic</td></tr><tr><td>20</td><td>25317501</td><td>2</td><td>0.8</td><td>0</td><td>0</td><td>0</td><td>intergenic</td></tr><tr><td>20</td><td>25318166</td><td>3</td><td>0.77</td><td>0</td><td>0</td><td>0</td><td>intergenic</td></tr><tr><td>20</td><td>25319515</td><td>3</td><td>0.77</td><td>0</td><td>0</td><td>0</td><td>intergenic</td></tr><tr><td>20</td><td>25320693</td><td>2</td><td>0.8</td><td>0</td><td>0</td><td>0</td><td>intergenic</td></tr><tr><td>20</td><td>25320734</td><td>2</td><td>0.8</td><td>0</td><td>0</td><td>0</td><td>intergenic</td></tr><tr><td>20</td><td>25321524</td><td>2</td><td>0.8</td><td>0</td><td>0</td><td>0</td><td>intergenic</td></tr><tr><td>20</td><td>25321720</td><td>2</td><td>0.8</td><td>0</td><td>0</td><td>0</td><td>intergenic</td></tr><tr><td>20</td><td>25321795</td><td>2</td><td>0.8</td><td>0</td><td>0</td><td>0</td><td>intergenic</td></tr><tr><td>20</td><td>25323147</td><td>3</td><td>0.77</td><td>0</td><td>0</td><td>0</td><td>intergenic</td></tr><tr><td>20</td><td>26912317</td><td>3</td><td>0.77</td><td>0</td><td>0</td><td>0</td><td>intergenic</td></tr><tr><td>20</td><td>27147204</td><td>3</td><td>0.77</td><td>0</td><td>0</td><td>0</td><td>intergenic</td></tr><tr><td>20</td><td>39872347</td><td>11</td><td>0.8</td><td>0.01</td><td>0.03</td><td>0</td><td>downstream SLC45A2<br>downstream RXFP3</td></tr><tr><td>22</td><td>35680794</td><td>3</td><td>0.77</td><td>0</td><td>0</td><td>0</td><td>upstream bta-mir</td></tr><tr><td>22</td><td>35970263</td><td>5</td><td>0.9</td><td>0.14</td><td>0.02</td><td>0.01</td><td>intergenic</td></tr><tr><td>22</td><td>35973490</td><td>8</td><td>0.93</td><td>0.23</td><td>0.03</td><td>0.02</td><td>intergenic</td></tr><tr><td>22</td><td>36011838</td><td>42</td><td>0.87</td><td>0.24</td><td>0.03</td><td>0.02</td><td>intron MAGI1</td></tr><tr><td>22</td><td>36017262</td><td>15</td><td>0.87</td><td>0.11</td><td>0.03</td><td>0.01</td><td>intron MAGI1</td></tr><tr><td>22</td><td>36019024</td><td>25</td><td>0.9</td><td>0.24</td><td>0.03</td><td>0.02</td><td>intron MAGI1</td></tr><tr><td>22</td><td>36022130</td><td>48</td><td>0.87</td><td>0.26</td><td>0.03</td><td>0.02</td><td>intron MAGI1</td></tr><tr><td>22</td><td>36022963</td><td>45</td><td>0.9</td><td>0.26</td><td>0.05</td><td>0.02</td><td>intron MAGI1</td></tr><tr><td>27</td><td>4220866</td><td>10</td><td>0.9</td><td>0.07</td><td>0.06</td><td>0.02</td><td>intergenic</td></tr></tbody></table>
<table><caption>Table S6: Candidate polymorphisms within shared sweeps.</caption><thead><tr><th>chr</th><th>pos</th><th>major</th><th>minor</th><th>Jersey</th><th>Angus</th><th>Fleckvieh</th><th>Holstein</th><th>yak</th><th>buffalo</th><th>sheep</th><th># yak</th><th>annotation</th></tr></thead><tbody><tr><td>1</td><td>1784083</td><td>G</td><td>A</td><td>0</td><td>0</td><td>0.01</td><td>0</td><td>A</td><td>-</td><td>A</td><td>8</td><td>intergenic</td></tr><tr><td>1</td><td>1802301</td><td>G</td><td>-</td><td>0</td><td>0</td><td>0</td><td>0</td><td>T</td><td>T</td><td>T</td><td>9</td><td>intergenic</td></tr><tr><td>1</td><td>1808680</td><td>C</td><td>-</td><td>0</td><td>0</td><td>0</td><td>0</td><td>T</td><td>-</td><td>T</td><td>8</td><td>intergenic</td></tr><tr><td>1</td><td>107502158</td><td>G</td><td>-</td><td>0</td><td>0</td><td>0</td><td>0</td><td>A</td><td>A</td><td>A</td><td>8</td><td>intron PPM1L</td></tr><tr><td>1</td><td>107502160</td><td>C</td><td>-</td><td>0</td><td>0</td><td>0</td><td>0</td><td>A</td><td>A</td><td>A</td><td>8</td><td>intron PPM1L</td></tr><tr><td>1</td><td>107550839</td><td>T</td><td>-</td><td>0</td><td>0</td><td>0</td><td>0</td><td>C</td><td>C</td><td>-</td><td>8</td><td>intron PPM1L</td></tr><tr><td>1</td><td>107576115</td><td>A</td><td>-</td><td>0</td><td>0</td><td>0</td><td>0</td><td>C</td><td>C</td><td>C</td><td>10</td><td>upstream PPM1L</td></tr><tr><td>1</td><td>107583661</td><td>G</td><td>-</td><td>0</td><td>0</td><td>0</td><td>0</td><td>A</td><td>A</td><td>A</td><td>10</td><td>intergenic</td></tr><tr><td>1</td><td>107585004</td><td>A</td><td>-</td><td>0</td><td>0</td><td>0</td><td>0</td><td>G</td><td>G</td><td>G</td><td>8</td><td>intergenic</td></tr><tr><td>1</td><td>107645676</td><td>G</td><td>C</td><td>0</td><td>0</td><td>0</td><td>0.02</td><td>C</td><td>C</td><td>-</td><td>8</td><td>intergenic</td></tr><tr><td>1</td><td>107646230</td><td>T</td><td>C</td><td>0</td><td>0</td><td>0</td><td>0.01</td><td>C</td><td>C</td><td>C</td><td>10</td><td>intergenic</td></tr><tr><td>1</td><td>107688061</td><td>G</td><td>-</td><td>0</td><td>0</td><td>0</td><td>0</td><td>A</td><td>-</td><td>A</td><td>8</td><td>synonymous ARL14</td></tr><tr><td>1</td><td>107688282</td><td>C</td><td>-</td><td>0</td><td>0</td><td>0</td><td>0</td><td>T</td><td>-</td><td>T</td><td>8</td><td>upstream ARL14</td></tr><tr><td>1</td><td>107713763</td><td>A</td><td>-</td><td>0</td><td>0</td><td>0</td><td>0</td><td>G</td><td>G</td><td>G</td><td>10</td><td>intergenic</td></tr><tr><td>1</td><td>107719968</td><td>G</td><td>A</td><td>0.03</td><td>0</td><td>0.02</td><td>0.04</td><td>A</td><td>-</td><td>A</td><td>8</td><td>intergenic</td></tr><tr><td>1</td><td>107720085</td><td>C</td><td>-</td><td>0</td><td>0</td><td>0</td><td>0</td><td>A</td><td>A</td><td>A</td><td>9</td><td>intergenic</td></tr><tr><td>1</td><td>107720131</td><td>T</td><td>-</td><td>0</td><td>0</td><td>0</td><td>0</td><td>C</td><td>C</td><td>C</td><td>8</td><td>intergenic</td></tr><tr><td>1</td><td>107736632</td><td>C</td><td>A</td><td>0</td><td>0.01</td><td>0.02</td><td>0.06</td><td>A</td><td>A</td><td>A</td><td>8</td><td>intergenic</td></tr><tr><td>5</td><td>68689789</td><td>T</td><td>A</td><td>0</td><td>0.02</td><td>0.02</td><td>0</td><td>A</td><td>-</td><td>A</td><td>8</td><td>intergenic</td></tr><tr><td>5</td><td>68706351</td><td>T</td><td>-</td><td>0</td><td>0</td><td>0</td><td>0</td><td>C</td><td>-</td><td>C</td><td>8</td><td>intron SLC41A2</td></tr><tr><td>5</td><td>68712273</td><td>A</td><td>C</td><td>0</td><td>0.02</td><td>0.09</td><td>0</td><td>C</td><td>-</td><td>C</td><td>8</td><td>intron SLC41A2</td></tr><tr><td>5</td><td>68758455</td><td>A</td><td>-</td><td>0</td><td>0</td><td>0</td><td>0</td><td>G</td><td>-</td><td>G</td><td>8</td><td>intron SLC41A2</td></tr></tbody></table>

*Continued on next page*
<table><thead><tr><th>chr</th><th>pos</th><th>major</th><th>minor</th><th>Jersey</th><th>Angus</th><th>Fleckvieh</th><th>Holstein</th><th>yak</th><th>buffalo</th><th>ovis</th><th># yak</th><th>annotation</th></tr></thead><tbody><tr><td>5</td><td>68783378</td><td>C</td><td>T</td><td>0</td><td>0</td><td>0.07</td><td>0</td><td>T</td><td>T</td><td>T</td><td>8</td><td>splice region SLC41A2, intron SLC41A2</td></tr><tr><td>7</td><td>4574596</td><td>T</td><td>C</td><td>0.03</td><td>0</td><td>0</td><td>0.04</td><td>C</td><td>-</td><td>C</td><td>8</td><td>3 prime UTR FKBP8</td></tr><tr><td>7</td><td>4671036</td><td>T</td><td>C</td><td>0</td><td>0</td><td>0.12</td><td>0</td><td>C</td><td>-</td><td>C</td><td>8</td><td>missense ISYNA1, downstream SSBP4</td></tr><tr><td>7</td><td>4676230</td><td>A</td><td>-</td><td>0</td><td>0</td><td>0</td><td>0</td><td>G</td><td>-</td><td>G</td><td>8</td><td>downstream ISYNA1, intron SSBP4</td></tr><tr><td>7</td><td>4676874</td><td>T</td><td>-</td><td>0</td><td>0</td><td>0</td><td>0</td><td>C</td><td>-</td><td>C</td><td>8</td><td>downstream ISYNA1, missense SSBP4</td></tr><tr><td>10</td><td>59182584</td><td>G</td><td>A</td><td>0.03</td><td>0.06</td><td>0</td><td>0.1</td><td>A</td><td>-</td><td>A</td><td>8</td><td>intergenic</td></tr><tr><td>10</td><td>59189890</td><td>C</td><td>-</td><td>0</td><td>0</td><td>0</td><td>0</td><td>T</td><td>T</td><td>T</td><td>8</td><td>intergenic</td></tr><tr><td>10</td><td>59192974</td><td>C</td><td>T</td><td>0.03</td><td>0.06</td><td>0</td><td>0.04</td><td>T</td><td>T</td><td>-</td><td>8</td><td>intergenic</td></tr><tr><td>10</td><td>59193189</td><td>T</td><td>-</td><td>0</td><td>0</td><td>0</td><td>0</td><td>C</td><td>-</td><td>C</td><td>8</td><td>intergenic</td></tr><tr><td>10</td><td>59193664</td><td>A</td><td>-</td><td>0</td><td>0</td><td>0</td><td>0</td><td>G</td><td>-</td><td>G</td><td>8</td><td>intergenic</td></tr><tr><td>10</td><td>59202522</td><td>G</td><td>-</td><td>0</td><td>0</td><td>0</td><td>0</td><td>A</td><td>A</td><td>A</td><td>9</td><td>intergenic</td></tr><tr><td>10</td><td>59203717</td><td>G</td><td>-</td><td>0</td><td>0</td><td>0</td><td>0</td><td>A</td><td>-</td><td>A</td><td>8</td><td>intergenic</td></tr><tr><td>10</td><td>59204098</td><td>G</td><td>-</td><td>0</td><td>0</td><td>0</td><td>0</td><td>A</td><td>-</td><td>A</td><td>8</td><td>intergenic</td></tr><tr><td>10</td><td>59205092</td><td>T</td><td>-</td><td>0</td><td>0</td><td>0</td><td>0</td><td>C</td><td>C</td><td>-</td><td>8</td><td>intergenic</td></tr><tr><td>10</td><td>59206693</td><td>G</td><td>-</td><td>0</td><td>0</td><td>0</td><td>0</td><td>T</td><td>-</td><td>T</td><td>8</td><td>intergenic</td></tr><tr><td>10</td><td>59248594</td><td>G</td><td>-</td><td>0</td><td>0</td><td>0</td><td>0</td><td>A</td><td>A</td><td>A</td><td>8</td><td>intron CYP19A1</td></tr><tr><td>10</td><td>59274258</td><td>T</td><td>G</td><td>0.03</td><td>0.06</td><td>0</td><td>0.08</td><td>G</td><td>A</td><td>A</td><td>8</td><td>intron CYP19A1</td></tr><tr><td>10</td><td>59308568</td><td>G</td><td>-</td><td>0</td><td>0</td><td>0</td><td>0</td><td>T</td><td>T</td><td>T</td><td>8</td><td>intergenic</td></tr><tr><td>10</td><td>59314493</td><td>A</td><td>-</td><td>0</td><td>0</td><td>0</td><td>0</td><td>T</td><td>T</td><td>-</td><td>8</td><td>intergenic</td></tr><tr><td>10</td><td>59319111</td><td>G</td><td>-</td><td>0</td><td>0</td><td>0</td><td>0</td><td>A</td><td>-</td><td>A</td><td>8</td><td>intergenic</td></tr></tbody></table>

*Continued on next page*
<table><thead><tr><th>chr</th><th>pos</th><th>major</th><th>minor</th><th>Jersey</th><th>Angus</th><th>Fleckvieh</th><th>Holstein</th><th>yak</th><th>buffalo</th><th>ovis</th><th># yak</th><th>annotation</th></tr></thead><tbody><tr><td>10</td><td>59326268</td><td>C</td><td>T</td><td>0</td><td>0</td><td>0</td><td>0</td><td>T</td><td>-</td><td>T</td><td>8</td><td>intergenic</td></tr><tr><td>10</td><td>59330483</td><td>G</td><td>-</td><td>0</td><td>0</td><td>0</td><td>0</td><td>A</td><td>-</td><td>A</td><td>8</td><td>intergenic</td></tr><tr><td>16</td><td>44674060</td><td>A</td><td>-</td><td>0</td><td>0</td><td>0</td><td>0</td><td>G</td><td>-</td><td>G</td><td>8</td><td>downstream CLSTN1, downstream PK3CD</td></tr><tr><td>16</td><td>44679395</td><td>T</td><td>-</td><td>0</td><td>0</td><td>0</td><td>0</td><td>A</td><td>-</td><td>A</td><td>8</td><td>missense PK3CD</td></tr><tr><td>16</td><td>44680925</td><td>A</td><td>-</td><td>0</td><td>0</td><td>0</td><td>0</td><td>G</td><td>G</td><td>G</td><td>10</td><td>missense PK3CD</td></tr><tr><td>16</td><td>44681107</td><td>T</td><td>-</td><td>0</td><td>0</td><td>0</td><td>0</td><td>C</td><td>C</td><td>C</td><td>9</td><td>intron PK3CD, upstream U6</td></tr><tr><td>16</td><td>44681560</td><td>T</td><td>-</td><td>0</td><td>0</td><td>0</td><td>0</td><td>C</td><td>-</td><td>C</td><td>8</td><td>intron PK3CD, upstream U6</td></tr><tr><td>16</td><td>44682130</td><td>A</td><td>G</td><td>0</td><td>0</td><td>0.02</td><td>0</td><td>G</td><td>G</td><td>G</td><td>9</td><td>intron PK3CD, upstream U6</td></tr><tr><td>16</td><td>44683730</td><td>A</td><td>-</td><td>0</td><td>0</td><td>0</td><td>0</td><td>G</td><td>-</td><td>G</td><td>8</td><td>synonymous PK3CD, upstream U6</td></tr><tr><td>16</td><td>44686127</td><td>A</td><td>-</td><td>0</td><td>0</td><td>0</td><td>0</td><td>G</td><td>-</td><td>G</td><td>8</td><td>intron PK3CD, downstream U6</td></tr><tr><td>16</td><td>44725237</td><td>T</td><td>G</td><td>0</td><td>0</td><td>0</td><td>0.06</td><td>C</td><td>C</td><td>C</td><td>9</td><td>upstream 5S_rRNA</td></tr><tr><td>16</td><td>44753716</td><td>C</td><td>-</td><td>0</td><td>0</td><td>0</td><td>0</td><td>T</td><td>-</td><td>T</td><td>8</td><td>intergenic</td></tr><tr><td>16</td><td>44777756</td><td>T</td><td>-</td><td>0</td><td>0</td><td>0</td><td>0</td><td>G</td><td>-</td><td>G</td><td>8</td><td>synonymous TMEM201</td></tr><tr><td>16</td><td>44806651</td><td>G</td><td>-</td><td>0</td><td>0</td><td>0</td><td>0</td><td>A</td><td>A</td><td>A</td><td>9</td><td>intron SLC25A33</td></tr><tr><td>16</td><td>44832764</td><td>A</td><td>-</td><td>0</td><td>0</td><td>0</td><td>0</td><td>G</td><td>G</td><td>G</td><td>8</td><td>upstream SLC25A33</td></tr><tr><td>16</td><td>44837088</td><td>G</td><td>-</td><td>0</td><td>0</td><td>0</td><td>0</td><td>A</td><td>-</td><td>A</td><td>8</td><td>intergenic</td></tr><tr><td>16</td><td>44879916</td><td>T</td><td>C</td><td>0.03</td><td>0.05</td><td>0.12</td><td>0.02</td><td>C</td><td>C</td><td>C</td><td>9</td><td>intergenic</td></tr><tr><td>16</td><td>44897589</td><td>A</td><td>-</td><td>0</td><td>0</td><td>0</td><td>0</td><td>G</td><td>-</td><td>G</td><td>8</td><td>intergenic</td></tr><tr><td>16</td><td>44899871</td><td>A</td><td>-</td><td>0</td><td>0</td><td>0</td><td>0</td><td>G</td><td>-</td><td>G</td><td>8</td><td>intergenic</td></tr><tr><td>16</td><td>44917105</td><td>T</td><td>-</td><td>0</td><td>0</td><td>0</td><td>0</td><td>C</td><td>-</td><td>C</td><td>9</td><td>intergenic</td></tr><tr><td>16</td><td>44952201</td><td>A</td><td>-</td><td>0</td><td>0</td><td>0</td><td>0</td><td>G</td><td>G</td><td>G</td><td>9</td><td>intergenic</td></tr></tbody></table>

*Continued on next page*
<table><thead><tr><th>chr</th><th>pos</th><th>major</th><th>minor</th><th>Jersey</th><th>Angus</th><th>Fleckvieh</th><th>Holstein</th><th>yak</th><th>buffalo</th><th>ovis</th><th># yak</th><th>annotation</th></tr></thead><tbody><tr><td>16</td><td>45663151</td><td>C</td><td>-</td><td>0</td><td>0</td><td>0</td><td>0</td><td>A</td><td>-</td><td>A</td><td>9</td><td>intron RERE</td></tr><tr><td>16</td><td>45663204</td><td>T</td><td>-</td><td>0</td><td>0</td><td>0</td><td>0</td><td>C</td><td>-</td><td>C</td><td>8</td><td>intron RERE</td></tr><tr><td>16</td><td>45663289</td><td>C</td><td>-</td><td>0</td><td>0</td><td>0</td><td>0</td><td>T</td><td>-</td><td>T</td><td>8</td><td>intron RERE</td></tr><tr><td>16</td><td>45680349</td><td>T</td><td>-</td><td>0</td><td>0</td><td>0</td><td>0</td><td>C</td><td>-</td><td>C</td><td>8</td><td>intron RERE</td></tr><tr><td>16</td><td>45707230</td><td>A</td><td>-</td><td>0</td><td>0</td><td>0</td><td>0</td><td>G</td><td>G</td><td>G</td><td>8</td><td>intron RERE</td></tr><tr><td>16</td><td>45708963</td><td>T</td><td>C</td><td>0.03</td><td>0.04</td><td>0.05</td><td>0.11</td><td>C</td><td>C</td><td>C</td><td>9</td><td>intron RERE</td></tr><tr><td>16</td><td>45711415</td><td>T</td><td>-</td><td>0</td><td>0</td><td>0</td><td>0</td><td>C</td><td>-</td><td>C</td><td>8</td><td>splice region RERE, synonymous RERE</td></tr><tr><td>16</td><td>45729243</td><td>G</td><td>-</td><td>0</td><td>0</td><td>0</td><td>0</td><td>A</td><td>-</td><td>A</td><td>8</td><td>intron RERE</td></tr><tr><td>16</td><td>45762651</td><td>A</td><td>-</td><td>0</td><td>0</td><td>0</td><td>0</td><td>G</td><td>-</td><td>G</td><td>8</td><td>intron RERE</td></tr><tr><td>16</td><td>45773132</td><td>A</td><td>G</td><td>0.03</td><td>0.03</td><td>0.05</td><td>0.13</td><td>G</td><td>G</td><td>G</td><td>8</td><td>intron RERE</td></tr><tr><td>16</td><td>45773382</td><td>A</td><td>-</td><td>0</td><td>0</td><td>0</td><td>0</td><td>T</td><td>-</td><td>T</td><td>8</td><td>intron RERE</td></tr><tr><td>16</td><td>45782314</td><td>G</td><td>-</td><td>0</td><td>0</td><td>0</td><td>0</td><td>A</td><td>-</td><td>A</td><td>9</td><td>intron RERE</td></tr><tr><td>16</td><td>45782934</td><td>G</td><td>-</td><td>0</td><td>0</td><td>0</td><td>0</td><td>A</td><td>A</td><td>A</td><td>10</td><td>intron RERE</td></tr><tr><td>16</td><td>45784480</td><td>T</td><td>-</td><td>0</td><td>0</td><td>0</td><td>0</td><td>C</td><td>-</td><td>C</td><td>8</td><td>intron RERE</td></tr><tr><td>16</td><td>45817029</td><td>G</td><td>C</td><td>0</td><td>0.04</td><td>0.05</td><td>0.13</td><td>C</td><td>C</td><td>C</td><td>9</td><td>intron RERE</td></tr><tr><td>16</td><td>45817370</td><td>C</td><td>-</td><td>0</td><td>0</td><td>0</td><td>0</td><td>T</td><td>-</td><td>T</td><td>8</td><td>intron RERE</td></tr><tr><td>16</td><td>45828780</td><td>T</td><td>-</td><td>0</td><td>0</td><td>0</td><td>0</td><td>C</td><td>-</td><td>C</td><td>8</td><td>intron RERE</td></tr><tr><td>16</td><td>45857131</td><td>G</td><td>A</td><td>0</td><td>0.02</td><td>0.01</td><td>0.11</td><td>A</td><td>A</td><td>A</td><td>8</td><td>intron RERE</td></tr><tr><td>16</td><td>45859900</td><td>C</td><td>G</td><td>0</td><td>0.02</td><td>0.02</td><td>0.12</td><td>G</td><td>-</td><td>G</td><td>8</td><td>intron RERE</td></tr><tr><td>16</td><td>45880729</td><td>T</td><td>-</td><td>0</td><td>0</td><td>0</td><td>0</td><td>C</td><td>-</td><td>C</td><td>8</td><td>downstream RERE</td></tr><tr><td>16</td><td>45897447</td><td>T</td><td>-</td><td>0</td><td>0</td><td>0</td><td>0</td><td>C</td><td>-</td><td>C</td><td>8</td><td>synonymous SLC45A1</td></tr></tbody></table>
<table><caption>Table S7: Genome regions exhibiting high levels of haplotype differentiation between breeds based on the hapFLK statistic.</caption><thead><tr><th>Id</th><th>Chr</th><th>begin (Mb)</th><th>end (Mb)</th><th>p-value</th><th>breed(s)</th><th>Candidate genes</th><th>References</th><th>Annotation</th></tr></thead><tbody><tr><td>1</td><td>1</td><td>26.884</td><td>26.926</td><td>1.4 10<sup>-6</sup></td><td>ANG</td><td>5Kb upstream of ROBO1 promoter region</td><td>[4, 5]</td><td>Follicular development</td></tr><tr><td>2</td><td>1</td><td>76.641</td><td>76.649</td><td>3.8 10<sup>-5</sup></td><td></td><td>CCDC50, OSTN</td><td></td><td></td></tr><tr><td>3</td><td>1</td><td>151.721</td><td>151.723</td><td>7.0 10<sup>-5</sup></td><td>HOL</td><td>KCNJ6 (syn. GIRK2)</td><td></td><td>rs41755372, weaver gene (mice)</td></tr><tr><td>4</td><td>2</td><td>128.400</td><td>128.450</td><td>1.9 10<sup>-5</sup></td><td>FLE, (JER, HOL)</td><td>RUNX3, RHCE region,</td><td>[31, 27]</td><td>Bone development (chondrocytes)</td></tr><tr><td>5</td><td>3</td><td>39.591</td><td>39.594</td><td>3.1 10<sup>-5</sup></td><td></td><td></td><td></td><td></td></tr><tr><td>6</td><td>3</td><td>53.727</td><td>54.414</td><td>1.5 10<sup>-5</sup></td><td>ANG, HOL, JER</td><td>LRRC8C (FAD158), LLRC8D, GBP5, GBP6</td><td>[33, 14, 9]</td><td>Adiposity (LRRC8C), Meat tenderness (GBP5), Embryo implantation (GBP5)</td></tr><tr><td>7</td><td>3</td><td>86.508</td><td>86.786</td><td>2.2 10<sup>-5</sup></td><td></td><td>HOOK1, Cytochrome P450 (LOC521656)</td><td></td><td></td></tr><tr><td>8</td><td>4</td><td>38.431</td><td>38.444</td><td>2.7 10<sup>-5</sup></td><td>ANG</td><td>CACNA2D1</td><td></td><td></td></tr><tr><td>9</td><td>4</td><td>41.745</td><td>41.749</td><td>6.5 10<sup>-5</sup></td><td></td><td>GNAI1, MAGI2</td><td></td><td></td></tr><tr><td>10</td><td>4</td><td>50.029</td><td>50.333</td><td>2.0 10<sup>-6</sup></td><td>FLE, JER</td><td>STARD3NL</td><td>[22]</td><td>Bone Mineral Density</td></tr><tr><td>11</td><td>4</td><td>76.939</td><td>76.946</td><td>7.6 10<sup>-5</sup></td><td></td><td>ADCY1, MYO1G</td><td></td><td></td></tr><tr><td>12</td><td>4</td><td>77.578</td><td>77.644</td><td>5.1 10<sup>-6</sup></td><td></td><td>NUDCD3, NPC1L1</td><td></td><td></td></tr><tr><td>13</td><td>4</td><td>94.791</td><td>94.797</td><td>2.7 10<sup>-5</sup></td><td></td><td>TMEM209</td><td></td><td></td></tr><tr><td>14</td><td>4</td><td>116.144</td><td>116.151</td><td>3.9 10<sup>-5</sup></td><td></td><td></td><td></td><td></td></tr></tbody></table>

*Continued on next page*
<table><thead><tr><th>Id</th><th>Chr</th><th>begin (Mb)</th><th>end (Mb)</th><th>p-value</th><th>breed(s)</th><th>Candidate genes</th><th>References</th><th>Annotation</th></tr></thead><tbody><tr><td>15</td><td>6</td><td>38.622</td><td>38.726</td><td>5.3 10<sup>-5</sup></td><td>FLE, (JER)</td><td>NCAPG, LCORL</td><td>[6, 3]</td><td>Fetal Growth</td></tr><tr><td>16</td><td>6</td><td>66.810</td><td>66.810</td><td>9.7 10<sup>-5</sup></td><td></td><td></td><td></td><td></td></tr><tr><td>17</td><td>6</td><td>70.332</td><td>71.607</td><td>1.2 10<sup>-12</sup></td><td>FLE</td><td>KIT</td><td>[10, 16, 7]</td><td>Coat color (spotting)</td></tr><tr><td>18</td><td>6</td><td>72.906</td><td>72.927</td><td>7.0 10<sup>-5</sup></td><td></td><td></td><td></td><td></td></tr><tr><td>19</td><td>6</td><td>74.542</td><td>74.609</td><td>3.1 10<sup>-6</sup></td><td></td><td></td><td></td><td></td></tr><tr><td>20</td><td>6</td><td>84.327</td><td>84.348</td><td>6.3 10<sup>-5</sup></td><td></td><td></td><td></td><td></td></tr><tr><td>21</td><td>6</td><td>90.578</td><td>90.592</td><td>3.6 10<sup>-5</sup></td><td></td><td></td><td></td><td></td></tr><tr><td>22</td><td>7</td><td>22.295</td><td>22.323</td><td>3.3 10<sup>-5</sup></td><td></td><td>GNG7</td><td></td><td></td></tr><tr><td>23</td><td>7</td><td>39.013</td><td>39.321</td><td>1.0 10<sup>-6</sup></td><td>(JER)</td><td>COMMD10, ARL10, NOP16, HIGD2A, CLTB, FAF2</td><td></td><td></td></tr><tr><td>24</td><td>7</td><td>40.958</td><td>40.972</td><td>1.4 10<sup>-6</sup></td><td></td><td>ZNF154</td><td></td><td></td></tr><tr><td>25</td><td>7</td><td>41.744</td><td>41.768</td><td>1.0 10<sup>-5</sup></td><td></td><td>TRIM41</td><td></td><td></td></tr><tr><td>26</td><td>7</td><td>43.473</td><td>43.474</td><td>7.0 10<sup>-5</sup></td><td></td><td>OR cluster</td><td></td><td></td></tr><tr><td>27</td><td>7</td><td>46.029</td><td>46.040</td><td>4.7 10<sup>-5</sup></td><td>FLE</td><td>GDF9, SHROOM1, SOWAHA</td><td>[18]</td><td>Folliculogenesis</td></tr><tr><td>28</td><td>7</td><td>46.998</td><td>47.091</td><td>4.2 10<sup>-6</sup></td><td></td><td></td><td></td><td></td></tr><tr><td>29</td><td>9</td><td>41.523</td><td>41.529</td><td>6.7 10<sup>-5</sup></td><td></td><td>CEP57L1</td><td></td><td></td></tr><tr><td>30</td><td>9</td><td>44.560</td><td>44.640</td><td>4.6 10<sup>-6</sup></td><td></td><td>PRDM1</td><td></td><td></td></tr><tr><td>31</td><td>10</td><td>5.736</td><td>5.782</td><td>4.9 10<sup>-6</sup></td><td>ANG, FLE</td><td></td><td></td><td></td></tr><tr><td>32</td><td>11</td><td>16.637</td><td>16.698</td><td>1.6 10<sup>-5</sup></td><td></td><td></td><td></td><td></td></tr><tr><td>33</td><td>11</td><td>25.631</td><td>25.650</td><td>2.1 10<sup>-5</sup></td><td>HOL</td><td>THADA</td><td>[2]</td><td>rs137792035 (intronic)</td></tr><tr><td>34</td><td>13</td><td>47.491</td><td>47.503</td><td>5.7 10<sup>-5</sup></td><td></td><td>RASSF2</td><td>[26]</td><td>Bone remodelling</td></tr><tr><td>35</td><td>13</td><td>63.879</td><td>64.545</td><td>2.0 10<sup>-5</sup></td><td>ANG</td><td>ASIP</td><td></td><td></td></tr></tbody></table>

*Continued on next page*
<table><thead><tr><th>Id</th><th>Chr</th><th>begin (Mb)</th><th>end (Mb)</th><th>p-value</th><th>breed(s)</th><th>Candidate genes</th><th>References</th><th>Annotation</th></tr></thead><tbody><tr><td>36</td><td>13</td><td>67.391</td><td>67.410</td><td>1.7 10<sup>-5</sup></td><td>FLE</td><td>CTNNBL1</td><td>[17]</td><td>rs442482970, rs208147519, rs109630255 Body Mass Index (human)</td></tr><tr><td>37</td><td>14</td><td>24.937</td><td>25.070</td><td>2.1 10<sup>-7</sup></td><td>HOL, ANG</td><td>PLAG1</td><td>[19, 13]</td><td>Stature</td></tr><tr><td>38</td><td>14</td><td>25.958</td><td>25.959</td><td>7.7 10<sup>-5</sup></td><td></td><td></td><td></td><td></td></tr><tr><td>39</td><td>14</td><td>58.185</td><td>58.186</td><td>7.8 10<sup>-5</sup></td><td></td><td>EMC2, EIF3E</td><td></td><td></td></tr><tr><td>40</td><td>14</td><td>61.828</td><td>61.840</td><td>5.0 10<sup>-5</sup></td><td>JER</td><td>ZFPM2 (human)</td><td></td><td>rs381703521, rs379012140 (conserved element, intron of human gene)</td></tr><tr><td>41</td><td>16</td><td>47.039</td><td>47.040</td><td>6.1 10<sup>-5</sup></td><td>ANG</td><td>CAMTA1</td><td>[24]</td><td>QTLs for Fat thickness, Car-cass Weight</td></tr><tr><td>42</td><td>17</td><td>29.453</td><td>29.457</td><td>5.5 10<sup>-5</sup></td><td>FLE, JER</td><td></td><td>[23]</td><td>QTLs for feed efficiency traits</td></tr><tr><td>43</td><td>18</td><td>14.305</td><td>14.872</td><td>2.1 10<sup>-8</sup></td><td>ANG, HOL, FLE</td><td>MC1R</td><td></td><td>Coat color (red and black)</td></tr><tr><td>44</td><td>18</td><td>52.664</td><td>52.683</td><td>2.5 10<sup>-5</sup></td><td>JER, FLE</td><td>CEACAM20</td><td></td><td></td></tr><tr><td>45</td><td>18</td><td>54.411</td><td>54.432</td><td>6.2 10<sup>-5</sup></td><td>HOL, JER</td><td>ARHGAP35</td><td></td><td></td></tr><tr><td>46</td><td>20</td><td>23.110</td><td>23.125</td><td>4.1 10<sup>-5</sup></td><td>JER</td><td>ANKRD55</td><td>[28, 1, 15]</td><td>Region associated with autoimmune diseases in Humans</td></tr><tr><td>47</td><td>20</td><td>24.855</td><td>28.602</td><td>5.1 10<sup>-5</sup></td><td>JER</td><td>ARL15</td><td>[21]</td><td>adiponectin regulation</td></tr><tr><td>48</td><td>20</td><td>32.342</td><td>32.352</td><td>8.1 10<sup>-5</sup></td><td>ANG</td><td>140 Kb upstream of GHR</td><td></td><td>Growth, stature</td></tr><tr><td>49</td><td>20</td><td>38.945</td><td>38.946</td><td>8.6 10<sup>-5</sup></td><td>JER, FLE</td><td>5Kb upstream of PRLR promoter region</td><td></td><td></td></tr><tr><td>50</td><td>20</td><td>50.108</td><td>50.142</td><td>2.7 10<sup>-6</sup></td><td>JER</td><td></td><td></td><td></td></tr><tr><td>51</td><td>24</td><td>14.021</td><td>14.023</td><td>8.2 10<sup>-5</sup></td><td>FLE</td><td></td><td></td><td></td></tr></tbody></table>

*Continued on next page*
<table><thead><tr><th>Id</th><th>Chr</th><th>begin (Mb)</th><th>end (Mb)</th><th>p-value</th><th>breed(s)</th><th>Candidate genes</th><th>References</th><th>Annotation</th></tr></thead><tbody><tr><td>52</td><td>24</td><td>62.218</td><td>62.223</td><td>5.6 10<sup>-5</sup></td><td></td><td></td><td></td><td></td></tr><tr><td>53</td><td>27</td><td>13.125</td><td>13.166</td><td>2.0 10<sup>-5</sup></td><td>FLE</td><td>WWC2</td><td></td><td></td></tr><tr><td>54</td><td>27</td><td>33.749</td><td>33.762</td><td>1.1 10<sup>-5</sup></td><td>JER</td><td>PLEKHA2</td><td></td><td></td></tr><tr><td>55</td><td>28</td><td>3.989</td><td>3.998</td><td>7.1 10<sup>-5</sup></td><td>FLE</td><td></td><td></td><td></td></tr><tr><td>56</td><td>29</td><td>38.204</td><td>38.280</td><td>2.5 10<sup>-6</sup></td><td>FLE</td><td>PAG12</td><td>[29, 30]</td><td>Fertility</td></tr><tr><td>57</td><td>29</td><td>51.104</td><td>51.190</td><td>1.3 10<sup>-5</sup></td><td>ANG</td><td>BRSK2</td><td></td><td></td></tr></tbody></table>
<figure><img src="image_13.png" /><figcaption>Figure S1: Histogram of genomic relationship coefficients in each breed. Diagonal terms correspond to within-animal inbreeding values and extra-diagonal terms correspond to between animal co-ancestries.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>Figure S2: Neutral allele frequency spectrum estimated in each breed. The Jersey breed is shown on a different box because it has a different sample size, which implies a different range of values on the x axis.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Figure S3: Distribution of the FLK statistic and theoretical $\chi^2(3)$ distribution.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>Figure S4: Relationship between hapFLK observed quantiles and $\chi^2(42)$ theoretical quantiles</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Figure S5: Distribution of hapFLK (left) and corresponding pp-plot (right) under neutral simulations.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>Figure S6: Proportion of Fleckvieh iHS p-values (obtained from [20]) in sweep windows vs other windows in the genome. On the left panel, the sweep windows considered were those detected in Fleckvieh, independently of what happened in other breeds. On the right panel, the sweep windows considered were those detected only in Holstein.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Figure S7: Distribution of the length of selection signatures detected with hapFLK</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>Figure S8: PP plot of the FLK test applied to SNPs of the Illumina 800K chip (left) or to all SNPs of the 1000 bull genomes project (right)</figcaption></figure>
<figure><img src="image_21.png" /><figcaption>Figure S9: Distribution of *T*<sub>OR</sub> in four populations in different genome regions. (Red) selected in the breed, (Green) neutral in all breeds and (Blue) selected in another breed.</figcaption></figure>
<figure><img src="image_22.png" /><figcaption>Figure S10: Proportion of Fleckvieh CLR p-values (obtained from [20]) in sweep windows vs other windows in the genome. On the left panel, the sweep windows considered were those specific to Fleckvieh. On the right panel, the sweep windows considered were those specific to Holstein.</figcaption></figure>
<figure><img src="image_23.png" /><figcaption>Figure S11: Haplotype cluster frequencies around the hapFLK selection signature near gene ROBO1 on BTA1. Haplotypes are clustered using the Scheet and Stephens [25] model into 15 coloured template haplotypes.</figcaption></figure>
<figure><img src="image_24.png" /><figcaption>Figure S12: FLK profile and allele frequencies on chromosome 20 around the ARL15 candidate gene. Mutations in an intro of ARL15 exhibiting FLK p-values < 10⁻⁵ are highlighted in red.</figcaption></figure>
<figure><img src="image_25.png" /><figcaption>Figure S13: FLK profile in the sweep region including SLC27A6 and RXFP3. The candidate mutation reported in Table ?? is shown in red together with its allele frequency in the selected population (Jersey).</figcaption></figure>
<figure><img src="image_26.png" /><figcaption>Figure S14: Allele frequencies in the sweep region including MAGI1. For SNPs where the ancestral allele is known (in red), the frequency is that of the derived allele. For other SNPs (in black) the frequency is that of the minor allele (among all breeds).</figcaption></figure>
<figure><img src="image_27.png" /><figcaption>Figure S15: FLK profile in the sweep region including MAGI1. The candidate mutations reported in Table ?? are shown in red together with their allele frequency in the selected population (Jersey).</figcaption></figure>
<figure><img src="image_28.png" /><figcaption>Figure S16: Haplotype diversity in regions detected with hapFLK and HMM. Haplotypes are clustered using the [25] model into 15 coloured template haplotypes. Each panel represents the template frequency in each breed.</figcaption></figure>
References

[1] Alloza I, Otaegui D, de Lapuente AL, Antiguedad A, Varade J, Nunez C, Arroyo R, Urcelay E, Fernandez O, Leyva L, et al., 2012. ANKRD55 and DHCR7 are novel multiple sclerosis risk loci. *Genes Immun* 13:253-7.

[2] Boitard S, Rocha D, 2013. Detection of signatures of selective sweeps in the blonde d'aquitaine cattle breed. *Animal Genetics* 44:579-583.

[3] Bongiorni S, Mancini G, Chillemi G, Pariset L, Valentini A, 2012. Identification of a short region on chromosome 6 affecting direct calving ease in piedmontese cattle breed. *PLoS ONE* 7:e50137.

[4] Dickinson RE, Duncan WC, 2010. The SLIT-ROBO pathway: a regulator of cell function with implications for the reproductive system. *Reproduction* 139:697-704.

[5] Dickinson RE, Hryhorskyj L, Tremewan H, Hogg K, Thomson AA, McNeilly AS, Duncan WC, 2010. Involvement of the SLIT/ROBO pathway in follicle development in the fetal ovary. *Reproduction* 139:395-407.

[6] Eberlein A, Takasuga A, Setoguchi K, Pfuhl R, Flisikowski K, Fries R, Klopp N, Furbass R, Weikard R, Kuhn C, 2009. Dissection of genetic factors modulating fetal growth in cattle indicates a substantial role of the non-SMC condensin I complex, subunit G (NCAPG) gene. *Genetics* 183:951-64.

[7] Fontanesi L, Tazzoli M, Russo V, Beever J, 2010. Genetic heterogeneity at the bovine KIT gene in cattle breeds carrying different putative alleles at the spotting locus. *Anim Genet* 41:295-303.

[8] Hata S, Abe M, Suzuki H, Kitamura F, Toyama-Sorimachi N, Abe K, Sakimura K, Sorimachi H, 2010. Calpain 8/ncl-2 and calpain 9/ncl-4 constitute an active protease complex, g-calpain, involved in gastric mucosal defense. *PLoS genetics* 6:e1001040.

[9] Hayashi T, Nozaki Y, Nishizuka M, Ikawa M, Osada S, Imagawa M, 2011. Factor for adipocyte differentiation 158 gene disruption prevents the body weight gain and insulin resistance induced by a high-fat diet. *Biol Pharm Bull* 34:1257-63.

[10] Hayes BJ, Pryce J, Chamberlain AJ, Bowman PJ, Goddard ME, 2010. Genetic architecture of complex traits and accuracy of genomic prediction: coat colour, milk-fat percentage, and type in Holstein cattle as contrasting model traits. *PLoS Genet* 6:e1001139.
[11] Hrabchak C, Varmuza S, 2004. Identification of the spermatogenic zip protein spz1 as a putative protein phosphatase-1 (pp1) regulatory protein that specifically binds the pp1cγ2 splice variant in mouse testis. *Journal of Biological Chemistry* 279:37079-37086.

[12] Karim L, Takeda H, Lin L, Druet T, Arias JA, Baurain D, Cambisano N, Davis SR, Farnir F, Grisart B, et al., 2011. Variants modulating the expression of a chromosome domain encompassing plag1 influence bovine stature. *Nature genetics* 43:405-413.

[13] Karim L, Takeda H, Lin L, Druet T, Arias JA, Baurain D, Cambisano N, Davis SR, Farnir F, Grisart B, et al., 2011. Variants modulating the expression of a chromosome domain encompassing PLAG1 influence bovine stature. *Nat Genet* 43:405-13.

[14] Klein C, Bauersachs S, Ulbrich SE, Einspanier R, Meyer HH, Schmidt SE, Reichenbach HD, Vermehren M, Sinowatz F, Blum H, et al., 2006. Monozygotic twin model reveals novel embryo-induced transcriptome changes of bovine endometrium in the preattachment period. *Biol Reprod* 74:253-64.

[15] Lill CM, Schjeide BM, Graetz C, Liu T, Damotte V, Akkad DA, Blaschke P, Gerdes LA, Kroner A, Luessi F, et al., 2013. Genome-wide significant association of ANKRD55 rs6859219 and multiple sclerosis risk. *J Med Genet* 50:140-3.

[16] Liu L, Harris B, Keehan M, Zhang Y, 2009. Genome scan of pigmentation traits in Friesian-Jersey crossbred cattle. *J Genet Genomics* 36:661-6.

[17] Liu YJ, Liu XG, Wang L, Dina C, Yan H, Liu JF, Levy S, Papasian CJ, Drees BM, Hamilton JJ, et al., 2008. Genome-wide association scans identified ctnnbl1 as a novel gene for obesity. *Human Molecular Genetics* 17:1803-1813.

[18] Otsuka F, McTavish KJ, Shimasaki S, 2011. Integral role of GDF-9 and BMP-15 in ovarian function. *Mol Reprod Dev* 78:9-21.

[19] Pausch H, Flisikowski K, Jung S, Emmerling R, Edel C, Gotz KU, Fries R, 2011. Genome-wide association study identifies two major loci affecting calving ease and growth-related traits in cattle. *Genetics* 187:289-97.

[20] Qanbari S, Pausch H, Jansen S, Somel M, Strom TM, Fries R, Nielsen R, Simianer H, 2014. Classic selective sweeps revealed by massive sequencing in cattle. *PLoS Genet* 10:e1004148.
[21] Richards JB, Waterworth D, O'Rahilly S, Hivert MF, Loos RJ, Perry JR, Tanaka T, Timpson NJ, Semple RK, Soranzo N, et al., 2009. A genome-wide association study reveals variants in ARL15 that influence adiponectin levels. *PLoS Genet* 5:e1000768.

[22] Rivadeneira F, Styrkarsdottir U, Estrada K, Halldorsson BV, Hsu YH, Richards JB, Zillikens MC, Kavvoura FK, Amin N, Aulchenko YS, et al., 2009. Twenty bone-mineral-density loci identified by large-scale meta-analysis of genome-wide association studies. *Nat Genet* 41:1199-206.

[23] Rolf MM, Taylor JF, Schnabel RD, McKay SD, McClure MC, Northcutt SL, Kerley MS, Weaber RL, 2012. Genome-wide association analysis for feed efficiency in angus cattle. *Animal Genetics* 43:367-374.

[24] Saatchi M, Schnabel RD, Taylor JF, Garrick DJ, 2014. Large-effect pleiotropic or closely linked qtl segregate within and across ten us cattle breeds. *BMC Genomics* 15:442.

[25] Scheet P, Stephens M, 2006. A fast and flexible statistical model for large-scale population genotype data: applications to inferring missing genotypes and haplotypic phase. *The American Journal of Human Genetics* 78:629-644.

[26] Song H, Kim H, Lee K, Lee DH, Kim TS, Song JY, Lee D, Choi D, Ko CY, Kim HS, et al., 2012. Ablation of Rassf2 induces bone defects and subsequent haematopoietic anomalies in mice. *EMBO J* 31:1147-59.

[27] Soung do Y, Dong Y, Wang Y, Zuscik MJ, Schwarz EM, O'Keefe RJ, Drissi H, 2007. Runx3/AML2/Cbfa3 regulates early and late chondrocyte differentiation. *J Bone Miner Res* 22:1260-70.

[28] Stahl EA, Raychaudhuri S, Remmers EF, Xie G, Eyre S, Thomson BP, Li Y, Kurreeman FA, Zhernakova A, Hinks A, et al., 2010. Genome-wide association study meta-analysis identifies seven new rheumatoid arthritis risk loci. *Nat Genet* 42:508-14.

[29] Telugu BP, Walker AM, Green JA, 2009. Characterization of the bovine pregnancy-associated glycoprotein gene family-analysis of gene sequences, regulatory regions within the promoter and expression of selected genes. *BMC Genomics* 10:185.

[30] Thompson IM, Cerri RL, Kim IH, Ealy AD, Hansen PJ, Staples CR, Thatcher WW, 2012. Effects of lactation and pregnancy on metabolic and hormonal responses and expression of selected conceptus and endometrial genes of Holstein dairy cattle. *J Dairy Sci* 95:5645-56.
[31] Yoshida CA, Yamamoto H, Fujita T, Furuichi T, Ito K, Inoue K, Ya-
mana K, Zanma A, Takada K, Ito Y, et al., 2004. Runx2 and Runx3 are
essential for chondrocyte maturation, and Runx2 regulates limb growth
through induction of Indian hedgehog. *Genes Dev* 18:952-63.

[32] Zhang L, Liu J, Zhao F, Ren H, Xu L, Lu J, Zhang S, Zhang X, Wei C, Lu G, et al., 2013. Genome-wide association studies for growth and meat production traits in sheep. *PloS one* 8:e66569.

[33] Zhao C, Tian F, Yu Y, Luo J, Hu Q, Bequette BJ, Baldwin Vi RL, Liu G, Zan L, Scott Updike M, et al., 2012. Muscle transcriptomic analyses in Angus cattle with divergent tenderness. *Mol Biol Rep* 39:4185-93.

[34] Zuo B, Zhu J, Li J, Wang C, Zhao X, Cai G, Li Z, Peng J, Wang P, Shen C, et al., 2015. microrna-103a functions as a mechanosensitive microrna to inhibit bone formation through targeting runx2. *Journal of Bone and Mineral Research* 30:330-345.
<table><caption>Table 3: Private hard sweep regions including candidate mutations. Horizontal lines are used to group closely related sweep windows, which might result from the same selection event.</caption><thead><tr><th>chr</th><th>start (Mb)</th><th>end (Mb)</th><th>sel pop</th><th>nb mut</th><th>pval FLK</th><th>sel freq</th><th>genes</th></tr></thead><tbody><tr><td>6</td><td>71.440</td><td>71.560</td><td>F</td><td>2</td><td>10<sup>-5</sup></td><td>0.92</td><td>intergenic</td></tr><tr><td>7</td><td>25.430</td><td>26.000</td><td>J</td><td>8</td><td>3.10<sup>-5</sup></td><td>0.77</td><td>CHSY3, KIAA1024L, ADAMTS19</td></tr><tr><td>7</td><td>26.280</td><td>27.060</td><td>J</td><td>16</td><td>10<sup>-5</sup></td><td>0.83</td><td>SLC27A6, FBN2, SLC12A2</td></tr><tr><td>14</td><td>24.810</td><td>25.080</td><td>H,A</td><td>7</td><td>3.10<sup>-5</sup></td><td>0.04</td><td>LYN, RPS20, PLAG1, CHCHD7</td></tr><tr><td>18</td><td>14.760</td><td>14.960</td><td>F</td><td>3</td><td>10<sup>-5</sup></td><td>0.95</td><td>ENSBTAG00000039031, MOS</td></tr><tr><td>20</td><td>24.230</td><td>24.740</td><td>J</td><td>1</td><td>3.10<sup>-5</sup></td><td>0.77</td><td>TUBB6, TUBB3, DEF8, LOC532875</td></tr><tr><td>20</td><td>25.030</td><td>25.390</td><td>J</td><td>21</td><td>2.10<sup>-5</sup></td><td>0.80</td><td>DBNDD1, GAS8, SHCBP1, MC1R</td></tr><tr><td>20</td><td>26.640</td><td>27.230</td><td>J</td><td>2</td><td>3.10<sup>-5</sup></td><td>0.77</td><td>LOC530348, SNX18, COX8A, LOC783202, HSPB3</td></tr><tr><td>20</td><td>39.830</td><td>39.970</td><td>J</td><td>1</td><td>10<sup>-5</sup></td><td>0.80</td><td>intergenic</td></tr><tr><td>22</td><td>35.680</td><td>35.790</td><td>J</td><td>1</td><td>3.10<sup>-5</sup></td><td>0.77</td><td>SLC45A2, RXFP3</td></tr><tr><td>22</td><td>35.960</td><td>36.030</td><td>J</td><td>2</td><td>5.10<sup>-5</sup></td><td>0.90</td><td>intergenic</td></tr><tr><td>27</td><td>4.140</td><td>4.230</td><td>J</td><td>1</td><td>10<sup>-4</sup></td><td>0.90</td><td>MAGI1</td></tr></tbody></table>