ARTICLE

DOI: 10.1038/s41467-017-02809-1

OPEN

# Sheep genome functional annotation reveals proximal regulatory elements contributed to the evolution of modern breeds

Marina Naval-Sanchez<fnref n="1" />, Quan Nguyen<fnref n="1" />, Sean McWilliam<fnref n="1" />, Laercio R. Porto-Neto<fnref n="1" />, Ross Tellam<fnref n="1" />, Tony Vuocolo<fnref n="1" />, Antonio Reverter<fnref n="1" />, Miguel Perez-Enciso²,³, Rudiger Brauning⁴, Shannon Clarke⁴, Alan McCulloch⁴, Wahid Zamani⁵, Saeid Naderi⁶, Hamid Reza Rezaei⁷, Francois Pompanon⁶,⁸, Pierre Taberlet⁸, Kim C. Worley⁹, Richard A. Gibbs⁹, Donna M. Muzny⁹, Shalini N. Jhangiani⁹, Noelle Cockett¹⁰, Hans Daetwyler¹¹,¹² & James Kijas<fnref n="1" />

Domestication fundamentally reshaped animal morphology, physiology and behaviour, offering the opportunity to investigate the molecular processes driving evolutionary change. Here we assess sheep domestication and artificial selection by comparing genome sequence from 43 modern breeds (Ovis aries) and their Asian mouflon ancestor (O. orientalis) to identify selection sweeps. Next, we provide a comparative functional annotation of the sheep genome, validated using experimental ChIP-Seq of sheep tissue. Using these annotations, we evaluate the impact of selection and domestication on regulatory sequences and find that sweeps are significantly enriched for protein coding genes, proximal regulatory elements of genes and genome features associated with active transcription. Finally, we find individual sites displaying strong allele frequency divergence are enriched for the same regulatory features. Our data demonstrate that remodelling of gene expression is likely to have been one of the evolutionary forces that drove phenotypic diversification of this common livestock species.

<footnote n="1">¹CSIRO Agriculture and Food, 306 Carmody Road, St. Lucia 4067 QLD, Australia. ²Centre for Research in Agricultural Genomics (CRAG), Bellaterra 08193, Spain. ³ICREA, Carrer de Lluís Companys 23, Barcelona 08010, Spain. ⁴AgResearch Ltd, Invermay Agricultural Centre, Private Bag 50034, Mosgiel 9053 Otago, New Zealand. ⁵Department of Environmental Sciences, Faculty of Natural Resources and Marine Sciences, Tarbiat Modares University, Noor, Mazandaran 46414-356, Iran. ⁶Natural Resources Faculty, University of Guilan, Guilan 41335-1914, Iran. ⁷Environmental Science Department, Gorgan University of Agricultural Sciences and Natural resources, Gorgan 49138-15739, Iran. ⁸Laboratoire d'Ecologie Alpine, Universite Grenoble Alpes, Grenoble 38041, France. ⁹Human Genome Sequencing Center, Department of Human and Molecular Genetics, Baylor College of Medicine, Houston, TX 77030, USA. ¹⁰Utah State University, 1435 Old Main Hill, Logan, UT 84322-1435-1435, USA. ¹¹Department of Economic Development, Jobs, Transport and Resources, Bundoora 3083 VIC, Australia. ¹²School of Applied Systems Biology, La Trobe University, Bundoora 3083 VIC, Australia. Correspondence and requests for materials should be addressed to J.K. (email: James.Kijas@csiro.au)</footnote>
The domestication of plants and animals commenced around 10,000 years ago and precipitated enormous societal change by transitioning human kind from hunter-gatherers to agricultural settlers¹,². The impact on domesticated animals themselves has also been transformative, with radical phenotypic, morphological and behavioural changes occurring in comparison to their wild counterparts. Elucidating the molecular basis of these changes has the potential to assist the process of animal breeding and illuminate how genotypic variation influences phenotype. Analyses using population scale SNP array data³⁻⁵ and more recently whole-genome sequence data⁶⁻¹⁰ have identified regions and genes under selection in domestic species. To extend our understanding beyond key genes, multiple studies have sought to identify categories of sites or genomic features over-represented in sweep regions. These have included evolutionarily conserved elements and proximity to transcription start sites⁶,⁸,¹⁰; however, the approach has been severely hampered by the paucity of detailed functional annotation available for livestock animal species. This is starting to be addressed by the Functional Annotation of Animal Genomes (FAANG) consortium; however, human regulatory information remains by far the most abundant courtesy of the ENCODE¹¹ and Epigenetics Roadmap initiatives¹². In both, chromatin state and DNA accessibility assays have been used across a plethora of human cell lines and tissues to comprehensively map regulatory features. This has proven invaluable for investigating the link between genetic variation, gene regulation and complex human diseases¹³,¹⁴, as well as human adaptive traits¹⁵.

Exploiting the wealth of ENCODE data in non-human species is dependent on the strength of sequence and functional conservation. Encouragingly, sequence conservation of non-coding elements has been successfully used to identify functional regulatory sequences between diverse species¹⁶,¹⁷. Further, recent comparative genomic studies in mammals have shown that despite a significant divergence of Transcription Factor Binding Sites (TFBSS)¹⁸,¹⁹ and enhancers²⁰, there is a substantial core of regulatory elements which can be characterised based on sequence conservation²¹⁻²³. For example, Mouse ENCODE comparative analyses concluded that 44% of promoters and 40% of enhancers in mice that mapped to human retained functional and tissue specific activity²².

In this study, we investigate the impact of sheep domestication and subsequent artificial selection on genomic variation. First, we analyse genome sequences from a global distribution of 43 domestic sheep breeds and 17 Asiatic mouflon, its wild counterpart. Population diversity analysis identified 1420 sweep regions under positive selection in domestic sheep. Next, we use human genome regulatory information to generate a comparative functional annotation of the sheep genome, before assessing it against histone modification data we collected from sheep adipose tissue. Finally, we assess the strength of overlap between selection signatures and genomic regions implicated in regulation of gene expression. We find strong and significant overlap for multiple components of the proximal gene regulatory machinery, including promoters and chromatin states indicative of active transcription. These results provide both a high-resolution genomic view of positive selection and a first draft functional annotation of the sheep genome. Considered together, our results suggest modification to gene regulatory networks has been an important evolutionary driver of the phenotypic changes that distinguish domestic sheep from their wild ancestors.

## Results

**Genetic variation among domestic and wild sheep.** We sequenced the genome of 67 domestic sheep drawn from 43

phenotypically diverse breeds for comparison with 17 Asiatic Mouflon representing their wild ancestor (Fig. 1a, Supplementary Data 1). Sequencing to an average depth of 11.82× coverage was followed by alignment and variant calling to identify 28.1 million high-quality SNPs. Analysis using the protein coding gene set of the reference assembly OARv3.1²⁴ revealed <1% of SNP (218,762) were located in exons (Supplementary Table 1), and of these the majority (133,891) were synonymous variants (Supplementary Table 2). Despite sampling 43 geographically dispersed breeds, domestic sheep contained fewer private SNP (Fig. 1b, c, Supplementary Table 3, Supplementary Data 2), and lower nucleotide diversity ($\pi$) than wild sheep (domestic sheep $\pi$=0.16% per nucleotide, mouflon $\pi$=0.20%; Fig. 1d). Reduced diversity is an expected consequence of domestication, however in sheep the magnitude of change appears smaller than for some other domesticates⁸,²⁵. This suggests sheep domestication captured a wide sampling of diversity from Asiatic mouflon, consistent with the finding that many modern breeds retain high-effective population size⁴ and a diverse representation of mtDNA haplogroups²⁶. Three metrics that assess genetic relatedness revealed Mouflon and domestic breeds remain relatively closely related. We found 70% of SNP in domestic sheep were also polymorphic in Mouflon (14.3 out of 20.4 M SNP, Fig. 1b) and the allele frequencies of SNP were highly correlated (Fig. 1c). Finally, population differentiation between domestic and wild sheep was low for an inter-species comparison ($F_{ST}$ = 0.093). We next evaluated the relationship between individual genomes using Principal Component Analysis (PCA), identifying two distinct clusters of Mouflon from geographically separate areas of their range (Supplementary Fig. 1). PCA of domestic sheep genomes clustered animals according to their geographic origin of breed formation (Fig. 1e). PC1 captured an east-west cline, PC2 a north-south cline and breeds from the Middle East (the domestication centre) took a central position with near zero PC values, recapitulating a phylo-geographic structure consistent with earlier SNP array based investigations⁴,²⁷.

**Genomic regions under selection in domestic sheep.** We sought to identify genomic regions impacted by positive selection during the domestication of sheep and/or the artificial selection that followed. By sequencing a globally representative collection of phenotypically diverse breeds and searching for common sweep regions, we expect to enrich for selection events acting during domestication or prior to the global radiation of breeds. However caution is required, as sweep events common to divergent breeds may be the result of much more recent selection²⁸. Domestic sheep genomes were interrogated for shared regions characterised by both low diversity and high divergence compared to wild sheep. First, we estimated allele frequency divergence as $F_{ST}$ between species in 20 kb windows (Supplementary Fig. 2a). To assess the direction of selection effect, nucleotide diversity ($\pi$) was estimated separately within domestic and wild sheep, before the ratio of mouflon/domestic $\pi$ was plotted for each genomic bin⁹,²⁹ (Supplementary Fig. 2b). Joint evaluation of both metrics defined 1420 bins putatively under selection in domestic animals and a further 28 in Mouflon (P-adj <0.01, Fig. 2a, Supplementary Data 3–4). Physically co-located bins (within 50 kb) were joined to define 635 regions (average size 36.3 kb) containing 430 genes (Supplementary Data 5–6). Inspection revealed many genes have previously been shown to be influenced by selection in livestock or associated in sheep biology (Table 1, Supplementary Table 4). These include SOCS2 for animal weight and milk production³⁰, IITCH-ASIP for coat colour³¹ and VEGFA involved in early development and reproductive success⁴,³². Additional genes implicated in other domesticates, include SOX2 in rabbits⁸,
<figure><img src="image_1.png" /><figcaption>**Fig. 1** Genome diversity and relatedness of wild and domestic sheep. **a** The geographic distribution of 43 breeds sampled for whole-genome sequencing. Population based proportion of allele type is given for each major geographic group of breeds, after calling reference variants present in the Texel derived genome assembly OARv3.1. **b** Proportion and number of private and shared SNP for the collection of wild and domestic sheep genomes. **c** Reference allele frequency (RAF) correlation between domestic and mouflon populations. RAF domestic and mouflon represent the frequency of the reference allele in domestic and mouflon populations respectively. (Pearson correlation: 76.36% *p*-value <2.2E-16). **d** Distribution of nucleotide diversity (π) within species, estimated within 20 kb bins. The boxplot compares the bin mean and variance. The distribution means, evaluated using the Wilcoxon rank-sum test, are significantly different (*p*-value <2.2E-16). **e** Principal component (PC) analysis clustering individual domestic sheep coloured to reflect the geographic origin of breed development</figcaption></figure>

NR6A1 associated with vertebrae number in domestic pig<sup>7</sup>, TBX3 associated with pigmentation in horses<sup>33</sup> and ESR1 associated with prolificacy in pigs<sup>34</sup> (Supplementary Fig. 3a). One prominent example is *KITLG*, a protein involved in melanocyte development and pigmentation in a variety of mammalian species<sup>35-39</sup>. In sheep, selection sweeps flank either side of the *KITLG* coding exons suggesting they flag alleles with a gene regulatory function (Fig. 2b, c). Another example is the *NCAPG-LCORL* locus implicated in controlling height and stature, however it is unclear which gene exerts the effect<sup>7,40,41</sup>. Domestic sheep contain a selective sweep spanning the *LCORL* promoter region, while mouflon contain a sweep downstream of *NCAPG* (Supplementary Fig. 3b). Many of these candidate genes are likely to contribute to phenotypic variation, however in the absence of individually recorded trait data analysed using an approach, such as GWAS, establishing the direct link between these genes and

their functional consequence is difficult using selection sweep methodology alone.

**Biological processes affected by domestication and selection.**

To identify the biological processes influenced by domestication and artificial selection, we mapped the ovine sweep regions to the human genome and examined their ontological enrichment using GREAT<sup>42</sup> (see Methods section and Supplementary Data 7). This revealed a significant enrichment (Padj <0.05) hypergeometric and binomial test) for processes consistent with the evolution of domestic animals from their wild ancestors. One prominent example is dorsal/ventral patterning (Table 2, Supplementary Data 7) as reduction in body size was a quick and early response to domestication across species<sup>43</sup>. A second is regulation of lipid metabolism which may reflect human meditated alteration of
<figure><img src="image_2.png" /><figcaption>**Fig. 2** Genomic regions putatively under positive selection in sheep. **a** Population differentiation ($F_{ST}$) and relative nucleotide diversity between wild and domestic sheep ($\pi_{mouflon}/\pi_{domestic}$) was estimated in 20 kb genomic bins. A total of 1420 outlier bins exhibiting evidence for selection in domestic sheep genomes are indicated in red (corresponding to Z-test $P < 0.001$, where $F_{ST} > 0.156$ and In ratio $>0.672$). **b** Genome-wide distribution of relative nucleotide diversity. Positive values identify genomic bins with depressed diversity in domestic sheep compared with mouflon, consistent with positive selection sweeps involved with domestication and selection (Table 1, Supplementary Table 4). **c** Selective sweeps located either side of *KITLG* coding exons. Integrated Genome Visualisation (IGV) screenshot of chromosome 3:124236035–125400001 illustrates the reduction of SNP variation in 67 domestic genomes compared with 17 Mouflon. The two regions identified by our selection metrics are shown inside the dashed red boxes</figcaption></figure>
<table><caption>Table 1 Genes associated with selective sweeps in domestic animals reported in other studies</caption><thead><tr><th>Chr</th><th>Start</th><th>End</th><th>Avg rnk F<sub>ST</sub> and π</th><th>π</th><th>F<sub>ST</sub></th><th>Closest genes</th><th>Distance to closest gene</th><th>Function</th></tr></thead><tbody><tr><td>13</td><td>50310001</td><td>50650000</td><td>2</td><td>3.56</td><td>0.53</td><td>PANK2</td><td>0</td><td>Neurodegeneration</td></tr><tr><td>3</td><td>124790001</td><td>125080000</td><td>21.5</td><td>2.80</td><td>0.35</td><td>KITLG</td><td>27570</td><td>Coat colour</td></tr><tr><td>2</td><td>184990001</td><td>185040000</td><td>42.5</td><td>1.98</td><td>0.33</td><td>GLI2</td><td>139099</td><td>Growth</td></tr><tr><td>1</td><td>199640001</td><td>199790000</td><td>51</td><td>1.79</td><td>0.35</td><td>IGF2BP2</td><td>0</td><td>Adiposity</td></tr><tr><td>4</td><td>78890001</td><td>79000000</td><td>51.5</td><td>2.44</td><td>0.31</td><td>GLI3</td><td>29289</td><td>Pigmentation</td></tr><tr><td>17</td><td>59450001</td><td>59520000</td><td>62.5</td><td>2.32</td><td>0.30</td><td>TBX3</td><td>215541</td><td>Pigmentation</td></tr><tr><td>9</td><td>36140001</td><td>36240000</td><td>63</td><td>1.72</td><td>0.34</td><td>PLAG1</td><td>0</td><td>Fertility, stature</td></tr><tr><td>15</td><td>21900001</td><td>22010000</td><td>70.5</td><td>1.60</td><td>0.39</td><td>BCO2</td><td>0</td><td>Yellow-fat</td></tr><tr><td>9</td><td>13570001</td><td>13610000</td><td>82.5</td><td>1.85</td><td>0.30</td><td>DGAT1</td><td>0</td><td>Milk fat</td></tr><tr><td>6</td><td>37420001</td><td>37510000</td><td>144</td><td>1.56</td><td>0.29</td><td>LCORL</td><td>0</td><td>Weight/height</td></tr><tr><td>3</td><td>124640001</td><td>124680000</td><td>184</td><td>1.69</td><td>0.26</td><td>KITLG</td><td>0</td><td>Coat colour</td></tr><tr><td>8</td><td>75670001</td><td>75700000</td><td>228.5</td><td>1.97</td><td>0.24</td><td>ESR1</td><td>0</td><td>Litter size; prolifacy</td></tr><tr><td>6</td><td>38500001</td><td>38530000</td><td>245.5</td><td>1.62</td><td>0.25</td><td>LCORL</td><td>1047670</td><td>Weight/height</td></tr><tr><td>5</td><td>19460001</td><td>19480000</td><td>312.5</td><td>1.23</td><td>0.27</td><td>IRF1</td><td>5040</td><td>Immune function</td></tr><tr><td>3</td><td>129730001</td><td>129750000</td><td>359</td><td>0.94</td><td>0.31</td><td>SOC52</td><td>7494</td><td>Weight and milk production</td></tr><tr><td>1</td><td>100990001</td><td>101010000</td><td>383.5</td><td>1.26</td><td>0.24</td><td>TCHH</td><td>0</td><td>Hair</td></tr><tr><td>4</td><td>78800001</td><td>78820000</td><td>386.5</td><td>1.15</td><td>0.26</td><td>GLI3</td><td>209289</td><td>Pigmentation</td></tr><tr><td>5</td><td>85880001</td><td>85900000</td><td>405</td><td>1.11</td><td>0.26</td><td>MEF2C</td><td>173852</td><td>Skeletal muscle development</td></tr><tr><td>2</td><td>105830001</td><td>105850000</td><td>417.5</td><td>1.27</td><td>0.23</td><td>HAND2</td><td>235173</td><td>Limb development</td></tr><tr><td>20</td><td>17530001</td><td>17550000</td><td>434.5</td><td>1.207</td><td>0.24</td><td>VEGFA</td><td>147890</td><td>Reproduction</td></tr><tr><td>13</td><td>63380001</td><td>63400000</td><td>469</td><td>1.237</td><td>0.22</td><td>ITCH</td><td>0</td><td>Coat colour</td></tr><tr><td>1</td><td>203460001</td><td>203480000</td><td>481</td><td>0.977</td><td>0.25</td><td>SOX2</td><td>52282</td><td>Stem cell maintenance</td></tr><tr><td>6</td><td>35290001</td><td>35310000</td><td>486.5</td><td>1.05</td><td>0.24</td><td>GPRIN3</td><td>201294</td><td>Brain development</td></tr><tr><td>6</td><td>70290001</td><td>70310000</td><td>492.5</td><td>1.24</td><td>0.21</td><td>KIT</td><td>55390</td><td>Coat colour</td></tr></tbody></table>

<table><caption>Table 2 Gene ontology biological process—top 15 enrichment results for identified selective sweeps (Great v.3.0)⁴²</caption><caption>GO biological process</caption><thead><tr><th>Term name</th><th>Binom raw P-value</th><th>Binom FDR Q-Val</th><th>Binom fold enrichment</th><th>Binom observed region hits</th><th>Binom region set coverage</th></tr></thead><tbody><tr><td>Primary alcohol catabolic process</td><td>1.12E-27</td><td>2.34E-24</td><td>23.4</td><td>27</td><td>0.020</td></tr><tr><td>Regionalisation</td><td>3.22E-19</td><td>9.62E-17</td><td>2.3</td><td>135</td><td>0.103</td></tr><tr><td>Regulation of MAPK cascade</td><td>5.84E-19</td><td>1.60E-16</td><td>2.1</td><td>173</td><td>0.132</td></tr><tr><td>Development of primary male sexual characteristics</td><td>6.82E-19</td><td>1.83E-16</td><td>3.1</td><td>83</td><td>0.063</td></tr><tr><td>Gland development</td><td>8.90E-19</td><td>2.21E-16</td><td>2.3</td><td>138</td><td>0.105</td></tr><tr><td>Male gonad development</td><td>2.34E-18</td><td>5.20E-16</td><td>3.3</td><td>73</td><td>0.055</td></tr><tr><td>Male sex differentiation</td><td>4.05E-18</td><td>7.98E-16</td><td>3.0</td><td>83</td><td>0.063</td></tr><tr><td>Dorsal/ventral pattern formation</td><td>9.33E-18</td><td>1.68E-15</td><td>3.5</td><td>67</td><td>0.051</td></tr><tr><td>Regulation of lipid metabolic process</td><td>2.58E-17</td><td>4.09E-15</td><td>2.7</td><td>92</td><td>0.070</td></tr><tr><td>Negative regulation of cell cycle</td><td>6.13E-17</td><td>9.28E-15</td><td>2.6</td><td>97</td><td>0.074</td></tr><tr><td>Embryonic hemopoiesis</td><td>1.30E-16</td><td>1.80E-14</td><td>6.0</td><td>35</td><td>0.027</td></tr><tr><td>Odontogenesis</td><td>8.84E-16</td><td>1.04E-13</td><td>3.2</td><td>66</td><td>0.050</td></tr><tr><td>Branching morphogenesis of an epithelial tube</td><td>1.43E-15</td><td>1.61E-13</td><td>2.63</td><td>91</td><td>0.069</td></tr><tr><td>Negative regulation of cellular protein metabolic process</td><td>2.00E-14</td><td>1.65E-12</td><td>2.1</td><td>127</td><td>0.097</td></tr><tr><td>Development of primary sexual characteristics</td><td>3.83E-14</td><td>2.96E-12</td><td>2.2</td><td>105</td><td>0.079</td></tr></tbody></table>

muscularity and fatness as sheep were actively managed as a food source⁴⁴. A third example involves sexual maturation, as humans increasingly exerted control over breeding, pressure to maintain beneficial sexual fitness traits in the wild were removed, reducing sexual dimorphism and altering the timing of reproduction⁴⁵. The residual consequences of these changes were identified in the data, as we identified significant biological processes associated with sexual differentiation, particularly in male sexual characteristics and differentiation (Table 2, Supplementary Data 7). We next evaluated sheep sweep regions against known mouse mutants and their phenotypic consequences. This identified strong enrichment for traits, including oogenesis, litter size, sexual maturation and ovary morphology (Supplementary Data 8).

**Preliminary functional annotation of the sheep genome.** In the absence of detailed functional annotation for any livestock genome, we used comparative genomics to predict ovine regulatory elements using human data¹¹,¹². Reciprocal liftOver was used to predict the ovine genome location of ENCODE promoters and enhancers, along with 12 chromatin states built using 127 diverse epigenomes¹² (Fig. 3a, see Methods section). ENCODE histone modification marks and transcription factor binding sites for proximal promoters and distal enhancers were predicted with similar mapping efficiency (Supplementary Table 5) while higher variability was observed across the set of 12 chromatin states. Approximately 40–50% of active state elements, or elements associated with expressed genes, were reciprocally mapped
<figure><img src="image_3.png" /><figcaption>**Fig. 3** Prediction and validation of sheep gene regulatory elements. **a** Human Epigenome Roadmap chromatin states and ENCODE features used in a reciprocal liftOver protocol to predict regulatory regions in sheep. **b** Overlap between chromosomal locations identified by experimental sheep tissue ChIP-Seq (H3K4me3 and H3K27ac) and regions predicted by liftOver. H3K4me3 peaks are compared against predicted Roadmap chromatin state TssA identifying active promoters (left) and unique H3K27ac peaks are compared to the Roadmap Enh state, excluding those regions overlapping with predicted TssAs, identifying distal enhancers (right). **c** The proportion of predicted promoters (left) and enhancers (right) recovered by ChIP-Seq using either H3K4me3 peaks (blue bar) or unique H3K27ac peaks (orange bar). This demonstrates the predicted annotations retain a level of specificity between components of the gene regulatory apparatus. **d** An example annotated region containing *LCORL* highlights the correspondence between ChIP-Seq data (Sheep H3K27ac) and liftOver predictions (TSSA, TSSAFlnk, ENCODE proximal, Enh, ENCODE distal)</figcaption></figure>
including active TSS promoter states (TssA 41%; TssAFlank 48.5%), and transcription at the 5' and 3' ends of genes (TxFlnk 40%). Mapping efficiency broadly reflected evolutionary conservation for each genomic feature as measured by enrichment for mammalian conserved elements¹².

**Human epigenetic data captures sheep regulatory elements.** To evaluate if our comparative genomics approach successfully predicted regulatory elements in sheep, we performed H3K4me3 and H3K27ac ChIP-seq analysis of sheep late gestation perirenal brown adipose tissue. H3K4me3 identifies active gene promoters⁴⁶,⁴⁷, and we identified a total of 16,098 regions. Comparing their genomic location against Roadmap active promoter state TssA predictions revealed an 85% recovery rate (Fig. 3b, 13,767/16,098). A large number of predicted promoters were not present, possibly reflecting our ChIP-Seq data originate from only one sheep tissue and timepoint. H3K27ac identifies both active promoters and active enhancers, and 35,366 regions were identified. Of these, 26,492 regions were uniquely identified by H3K27ac but not H3K4me3 and defined in this study as enhancers. Comparing the location of these 26,492 experimentally defined enhancers with predicted enhancer state Enh, excluding those overlapping the defined state TssA, revealed a recovery rate of 71% (Fig. 3b, 18,850/26,492). The recovery rate and significance for each Roadmap and ENCODE predicted genomic feature is given in Supplementary Fig. 4–7 and Supplementary Table 6. This revealed the predicted features were significantly enriched for sheep regulatory elements and retained specificity between different components of the gene regulatory machinery. For example, H3K4me3 marks recovered 85% of predicted promoters but only 29% of enhancers while H3K27ac unique marks recovered only 19% of predicted promoters but 71% of enhancers (Fig. 3c). Next, we performed H3K27me3 ChIP-Seq in sheep brown adipose tissue, a mark related to Polycomb repression. This yielded 31,942 elements. Comparison against the set of predicted regulatory features revealed strong enrichment only for repressive features, including repressed Polycomb and bivalent/poised TSSs (Supplementary Table 6). This provides confidence our comparative approach successfully recovered biologically meaningful genome feature annotations with the sensitivity to discriminate promoters, enhancers, and repressive regulatory elements. An example of a predicted functional annotation is given in Fig. 3d. Finally, we sought to explore how robust the reciprocal liftOver approach is to the combination of species in which it is applied. Similar results were obtained when using human Roadmap elements to predict promoters and enhancers in the mouse, cow and pig genomes, where in each case available histone modification mark data sets were strongly enriched for predicted elements (Supplementary Fig. 8, Supplementary Table 7).

Three additional properties of the predicted genome features were assessed to evaluate their quality. First, the distance from each feature to the nearest transcription start site (TSS) was assessed for both the sheep and human genome. This revealed that ENCODE predicted proximal promoter elements were consistently closer to TSS in the sheep genome than distal enhancers (Supplementary Fig. 9). Second, levels of nucleotide diversity within each predicted annotation feature were compared. This showed epigenomic modifications reporting proximal promoters contained lower diversity than those marking distal enhancers, consistent with findings for other species⁴⁸ (Supplementary Fig. 10). Third, nucleotide diversity was compared between elements active in either a ubiquitous or tissue specific manner. Predicted promoters active in many Roadmap epigenomes (>100 cell types; TssA $\pi < 0.06$; Enh $\pi < 0.15$) had higher evolutionary constraint compared with promoters expressed in a restricted

manner (<10 cell types, TssA $\pi < 0.09$; Enh $\pi < 0.12$, Supplementary Fig. 11). Taken together, the results from our sheep specific histone modification data and the analysis of nucleotide diversity validate that the collection of predicted annotated features are of sufficient utility to address evolutionary questions.

**Selection sweeps are enriched for proximal features.** A key objective of the study sought to investigate if ovine genomic regions putatively affected during domestication and artificial selection displayed an enrichment of functional domains associated with protein translation, transcription or the regulation of gene expression. As input we used the collection of genome features described above that included: (i) sheep H3K4me3 and H3K27ac ChIP-Seq peaks; (ii) Epigenome Roadmap predicted chromatin states; (iii) ENCODE predicted proximal and distal elements, and; (iv) gene model annotations from the reference assembly OARv3.1. To test for enrichment, a location overlap approach developed for interpretation of epigenomic data (LOLA)⁴⁹ was used to assess each genomic feature against the distribution of sweep regions (1420 bins, Fig. 2a). This revealed selection sweeps strongly co-localised with gene regulatory elements after correction for multiple testing (Fig. 4a). Three components of the proximal gene regulatory machinery were significant: enhancer elements physically associated with genes (EnhG, LogP = 12.83), actively transcribed states (Tx, LogP = 10.68) and proximal promoter states (TxFlnk, LogP = 8.03); (Fig. 4a, Supplementary Table 8–9). Next, we repeated the analysis after linking neighbouring outlier bins (within 50 kb) to define 635 selective sweep regions. This was a more permissive analysis as the genome fraction implicated with domestication and artificial selection was higher (Methods section). Sweep regions were significantly enriched for 14 genome features (Supplementary Fig. 12, Supplementary Table 10). Strikingly, the eight most strongly enriched were each components of the proximal gene regulatory machinery and active transcription including sheep ChIP-Seq modifications, Epigenome Roadmap states (TxFlnk, EnhG, Tx, TssA), proximal ENCODE features obtained by DNase I signal and transcription factor binding sites and translated exons (CDS, Supplementary Fig. 12).

**Site frequency analysis.** Domestication and selection are anticipated to have altered the allele frequency of loci mechanistically involved in the phenotypic and molecular changes that distinguish domesticates from their wild ancestors. We therefore applied site frequency analysis to identify nucleotides with divergent allele frequency when compared between wild and domestic sheep genomes. This exploits the availability of whole-genome sequence by assessing every site, and is independent from the selection sweep methodology used above. We first estimated allele frequencies within domestic and wild sheep using 14 M sites found to be segregating in both species. Comparison between species revealed the majority displayed highly correlated allele frequencies, with 8.9 M sites (63% of 14 M, Supplementary Data 9) having allele frequency difference ($\Delta$AF) <0.2. This is consistent with the behaviour of neutral loci between closely related species. Only a very small fraction of loci exhibited strongly divergent allele frequency. These are characteristic of sites under selection, however they can be difficult to distinguish from those impacted by genetic drift alone. A total of 9161 SNP (0.016%) had $\Delta$AF>0.8 and only 1059 sites had divergence approaching the maximum (>0.9, Supplementary Data 10). We expect this set of loci are likely to have played an outsized role in domestication and selection, prompting us to characterise them in more detail. To evaluate the impact of gene inactivation we searched these sites for nonsense mutations introducing stop
<figure><img src="image_4.png" /></figure>

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

<figure><img src="image_6.png" /><figcaption>**Fig. 4** Genomic feature enrichment in selection sweeps and differentiated sites. **a** Strength of enrichment for 29 genome features within 1420 sweep bins assessed by location overlap⁴⁹. Genome features were derived from four different sources. The significance threshold from multiple testing is represented by the horizontal line. **b** Intersection of delta allele frequency (ΔAF) with protein coding gene annotations from reference OARv3.1. The number of SNP in ΔAF bins is given at left, and the *M*-value (at right) was calculated by comparing the frequency of SNP in each genome feature and ΔAF bin with the corresponding frequency across all bins. **c** As for **b** using chromatin state annotations derived from the Roadmap dataset. Additional *M*-value results using predicted ENCODE marks are provided in Supplementary Fig. 13</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>**Fig. 5** Candidate causal missense mutation in *FBXL3*. a *FBXL3* gene representation in human hg19 coordinates chr13:77579389-77601337. b Multiple sequence alignment of SNP across representative vertebrates. c Allele frequency of the G (reference) or A (alternative) allele in domestic sheep (*O. aries*) (*n* = 67) or mouflon (*O. orientalis*) (*n* = 17). d *FBXL3* amino acid sequence with the substitution at residue 182 highlighted in red. Residue 358 (green) is associated with the mouse After hour (Afh) phenotype residue 364 with the overtime phenotype (Ovt)</figcaption></figure>

codons or that modify the splicing machinery. This revealed a single variant, located in a splice donor site of *MICAL3*, was present in sites with extreme ΔAF (>0.8) (Supplementary Data 9). Therefore, as in chicken⁶, rabbit⁸ and pigs⁷, our conclusion is that gene inactivation via nonsense (frame-changing) and splice site mutation has not played a major role during sheep domestication. We then performed enrichment analysis to evaluate the number of SNP across ΔAF bins compared to their location within annotated genome features. This approach has been used previously⁸,⁵⁰ and captures data from every segregating SNP available. Assessment of protein coding gene annotations derived from OARv3.1 revealed non-synonymous substitutions were strongly enriched for SNP in high-ΔAF bins, while no relationship was detected for introns, UTR or intergenic regions (Fig. 4b). Regulatory elements are embedded within the intergenic

genome feature; however, they likely represent an insufficiently large proportion to cause intergenic regions to display enrichment. Our result do demonstrate that selection within coding regions has clearly played a role in adaptive change during domestication and selection. Of 27 SNP with ΔAF>0.80, there is only one missense mutation where domestic sheep carry an allele different from other vertebrate genomes (Fig. 5). The variant is in *FBXL3* (F-Box And Leucine-Rich Repeat Protein 3) which is known to be involved in regulation of circadian rhythms⁵¹,⁵² and mice with missense mutations have an altered circadian period⁵³,⁵⁴. Many modern sheep breeds have elongated reproductive seasonality compared to wild sheep, and seasonality is photoperiod dependent. It is therefore possible the *FBXL3* mutation identified here plays a role in controlling seasonality in sheep.
In addition to analysis of protein coding gene annotation, the same site frequency approach was applied to the RoadMap and ENCODE predicted features (Fig. 4c, Supplementary Fig. 13). Striking enrichment was detected for chromatin states associated with both active and developmentally poised promoters (TxFlnk, TssBiv), as well as proximal enhancers (EnhG, Fig. 4c). Further, analysis of ENCODE derived features show enrichment for proximal elements that was absent for each distal category of regulatory element (Supplementary Fig. 13). We artificially altered exonic allele frequencies to search for any impact on other annotation features, and found modest enrichment that serves as a note of caution (Supplementary Fig. 14). None the less, the findings are remarkably similar to the feature enrichments present in the selective sweep analysis (Fig. 4a), thus providing evidence that domestication and artificial selection have prefer-entially impacted the proximal gene regulatory machinery.

Discussion

We explore the genomic consequence of domestication and subsequent artificial selection by sequencing 43 breeds of domestic sheep and comparing them to descendants of their wild ancestors. The data revealed both protein coding genes and proximal regulatory elements played major roles during the evolution of modern sheep, addressing a long standing question in evolutionary biology concerning the relative contribution played by regulatory versus protein coding sequence on pheno-typic variation⁵⁵,⁵⁶.

Genome sequencing clearly showed protein coding genes have been directly altered during the development of modern sheep. Genomic regions with evidence of positive selection were sig-nificantly enriched for coding exons, and missense (non-synon-yymous) substitutions were fourfold over represented in sites with extreme allele frequency differences between wild and domestic genomes. Sweep regions contained a collection of genes pre-viously implicated in the control of pigmentation, reproduction and stature and biological process enrichment was observed for sexual differentiation and altered timing of reproduction.

A more difficult task has been any detailed evaluation of the
contribution made by sequence elements that regulate gene
expression. To address this we built a comparative functional
annotation of sheep and validated its utility by comparison with
experimentally derived histone modification data identifying
active sheep promoters, enhancers and repressors. Our major
finding is multiple chromatin states and regulatory elements are
strongly enriched in both selective sweeps and sites with diver-
gent allele frequency. Involvement of the gene regulatory
machinery is consistent with the emerging view that complex trait
architecture is often governed by changes in gene expression. For
example, the majority of trait associated SNP in human disease
GWAS fall outside transcribed regions⁵⁷. Intersecting these
associated SNP with epigenomic marks reveals strong enrichment
within both ENCODE⁵⁷ and Roadmap genome features¹².
Similarly in livestock, testing production trait associated SNP has
revealed enrichment within non-coding elements neighbouring
genes, suggesting a role for variants that influence gene expres-sion⁵⁸. Next, evo-devo studies have put forward the relevance of
regulatory mutations affecting morphology, behaviour and
adaptation⁵⁹,⁶⁰. Seminal examples include *Drosophila* wing pig-
ment patterns⁶¹ and skeletal reduction in stickleback fish⁶².
Perhaps of most relevance to animal domestication, analysis of
modern rabbits suggests a strong enrichment has occurred at
conserved non-coding sites expected to be enriched for regulatory
function⁸. In addition, the earliest stages of trout domestication
are characterised by large, heritable changes to gene expression⁶³.
Our analysis, armed with the extensive epigenomic annotation

described here, finds evidence that remodelled gene regulation
has been a force during the profound morphological, behavioural
and metabolic changes occurring through the transition to a
domesticated state.

The richness of predicted epigenomic feature annotations allowed us to discriminate between the components of the gene regulatory apparatus that have been preferentially impacted. Perhaps surprisingly, active promoters and proximal enhancers were strongly enriched in both sweep regions and sites with divergent ΔAF, in contrast to each distal ENCODE mark or inactive Roadmap chromatin states, such as repressed Polycomb and constitutive heterochromatin. This suggests remodelling of control elements immediately adjacent to genes has been a dominant mechanism of change, however it is currently unknown why distal elements appear to have made a diminished contribution.

Most experiments into the consequences of animal domes-
tication and selection have relied on SNP array data and have
focussed on protein coding genes in the absence of functional
annotation³⁻⁵. Genome sequencing has dramatically improved
the resolution of selective sweep regions in this study of sheep and
other studies of pigs and chickens⁶,⁶⁴,⁶⁵. Further, our comparative
application of human epigenomic data has proven highly
insightful, however limitations remain. The reciprocal BLAST
methodology may have introduced bias, resulting in an over-
representation of evolutionarily conserved promoter elements
compared with enhancers. Indeed the recovery of enhancers
within sheep ChIP-Seq data was lower than for promoters,
however the difference was not large and enhancers have a higher
evolutionary turnover compared to proximal regulatory elements
or protein-coding sequences⁴⁸,⁶²,⁶⁶. Further, it is possible that
species specific and tissue specific regulatory elements are likely to
be under represented as a consequence of both our methodology
and the application of ChIP-Seq using only a single tissue. A
separate note of caution relates to the sampling of animals
sequenced in the study, which is likely to have an impact on the
final collection of sweeps identified as the choice of genomes used
has a measurable impact on allele frequencies (Supplementary
Fig. 15). A final challenge arises from linkage disequilibrium and
the overlapping nature of some genome annotation features,
which together blur the relative contribution of physically co-
located genomic features (Supplementary Fig. 14). A case in point
is the observation that synonymous sites displayed enrichment in
our site frequency analysis. It is not clear if this reflects a func-
tional role for many synonymous positions via cryptic effects on
transcription or mRNA transport as previously suggested⁶⁷ or if
their proximity to missense mutations has generated the enrich-
ment via linkage disequilibrium. Despite these limitations, the
collection of annotated genome features represents a useful
resource for the research community to interpret GWAS and
inform the design of SNP panels to increase power to deliver
genomic prediction. Knowledge concerning the elements likely to
regulate gene expression enhances our understanding of the sheep
genome, and suggests that remodelling gene expression has been
a key mechanism during the evolution of this important livestock
species.

Methods

Samples. A total of 70 animals were sampled from 43 domestic breeds and sub-jected to genome sequencing (Supplementary Data 1). These comprise 46 animals selected from an earlier SNP array based global survey of breed diversity⁴ and another six animals used for SNP discovery, construction of the SNP50 BeadChip and CNV detection. The final group of 18 individuals have not been examined before. Breeds were drawn from Asia (12), Africa (6), the Middle East (13), the Americas (8), the United Kingdom (8) and continental Europe (23). Whole-genome sequence data for 19 Asian mouflon (*Ovis orientalis*) was collected and made available by the NEXTGEN project (http://nextgen.epfl.ch/). Fastq files were
downloaded from the ENA public repository (http://www.ebi.ac.uk/ena/data/view/PRJEB3139) and processed as described below for the domestic sheep genomes.

**Genome sequencing, variant detection and annotation.** Paired-end short insert libraries were constructed using 5 µg of genomic DNA and sequenced on the Illumina HiSeq 2000 platform. Reads were mapped against the sheep reference assembly v3.1<sup>24</sup> using BWA aligner v0.7.12 (bwa aln + bwa sampe, default parameters)<sup>68</sup>. Animals were sequenced to an average median depth of 11.8 × (8.4–17.2×) (Supplementary Data 1). Duplicate reads were removed using Picard tools (http://broadinstitute.github.io/picard/), and local realignment around INDELS was performed using GATK v3.2.0<sup>69</sup>. Variant detection and SNP diversity analyses were performed using SAMTOOLS 1.2.1 mpileup and annotated using VCFTools v0.1.14<sup>70</sup>. After obtaining genotype calls for a total of 89 samples the following filters were applied using a combination of VCFtools<sup>70</sup> and in-house scripts: (i) SNP were retained in positions with read depth between 5× and twice the average depth per sample; (ii) minimum mapping quality of 30 and base quality of 20 were applied; (iii) SNP within 5 bp of INDELS were removed; (iv) for SNP pairs separated by <4 bp, the lower quality variant was excluded; (v) tri-allelic variants were removed; (vi) SNP called in <90% of animals were excluded and (vii) SNP displaying an excess of heterozygosity were excluded (--hwe 0.001). This defined a set of 28,100,631 SNP across domestic (67) and mouflon (17) genomes. A total of five low-coverage animals were excluded (3 domestic and 2 mouflon). PLINK v1.9 was used to perform genetic diversity estimates and PCA (https://www.cog-genomics.org/plink2)<sup>71</sup>. The variant effect predictor tool from ensembl (version 78) was used to identify 24 separate SNP classifications, including coding, missense and non-synonymous substitutions, intron and intergenic, in relation to the gene models annotated on reference assembly OARv3.1<sup>24</sup> (Supplementary Tables 1–2).

**Genome scanning for selection sweeps.** Population differentiation between wild and domestic populations were measured as $F_{ST}$ using the Weir and Cockerham method<sup>72</sup>. Average SNP $F_{ST}$ values were plotted in 20 kb genomic bins with a 10 kb step (Fig. 2a; Supplementary Fig. 2). Nucleotide diversity ($\pi$) was estimated for the same bins in both wild and domestic sheep genome collections as the number of heterozygous SNP by bin size. Genomic bins with fewer than 20 SNP were excluded. Positive selection is characterised by reduced $\pi$, prompting us to estimate a ratio (wild/ domestic) that identifies the direction of selection<sup>9,29</sup>. The ratio (wild/ domestic) is given in Fig. 2a and plotted in genomic order (Supplementary Fig. 2). The genome-wide threshold to declare significance of Z-transformed values ($F_{ST}$ > 0.156 and log $\pi$ ratio >0.672) represents Bonferroni adjusted $p$-value of <0.01.

**Biological process enrichment within selective sweep regions.** Locus-based Gene Ontology enrichment was performed using GREAT v3.0.0<sup>42</sup>. Genomic bins with evidence of positive selection were translated to human coordinates (GRC37/hg19) using UCSC's liftOver tool (minMatch = 0.1)<sup>73</sup>. The membership and frequency of gene regulatory domains present was then compared against a background set representing all genomic bins using two statistical tests. A binomial and a hypergeometric test assessed the enrichment of molecular function terms and biological process terms. We also used GREAT to interrogate the Mouse Genome Informatics (MGI) database phenotypes<sup>74</sup>.

**Building a comparative sheep functional annotation.** Our approach exploited the wealth of functional annotation data generated by the Epigenome Roadmap and ENCODE studies<sup>11,12</sup>. We performed reciprocal liftOver (minMatch = 0.1)<sup>73</sup>, meaning elements that mapped to sheep also needed to map in the reverse direction back to human with high quality. This bidirectional comparative mapping approach was applied to 12 chromatin states defined using 5 core histone modification marks, H3K4me3, H3K4me1, H3K36me3, H3K9me3 and H3K27me3. Mapping success is given in Supplementary Table 5. The same approach was applied to ENCODE marks derived from 94 cell types (https://www.encodeproject.org/data/annotations/v2/) with DNase-seq and TF ChIP-seq.

**ChIP-seq of ovine adipose tissue.** Chromatin immunoprecipitation followed by next generation sequencing (ChIP-Seq) of the histone chromatin modification H3K4me3, H3K27ac and H3K27me3 was undertaken using late gestation fetal perirenal adipose tissue (PRAT). This tissue was selected as it may be important in a production setting, where it influences lamb survival under some circumstances. All procedures involving animals were carried out with approval from the University of Adelaide Animal Ethics Committee. The adipose tissue was taken at 130 days post conception from three animals for H3K4me3 and two animals for H3K27ac and H3K27me3. Input nucleosomal DNA from the combined fetal samples were used as input controls. Nuclei were isolated from 1.2 g of frozen PRAT by pulverising the tissue under liquid nitrogen at 4 °C in Nuclei Buffer (0.3 M sucrose, 60 mM KCl, 15 mM NaCl, 5 mM MgCl₂, 0.1 mM EGTA, 15 mM Tris-HCl, pH 7.5) with the inclusion of protease inhibitors (0.5 mM DTT, 0.1 mM AEBSF, CompleteTM EDTA-free protease inhibitor cocktail (Roche Diagnostics GmbH, Mannheim Germany) using a hand held Dounce homogeniser. The filtered homogenate was centrifuged (3000×g, 5 min) and the pellet recovered and resuspended in ice-cold Nuclei Buffer containing 0.4% IGEPAL CA-630 (Sigma),

incubated on ice for 5 min and the suspension was then layered on a cushion of Nuclei Buffer containing 1.2 M sucrose. The sample was centrifuged at 10,000×g (25 min, 4 °C), the chromatin pellet recovered and resuspended in micrococcal nuclease digestion buffer. The isolated chromatin was treated with micrococcal nuclease (MNase) (0.2 U/µl) (New England BioLabs, Ipswich, MA, USA) and soluble chromatin recovered using mild sonication and centrifugation. Specific antibodies to H3K4me3 (Abcam ab8580) and H3K27ac (Abcam ab4729) and H3K27me3 (Millipore 07-449) (10 µg per incubation) were used for immunoprecipitation (4 °C, gentle mixing, 15 h). Immune complexes were isolated using Protein A-Sepharose (4 h, 4 °C) and washed sequentially with low salt buffers supplemented with 0.1% NP-40 and a final wash in 1xTE buffer. The chromatin was released from the Protein A Sepharose using 1% SDS, 0.1 M NaHCO₃, treated with proteinase K, extracted with TE-saturated phenol/chloroform (1:1) and DNA purified using Minielute DNA columns (Qiagen, Maryland, USA). DNA (20–100 ng) was quantified using Quant-ITTM PicoGreen and validated for purity and size using a High Sensitivity DNA chip (Agilent 2100 Bioanalyser; Agilent Technologies, Santa Clara, CA, USA). DNA samples (10 ng) isolated from the immunopurified nucleosomes and input nucleosomes were used to generate sequencing libraries for each sample using the Illumina TruSeq ChIP-Seq kit and these were sequenced using the Illumina HiSeq2000 (H3K27ac and H3K4me3) platform (Illumina, San Diego, CA, USA) and Illumina Genome Analyser II (GAII) for H3K27me3. Approximately 205 M 50 bp sequence reads per sample were produced. Reads were mapped to the unmasked ovine genome sequence (Ovis aries Oar_v3.1.74) using the NGS core tool mapping application in CLCBIO (mapping parameters: length fraction = 0.7; similarity fraction = 0.8; penalties, mismatch = 2, insertion = 3, deletion = 3). A range of 75–90% of reads was uniquely mapped and retained. Peak calling comparing the H3K4me3 or H3K27ac ChIP-Seq versus the input control was performed using MACS<sup>75</sup>. Only peaks found in both replicates per chromatin mark, either H3K4me3, H3K27ac or H3K27me3, were further considered. For H3K4me3 there were 16098 peaks with a FDR 0.1%, whereas for H3K27ac there were 35622 peaks with 1% FDR and for H3K27me3 31942 peaks with 5% FDR. The number of peaks is comparable to those previously reported for the same chromatin marks in cattle<sup>76</sup>.

**Validation of the predicted sheep epigenome using ChIP-seq.** First, we calculated the overlap between experimentally observed (ChIP-Seq) sheep enhancer and promoter elements with our collection of predicted epigenome features. We created 1000 randomisations for each genomic feature using bedtools shuffle (–noOverlapping)<sup>77</sup>. We then calculated overlap for each randomisation against the experimentally observed (ChIP-Seq) sheep enhancer and promoter elements. Finally, we counted how many times an equal or greater overlap observed in the original features were observed in the 1000 randomisations to estimate an empirical $p$-value per overlap and feature.

**Epigenome feature enrichment within selective sweep regions.** The Bioconductor package Locus Overlap Analysis (LOLA) was used to assess the relationship between selection sweeps and genome regions containing functional annotation<sup>49</sup>. This process used (i) a 'query set' comprising each genome feature derived from four sources as detailed in Fig. 3a; (ii) a 'reference set' of 1420 genomic bins with evidence of selection (Fig. 2a) and (iii) a 'Universe Set' containing all 20 kb genome-wide bins. A Fisher's exact test with false discovery rate correction was performed to assess the significance of overlap in each pairwise comparison<sup>49</sup>.

**Analysis of divergent allele frequencies.** Allele frequency (AF) was estimated for each SNP separately for domestic and wild sheep genomes using PLINK V1.9 (--freq –within)<sup>71</sup>. Next, the difference ($\Delta$AF) was used to identify SNP with divergent AF between species. SNP was allocated into $\Delta$AF bins and compared across functional annotation features by using bedtools intersect<sup>77</sup>. Log2 fold change of the observed SNP count for each genomic feature in each bin was compared against the expected SNP count ($M$-value; average across bins)<sup>8,50</sup>. Statistical significance of the deviations from the expected values was tested with a chi-squared ($\chi^2$) test.

**Data availability.** Genome sequence for each of 68 domestic sheep are available through the International Sheep Genomics Consortium (http://www.sheepmapp.org/) with Fastq files deposited in the SRA repository (https://www.ncbi.nlm.nih.gov/bioproject/160933) and SheepGenomes DB project (http://sheepgenomesdb.org/) via European Variant Archive (EVA) study PRJEB14685 (http://www.ebi.ac.uk/ena/data/view/PRJEB14685). Genome sequence for 19 mouflon (Ovis orientalis) were made available by the NEXTGEN project (http://nextgen.epfl.ch/) with Fastq files deposited in the ENA public repository as study number PRJEB3139 and secondary accession ERP001583 (https://www.ebi.ac.uk/ena/data/view/PRJEB3139). Fetal Perirenal adipose tissue H3K27ac, H3K4me3 ChIP-Seq data are accessible in the GEO series (GSE90812). Human epigenomic tracks converted to sheep coordinates are available through CSIRO Data Access Portal https://doi.org/10.4225/08/5a03a9c39a0ba.
Received: 2 February 2017 Accepted: 3 December 2017

## References

1. Diamond, J. Evolution, consequences and future of plant and animal domestication. *Nature* **418**, 700–707 (2002).

2. Larson, G. et al. Current perspectives and the future of domestication studies. *Proc. Natl Acad. Sci. USA* **111**, 6139–6146 (2014).

3. Ramey, H. R. et al. Detection of selective sweeps in cattle using genome-wide SNP data. *BMC Genom.* **14**, 382 (2013).

4. Kijas, J. W. et al. Genome-wide analysis of the world's sheep breeds reveals high levels of historic mixture and strong recent selection. *PLoS Biol.* **10**, e1001258 (2012).

5. Elferink, M. G. et al. Signatures of selection in the genomes of commercial and non-commercial chicken breeds. *PLoS ONE* **7**, e32720 (2012).

6. Rubin, C.-J. et al. Whole-genome resequencing reveals loci under selection during chicken domestication. *Nature* **464**, 587–591 (2010).

7. Rubin, C.-J. et al. Strong signatures of selection in the domestic pig genome. *Proc. Natl Acad. Sci. USA* **109**, 19529–19536 (2012).

8. Carneiro, M. et al. Rabbit genome analysis reveals a polygenic basis for phenotypic change during domestication. *Science* **345**, 1074–1079 (2014).

9. Qiu, Q. et al. Yak whole-genome resequencing reveals domestication signatures and prehistoric population expansions. *Nat. Commun.* **6**, 10283 (2015).

10. Schubert, M. et al. Prehistoric genomes reveal the genetic foundation and cost of horse domestication. *Proc. Natl Acad. Sci. USA* **111**, E5661 (2014).

11. ENCODE Project Consortium. An integrated encyclopedia of DNA elements in the human genome. *Nature* **489**, 57–74 (2012).

12. Roadmap Epigenomics Consortium. et al. Integrative analysis of 111 reference human epigenomes. *Nature* **518**, 317–330 (2015).

13. GTEx Consortium. et al. Genetic effects on gene expression across human tissues. *Nature* **550**, 204–213 (2017).

14. Ongen, H. et al. Estimating the causal tissues for complex traits and diseases. *Nat. Genet.* **49**, 1676–1683 (2017).

15. Fraser, H. B. Gene expression drives local adaptation in humans. *Genome Res.* **23**, 1089–1096 (2013).

16. Siepel, A. et al. Evolutionarily conserved elements in vertebrate, insect, worm, and yeast genomes. *Genome Res.* **15**, 1034 (2005).

17. Pennacchio, L. A. et al. In vivo enhancer analysis of human conserved non-coding sequences. *Nature* **444**, 499–502 (2006).

18. Schmidt, D. et al. Five vertebrate ChIP-seq reveals the evolutionary dynamics of transcription factor binding. *Science* **328**, 1036 (2010).

19. Stefflova, K. et al. Cooperativity and rapid evolution of cobound transcription factors in closely related mammals. *Cell* **154**, 530 (2013).

20. Villar, D., Flicek, P. & Odom, D. T. Evolution of transcription factor binding in metazoans—mechanisms and functional implications. *Nat. Rev. Genet.* **15**, 221–233 (2014).

21. Cheng, Y. et al. Principles of regulatory information conservation between mouse and human. *Nature* **515**, 371–375 (2014).

22. Yue, F. et al. A comparative encyclopedia of DNA elements in the mouse genome. *Nature* **515**, 355 (2014).

23. Vierstra, J. et al. Mouse regulatory DNA landscapes reveal global principles of cis-regulatory evolution. *Science* **346**, 1007–1012 (2014).

24. Jiang, Y. et al. The sheep genome illuminates biology of the rumen and lipid metabolism. *Science* **344**, 1168 (2014).

25. Groenen, M. A. M. et al. Analyses of pig genomes provide insight into porcine demography and evolution. *Nature* **491**, 393–398 (2012).

26. Meadows, J. R. S., Cemal, I., Karaca, O., Gootwine, E. & Kijas, J. W. Five ovine mitochondrial lineages identified from sheep breeds of the near East. *Genetics* **175**, 1371–1379 (2007).

27. Ciani, E. et al. Merino and Merino-derived sheep breeds: a genome-wide intercontinental study. *Genet. Sel. Evol. GSE* **47**, 64 (2015).

28. Girdland Flink, L. et al. Establishing the validity of domestication genes using DNA from ancient chickens. *Proc. Natl Acad. Sci. USA* **111**, 6184–6189 (2014).

29. Li, M. et al. Genomic analyses identify distinct patterns of selection in domesticated pigs and Tibetan wild boars. *Nat. Genet.* **45**, 1431–1438 (2013).

30. Rupp, R. et al. A point mutation in suppressor of cytokine signalling 2 (Socs2) increases the susceptibility to inflammation of the mammary gland while associated with higher body weight and size and higher milk production in a sheep model. *PLoS Genet.* **11**, e1005629 (2015).

31. Norris, B. J. & Whan, V. A. A gene duplication affecting expression of the ovine ASIP gene is responsible for white and black sheep. *Genome Res.* **18**, 1282–1293 (2008).

32. Brown, H. M., Fabre Nys, C., Cognié, J. & Scaramuzzi, R. J. Short oestrous cycles in sheep during anoestrus involve defects in progesterone biosynthesis and luteal neovascularisation. *Reproduction* **147**, 357–367 (2014).

33. Imsland, F. et al. Regulatory mutations in TBX3 disrupt asymmetric hair pigmentation that underlies Dun camouflage color in horses. *Nat. Genet.* **48**, 152–158 (2016).

34. Wang, C. et al. Genome-wide analysis reveals artificial selection on coat colour and reproductive traits in Chinese domestic pigs. *Mol. Ecol. Resour.* **15**, 414–424 (2015).

35. Wilkinson, S. et al. Signatures of diversifying selection in European pig breeds. *PLoS Genet.* **9**, e1003453 (2013).

36. Reissmann, M. & Ludwig, A. Pleiotropic effects of coat colour-associated mutations in humans, mice and other mammals. *Semin. Cell Dev. Biol.* **24**, 576–586 (2013).

37. Sulem, P. et al. Genetic determinants of hair, eye and skin pigmentation in Europeans. *Nat. Genet.* **39**, 1443–1452 (2007).

38. Mengel-From, J., Wong, T. H., Morling, N., Rees, J. L. & Jackson, I. J. Genetic determinants of hair and eye colours in the Scottish and Danish populations. *BMC Genet.* **10**, 88 (2009).

39. Guenther, C. A., Tasic, B., Luo, L., Bedell, M. A. & Kingsley, D. M. A molecular basis for classic blond hair color in Europeans. *Nat. Genet.* **46**, 748–752 (2014).

40. Vaysse, A. et al. Identification of genomic regions associated with phenotypic variation between dog breeds using selection mapping. *PLoS Genet.* **7**, e1002316 (2011).

41. Widmann, P. et al. A systems biology approach using metabolomic data reveals genes and pathways interacting to modulate divergent growth in cattle. *BMC Genom.* **14**, 798 (2013).

42. McLean, C. Y. et al. GREAT improves functional interpretation of cis-regulatory regions. *Nat. Biotechnol.* **28**, 495–501 (2010).

43. Zeder, M. A., Emshwiller, E., Smith, B. D. & Bradley, D. G. Documenting domestication: the intersection of genetics and archaeology. *Trends Genet.* **22**, 139–155 (2006).

44. Chessa, B. et al. Revealing the history of sheep domestication using retrovirus integrations. *Science* **324**, 532–536 (2009).

45. McPherson, F. J. & Chenoweth, P. J. Mammalian sexual dimorphism. *Anim. Reprod. Sci.* **131**, 109–122 (2012).

46. Cain, C. E., Blekhman, R., Marioni, J. C. & Gilad, Y. Gene expression differences among primates are associated with changes in a histone epigenetic modification. *Genetics* **187**, 1225–1234 (2011).

47. Santos-Rosa, H. et al. Active genes are tri-methylated at K4 of histone H3. *Nature* **419**, 407–411 (2002).

48. Yu, F. et al. Population genomic analysis of 962 whole genome sequences of humans reveals natural selection in non-coding regions. *PLoS ONE* **10**, e0121644 (2015).

49. Sheffield, N. C. & Bock, C. LOLA: enrichment analysis for genomic region sets and regulatory elements in R and Bioconductor. *Bioinformatics* **32**, 587 (2016).

50. Barrio, A. M. et al. The genetic basis for ecological adaptation of the Atlantic herring revealed by genome sequencing. *eLife* **5**, e12081 (2016).

51. Busino, L. et al. SCFFbxl3 controls the oscillation of the circadian clock by directing the degradation of cryptochrome proteins. *Science* **316**, 900–904 (2007).

52. Xing, W. et al. SCF(FBXL3) ubiquitin ligase targets cryptochromes at their cofactor pocket. *Nature* **496**, 64–68 (2013).

53. Siepka, S. M. et al. Circadian mutant Overtime reveals F-box protein FBXL3 regulation of cryptochrome and period gene expression. *Cell* **129**, 1011–1023 (2007).

54. Godinho, S. I. H. et al. The after-hours mutant reveals a role for Fbxl3 in determining mammalian circadian period. *Science* **316**, 897–900 (2007).

55. Hoekstra, H. E. & Coyne, J. A. The locus of evolution: evo devo and the genetics of adaptation. *Evolution* **61**, 995–1016 (2007).

56. Halligan, D. L. et al. Contributions of Protein-Coding and Regulatory Change to Adaptive Molecular Evolution in Murid Rodents. *PLoS Genet.* **9**, e1003995 (2013).

57. Schaub, M. A., Boyle, A. P., Kundaje, A., Batzoglou, S. & Snyder, M. Linking disease associations with regulatory information in the human genome. *Genome Res.* **22**, 1748–1759 (2012).

58. Koufariotis, L., Chen, Y.-P. P., Bolormaa, S. & Hayes, B. J. Regulatory and coding genome regions are enriched for trait associated variants in dairy and beef cattle. *BMC Genom.* **15**, 436 (2014).

59. Wray, G. A. The evolutionary significance of cis-regulatory mutations. *Nat. Rev. Genet.* **8**, 206–216 (2007).

60. Carroll, S. B. Evo-devo and an expanding evolutionary synthesis: a genetic theory of morphological evolution. *Cell* **134**, 25–36 (2008).

61. Gompel, N., Prud'homme, B., Wittkopp, P. J., Kassner, V. A. & Carroll, S. B. Chance caught on the wing: cis-regulatory evolution and the origin of pigment patterns in Drosophila. *Nature* **433**, 481–487 (2005).
62. Chan, Y. F. et al. Adaptive evolution of pelvic reduction in sticklebacks by recurrent deletion of a Pitx1 enhancer. *Science* **327**, 302–305 (2010).

63. Christie, M. R., Marine, M. L., Fox, S. E., French, R. A. & Blouin, M. S. A single generation of domestication heritably alters the expression of hundreds of genes. *Nat. Commun.* **7**, 10676 (2016).

64. Bosse, M. et al. Genomic analysis reveals selection for Asian genes in European pigs following human-mediated introgression. *Nat. Commun.* **5**, 4392 (2014).

65. Ai, H. et al. Adaptation and possible ancient interspecies introgression in pigs identified by whole-genome sequencing. *Nat. Genet.* **47**, 217–225 (2015).

66. Arbiza, L. et al. Genome-wide inference of natural selection on human transcription factor binding sites. *Nat. Genet.* **45**, 723–729 (2013).

67. Chamary, J. V., Parmley, J. L. & Hurst, L. D. Hearing silence: non-neutral evolution at synonymous sites in mammals. *Nat. Rev. Genet.* **7**, 98–108 (2006).

68. Li, H. & Durbin, R. Fast and accurate short read alignment with Burrows-Wheeler transform. *Bioinformatics* **25**, 1754 (2009).

69. DePristo, M. A. et al. A framework for variation discovery and genotyping using next-generation DNA sequencing data. *Nat. Genet.* **43**, 491–498 (2011).

70. Danecek, P. et al. The variant call format and VCFtools. *Bioinformatics* **27**, 2156–2158 (2011).

71. Chang, C. C. et al. Second-generation PLINK: rising to the challenge of larger and richer datasets. *GigaScience* **4**, 7 (2015).

72. Weir, B. S. & Cockerham, C. C. Estimating F-statistics for the analysis of population structure. *Evolution* **38**, 1358–1370 (1984).

73. Hinrichs, A. S. et al. The UCSC genome browser database: update 2006. *Nucleic Acids Res.* **34**, D590 (2006).

74. Blake, J. A. et al. Mouse genome database (MGD)-2017: community knowledge resource for the laboratory mouse. *Nucleic Acids Res.* **45**, D723–D729 (2017).

75. Zhang, Y. et al. Model-based analysis of ChIP-Seq (MACS). *Genome Biol.* **9**, R137 (2008).

76. Villar, D. et al. Enhancer evolution across 20 mammalian species. *Cell* **160**, 554–566 (2015).

77. Quinlan, A. R. BEDTools: the Swiss-army tool for genome feature analysis. *Curr. Protoc. Bioinforma* **47**, 11.12.1 (2014).

## Acknowledgements

M.N.-S. and Q.N. are funded by the CSIRO Science Excellence Research Office. M.P.-E. was funded by a CSIRO McMaster Visiting Fellowship. Mouflon genome sequences were generated by the NextGen Consortium, as part of the European Union's Seventh Framework Programme (FP7/2010-2014 grant 244356). We thank the Baylor College of Medicine Human Genome Sequencing Center production teams, including S.L. Lee, S. Dugan, S. Jhangiani, D. Bandaranaike, M. Batterton, M. Bellair, C. Bess, K. Blankenburg, H. Chao, S. Denson, H. Dinh, S. Elkadiri, Q. Fu, B. Hernandez, D. Kalra,

M. Javaid, J.C. Jayaseelan, S. Lee, M. Li, X. Liu, T. Matskevitch, M. Munidasa, R. Najjar, L. Nguyen, F. Ongeri, N. Osuji, L. Perales, L.-L. Pu, M. Puazo, S. Qi, J. Quiroz, R. Raj, J. Shafer, H. Shen, N. Tabassum, L.-Y. Tang, A. Taylor, G. Weissenberger, K. Wilczek-Boney, Y.-Q. Wu, Y. Xin, Y. Zhang, Y. Zhu and X. Zou.

## Author contributions

M.N.-S. and J.K. conceived the experiments, analysed the data and wrote the paper. Q.N., S.W., L.R.P.-N., A.R., M.P.-E, R.B, S.C, A.C. and H.D. performed data analysis. R.T. and T.V. contributed sheep histone modification data and analysis. K.C.W., R.A.G., D.M.M. and S.N.J. contributed sequencing of domestic sheep genomes. W.Z., S.N., H.R.R., F.P. and P.T. contributed genome sequence from Asian mouflon. J.K., N.C., K.C.W., R.B. and H.D. provided project management under the auspices of the International Sheep Genomics Consortium.

## Additional information

Supplementary Information accompanies this paper at https://doi.org/10.1038/s41467-017-02809-1 and Supplementary Data 9 can be found at https://doi.org/10.4225/08/5a78ff454afe2

**Competing interests:** The authors declare no competing financial interests.

**Reprints and permission** information is available online at http://npg.nature.com/reprintsandpermissions/

**Publisher's note:** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

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

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2018