Whole-genome resequencing of wild and domestic sheep identifies genes associated with morphological and agronomic traits

Understanding the genetic changes underlying phenotypic variation in sheep (Ovis aries) may facilitate our efforts towards further improvement. Here, we report the deep resequencing of 248 sheep including the wild ancestor (O. orientalis), landraces, and improved breeds. We explored the sheep variome and selection signatures. We detected genomic regions harboring genes associated with distinct morphological and agronomic traits, which may be past and potential future targets of domestication, breeding, and selection. Furthermore, we found non-synonymous mutations in a set of plausible candidate genes and significant differences in their allele frequency distributions across breeds. We identified PDGFD as a likely causal gene for fat deposition in the tails of sheep through transcriptome, RT-PCR, qPCR, and Western blot analyses. Our results provide insights into the demographic history of sheep and a valuable genomic resource for future genetic studies and improved genome-assisted breeding of sheep and other domestic animals.

S heep (Ovis aries) is an important livestock species, which has provided meat, wool, skin, and milk for humans since the Neolithic. Characterization of genome-wide sequence variation and identification of phenotype-associated functional variants are essential steps for guidance of genome-assisted breeding in the near future. The impact of domestication and subsequent selection on genomic variation has recently been investigated in sheep 1,2 , and a number of quantitative trait loci (QTLs) and functional genes have been associated with phenotypic traits 3 . However, most of these investigations focused on a few phenotypes and involved a limited number of molecular markers and breeds/populations. So far, whole-genome resequencing has allowed the identification of genomic variants involved in domestication and genetic improvement for several domestic plants (e.g., rice and soybean) 4,5 and animals (e.g., cattle and sheep) 1,2,6 .
The completion of a sheep reference genome 7 has allowed comparison of the genomes from a wide collection of phenotypically diverse authentic landraces and improved breeds of domestic sheep with their wild ancestors. In this study, we resequence the genomes of 16 Asiatic mouflon, 172 sheep from 36 landraces and 60 sheep from six improved breeds to a depth of 25.7× coverage. Tests for selective sweeps and genome-wide association studies (GWAS) identify a number of selected regions and genes potentially affected by domestication and associated with several important morphological and agronomic traits. Moreover, we conduct a survey of non-silent single-nucleotide polymorphisms (SNPs) and gene-containing copy number variations (CNVs), which are part of the selective signatures. These data provide a valuable genomic resource for facilitating future molecular-guided breeding and genetic improvement of domestic sheep, potentially valuable in the face of ongoing climate change and consequent impacts in agricultural practice. In addition, our findings contribute to further understanding of the demographic history of sheep and the molecular basis of distinct phenotypes in the species and other animals. . Using 10,007 homozygous reference loci on the Ovine HD BeadChip for 211 individuals, false-positive SNP calling rates of 6.38% and 5.37% were observed for GATK and SAMtools, respectively. After filtering, the false-positive rate for the SNP set identified by both methods was estimated to be 0.66%. Moreover, inspection of 68 randomly selected SNPs in candidate genes from 1,414 individuals of 21 breeds obtained by Sanger sequencing produced an overall validation rate of 95.69% (Supplementary Table 1  Patterns of variation. The 28.36 million SNPs were analyzed across the three groups of sheep (Asiatic mouflon, landraces, and improved breeds). A majority up to 23.27 million SNPs were observed in Asiatic mouflon at 7.77-9.16 million per individual (12.06 million to be unique for this group), followed by 14.38 million in landraces at 5.62-8.92 million per individual (1.06 million to be unique) and 14.01 million in improved breeds at 5.90-6.90 million per individual (1.08 million to be unique) (Fig. 2, Table 1, Supplementary Table 2, and Supplementary Data 3). Using the Asiatic mouflon reference genome (ftp://ftp. ebi.ac.uk/pub/databases/nextgen/ovis/assembly/mouflon.Oori1. PRJEB3141/) for SNP calling, we identified 28.75 million SNPs in Asiatic mouflon, which was higher than that based on the sheep reference genome Oar v.4.0 (23.27 million).

Results
We observed 12.09 million SNPs shared between landraces and improved breeds, which exceeded that shared between the mouflon and landraces (10.38 million) or between the mouflon and improved breeds (9.98 million) ( Fig. 2 and Table 1). Pairwise genome-wide F ST values also indicated that genomic differentiation between landraces and improved breeds (F ST = 0.032, P = 0.041) was less than that between Asiatic mouflon and landraces (F ST = 0.125, P = 0.015) or between Asiatic mouflon and improved breeds (F ST = 0.132, P = 0.016). The genomic diversities (π) in Asiatic mouflon, landraces, and improved breeds based on the SNPs with <10% missing data were 0.00127, 0.00113, and 0.00109, respectively. The distribution of SNPs at various regions near or within genes was similar in Asiatic mouflon, landraces, and improved breeds (Supplementary Table 1 Summary information of whole-genome variations  identified in Asiatic mouflon, landraces, and improved  breeds.   Variations  Asiatic mouflon  Landraces  Improved breeds   SNPs  23,269,423  14,382,975  14,008,509  Indels  3,501,571  3,481,234  2,743,640  Insertions  1,351, Table 2). The ratio of non-synonymous and synonymous substitutions in wild (0.55) and domestic sheep (0.53-0.54) was comparable.
To ascertain the effect of major evolutionary transitions (e.g., domestication and intensive artificial breeding) on CNVs and SVs 8 , these variants were pooled for Asiatic mouflon, landraces, and improved breeds separately ( Fig. 2 and Table 1). This yielded a high depth of coverage and information about the uniqueness and sharing of CNVs and SVs among the three groups ( Fig. 2 and Table 1). The abundance of CNVs per individual ranged from 443 to 804 (average of 563) for Asiatic mouflon, from 311 to 777 (average of 589) for landraces, and from 514 to 686 (average of 616) for improved breeds (Supplementary Data 4). The number of SVs per individual varied between 5,035 and 6,657 (mean = 5,874) in Asiatic mouflon, between 4,515 and 6,323 (mean = 5,393) in landraces, and between 4,863 and 6,203 (mean = 5,366) in improved breeds (Supplementary Data 5). In contrast to the abundance of SNPs identified in the wild species, we detected significant differences in the numbers of CNVs (Kruskal-Wallis, P = 0.047) and SVs (Kruskal-Wallis, P < 0.001) per individual among the three groups (Supplementary Data 4 and 5).
From a total of 6,929 common CNV regions (read-depth signal value <0.3 or >1.7 for all 248 individuals; see Online Methods), we found 946 functional genes overlapping with 1,999 CNV regions (Supplementary Data 9 and Supplementary Notes). The top 15 significant Gene Ontology (GO) terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways for the 6,220 unique CNV regions in the 232 domestic sheep were enriched for biological processes involved in binding of sperm to zona pellucida and cell-cell recognition as well as for pathways associated with neural system function and immune system response (Supplementary Data 10). GO and KEGG pathway analyses for the 402 unique CNV regions in the 16 Asiatic mouflon uncovered enriched GO terms associated with adhesion that play essential roles in cell shape, motility, and proliferation as well as pathways involved in metabolism, neural system, and focal adhesion (Supplementary Data 10).
Population structure, linkage disequilibrium, and demography. To understand the population structure and demographic history of Old World domestic sheep, we utilized the SNP dataset (Fig. 1a) in a number of contexts and analyses as follows. Using the Asiatic mouflon as an outgroup, we produced a phylogenetic tree that divided domestic populations into four subgroups of European, Middle Eastern, Asian, and two African lineages (i.e., the Dorper sheep and the 10 African landraces)  Neighbor-joining (NJ) tree of the 248 individuals constructed using the p-distances between individuals, with Asiatic mouflon as an outgroup. c Plots of principal components 1 and 2 for the 248 individuals. d Decay of linkage disequilibrium in the Asiatic mouflon, landraces, and improved breeds. e Neighbor-joining tree of five genetic groups based on the Reynolds genetic distances. Red numbers beside divergence nodes are bootstrap values based on 1,000 replications. A scale bar represents branch length in terms of percent divergences (%). Source data are provided as a Source Data file.
NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-16485-1 ARTICLE ( Fig. 1b and Supplementary Fig. 3). This geographic subdivision was confirmed by principal component analysis (PCA; Fig. 1c) and clustering analysis based on maximum likelihood estimation ( Supplementary Fig. 4). A neighbor-joining (NJ) tree constructed for the above four groups further revealed a close genetic affinity between East Asian and Middle Eastern populations, whereas the two African lineages showed larger genetic divergence from the other three subgroups (Fig. 1e). Nucleotide diversities (π) in European, Asian, African, and Middle Eastern populations based on the SNPs with <10% missing data were 0.00105, 0.00118, 0.00113, and 0.00114, respectively ( Supplementary Fig. 5a). Asiatic mouflon were more genetically similar to Middle Eastern sheep than to other domestic populations as measured by pairwise F ST (Supplementary Fig. 5b).
Linkage disequilibrium (LD, measured as r 2 ) decreased to half of its maximum value at 2.8 kb in Asiatic mouflon but at 12.  Table 3) showed that European populations had a higher level of LD (25.1 kb) than that in Asian populations (9.8 kb). Pairwise sequentially Markovian coalescent (PSMC) analysis revealed concordant demographic trajectories for wild and domestic sheep, with two expansions and two contractions in N e during the last one million years ( Supplementary Fig. 7). The estimated N e for Asiatic mouflon and domestic populations~50 generations ago 9 were 344.1 and 73.7-199.3, respectively, being inversely correlated with the extent of LD ( Supplementary Fig. 5c) as expected. These observations suggested that artificial selection and genetic isolation, leading to the formation of breeds, had stronger effects on LD and N e than on nucleotide diversity.
Genomic signatures related to domestication. To identify genomic regions influenced by domestication, we compared the genomes of 16 Asiatic mouflon and five old landrace populations representing different geographic and genetic origins: 5 Dutch Drenthe Heathen 10 , 10 East-Asian Hu 11 , 10 Central-Asian Altay 11 , one African Djallonké 12 , and one Middle Eastern Karakul sheep 13 . Using the cross-population composite likelihood ratio (XP-CLR) test, we scanned for genomic regions with extreme allele frequency differentiation. The top 1% XP-CLR values identified 302 putative selective sweeps in the five old landraces after annotation and removing repeats ( Fig. 3a and Supplementary Data 11). As genomic regions targeted by artificial selection may be expected to have decreased levels of genetic variation, we also measured and plotted nucleotide diversity (π) along their genomes. Selecting the windows with the top 1% diversity ratios, i.e., low diversity in the five old landraces but high in the Asiatic mouflon, we found 529 putative selective sweeps (Supplementary Data 12). Combination of the XP-CLR and π ratio analyses unveiled 144 putative selective regions covering or being near to 261 genes in the five old landraces (Supplementary Data 13). Additional analyses involving the integrated haplotype score (iHS) analysis (top 5% outliers) and the Hudson-Kreitman-Aguadé (HKA) test (χ 2 = 5.99, df = 2, P = 0.05) identified 899 and 1,503 putative selective sweeps, respectively (Supplementary Data 14 and 15). Sixty-five and 71 selected genes identified by both XP-CLR and π ratio analyses were also detected by the iHS and HKA analyses, respectively (Supplementary Data 16 and 17). A comparison of the domestication-associated selective sweeps and known QTLs 14 (permutation test, P < 0.001; Supplementary were also identified to be the targets of selection in the comparison of Asiatic mouflon with domestic sheep in two recent studies (Supplementary  Table 6 and Supplementary Data 21) and found that the variant allele frequencies of non-synonymous SNPs in five genes (PDE6B, BCO2, ADAMTSL3, NKX2-1, and an olfactory receptor 51A4-like gene LOC101108252) and the genotype pattern in LOC101108252 showed significant differences (Mann-Whitney, P < 0.01) between the Asiatic mouflon and the five old landraces (Fig. 3a, b and Supplementary Table 7). In the functional enrichment analysis of the 261 genes putatively influenced by domestication, we identified the top 15 over-represented GO terms and 12 KEGG pathways (Supplementary Data 22). Specifically, four biological process GO terms and one KEGG pathway were associated with biosynthesis. Five biological process GO terms and four KEGG pathways were associated with metabolic processes. One KEGG pathway was associated with olfactory transduction.
In the selective sweep analysis of CNVs, we identified 137 candidate selected CNVs associated with domestication (Supplementary Data 23). Annotation of the CNVs indicated the CNVs to be located in genes (Supplementary Table 8) or coincident with known QTLs 14 (Supplementary Data 24), which are functionally related to traits and biological processes such as follicular development and fertility (SLIT2) 15 , milk production (JAK2) 16 , wool production (KIF16B) 17 , adipogenesis (TCF7L1 and BCO2) 18,19 , and spleen size, oxygenated red blood cells and consequently high tolerance to hypoxia (PDE10A) 20 (Supplementary Notes). Also, we found divergent frequency distributions for seven deletions (overlapping with RFX3, AGMO, BCO2, LOC101112255, ADAMTSL3, and SGCZ) and three translocations (overlapping with GTF2I, CAMK4, and SGCZ) between the Asiatic mouflon and domestic sheep (Supplementary Table 9 and Supplementary Notes). Additionally, by comparing these 137 candidate domestication CNVs with the 144 domestication sweeps identified using both XP-CLR and π ratio analyses, we detected CNVs located within two selective sweeps and annotated three genes (i.e., BCO2, USP6NL, and LOC101112255; Supplementary Table 10).
Selective signatures during breeding and improvement. After domestication, selective signatures in sheep are expected to be engendered in different breeds through adaptation to a diverse range of environments and specialized production systems during breeding and improvement (Supplementary Fig. 8) 9,21 . In this context, we further compared the genomes of domestic breeds (i.e., the 36 landraces and six improved breeds; Supplementary Data 1) to detect signatures of positive selection during this process.
We calculated global F ST among the domestic breeds using a 50 kb sliding window and shift of 25 kb across genome, and identified 205 putatively selected genomic regions ( Fig. 3c Data 26). Annotation of these genes revealed functions associated with phenotypic and production traits including presence or absence of horns, pigmentation, reproduction, and body size (Supplementary Table 11 and Supplementary Data 27). We also observed genes functionally related to environmental adaptation, energy metabolism, and immune response, which may have been the targets of long-term natural selection 21,22 . Functional analysis of the 391 selected genes revealed significant enrichments for GO categories involved in four biological process categories including immune response, and immune system processes as well as 11 molecular function categories, such as cytokine activity and ATPase activity, coupled to transmembrane movement of ions, phosphorylative mechanism associated with energy metabolism, and immune function (Supplementary Data 28). The most significantly enriched pathway was cytokine-cytokine receptor interaction (Supplementary Data 28).
Notably, the selective region with the highest F ST value (F ST = 0.56) was located near the gene PAPPA2, which has been reported to be associated with fat deposition in humans 23 and has been identified as a candidate gene for milk, reproduction, and body size traits in cattle 24 . Comparison of allele frequencies at nonsynonymous SNPs in candidate selected genes revealed four (e.g., SPAG8, FAM184B, PDE6B, and PDGFD) with significantly differentiated allele frequencies among domestic sheep breeds (Supplementary Data 29).
Of these 391 selected genes, nine were also among the previously identified 59 most plausible domestication genes (Supplementary Table 6 and Supplementary Data 21) and they (RASGRF2, FBXO39, XAF1, GMDS, HOXA11, PDE6B, FLT1, EDN3, and RALY) (Supplementary Table 11) were linked to both domestication and breed-level genetic differentiation. In addition, 22 selected genes were confirmed to be under selection in our analyses for specific phenotypic traits such as reproduction, presence of horns, fat tail, wool fineness, nipple number, and ear size (see below; Supplementary Table 11). Moreover, 52 selected genes (e.g., SOCS2, EDA2R, PDE1B, PDGFD, HOXA10, etc.) have been shown to be under selection in sheep, humans, and other domestic animals in previous investigations (Supplementary Data 27). Additionally, a comparison of the selected genomic regions with previously reported QTLs 14 revealed 131 regions with high F ST values spanning the QTLs. These regions covered 19.97 Mb of the sheep genome and were found to be associated with morphological and production traits, such as reproductive seasonality, milk related traits, body weight, meat-related traits, teat number, tail fat deposition, and presence of horns (Supplementary Data 30).
Genetic mechanisms of the tail configuration trait. We implemented genome-wide selection tests between domestic breeds representing contrasting phenotypes for several traits that are relevant for sheep husbandry (Supplementary Table 12), such as morphological traits. Focusing on an iconic trait, tail configuration (Fig. 4a), we performed separate pairwise-population  1 1 21 3 1 4 1 5 1 6 1 7 1 8 1 9 2 0 2 1 2 2 2 3 2   We dissected the genomic architecture of the selected genes by calculating the frequency of the variant allele at non-synonymous SNPs. The frequencies of one variant allele each at genes PDGFD, XYLB, TSHR, and SGCZ as well as the genotype pattern located in the promoter region of PDGFD were different (Mann-Whitney, P < 0.001) between the fat-tailed (e.g., HDW, ALS, BSB, and TAN) and thin-tailed (e.g., SHE, Gotland sheep (GOT) and Finnsheep (FIN)) breeds (Fig. 4d, e, Supplementary Fig. 10 and Supplementary Data 39).
Remarkably, PDGFD was consistently selected by multiple comparisons (Supplementary Data 38). Transcriptome analysis among populations with different tail configurations (Supplementary Table 13) also identified PDGFD as significantly differentially expressed gene (log 2 (fold change) = 3.08; P adj = 0.045) between the fat-tailed/fat-rumped and the thin-tailed sheep (Supplementary Data 40). Furthermore, we detected four transcripts (i.e., transcripts I, II, III, and IV) of the PDGFD gene with the transcript I to be the most differentially expressed isoform between the thin-tailed and the fat-tailed/fat-rumped sheep (Fig. 4f), indicating its primary role in regulating fat deposition in tail. Furthermore, RT-PCR, qPCR, and western blot analyses demonstrated that gene expression level and protein level of PDGFD were consistently correlated negatively with fat deposition in sheep tail, with the highest level observed in the thin-tailed Merino sheep (MFW), followed sequentially by the small-tailed Han sheep (SXW), the large-tailed Han sheep (HDW), and fat-rumped Altay sheep (ALS) (Fig. 4g-j and Supplementary Fig. 11).
Selective and association signatures for other traits. Apart from tail configuration, we found many selected regions, novel functional genes, and non-synonymous SNPs related to the potentially selected genes, which may be responsible for traits such as reproduction, milk yield, wool fineness, meat production, and growth rate as well as for morphological traits including numbers of horns and nipples, pigmentation, and ear size ( To fine-map regions identified using selective sweep methodologies and search for direct evidence of genotype-phenotype associations, we performed GWAS for three quantitative traits (i.e., litter size and numbers of horns and nipples) with informative phenotypic records ( Supplementary Fig. 21). Using a panel of samples from multiple breeds and high-quality SNPs as well as the mixed linear model (MLM), we identified 600, 989, and 1969 significant GWAS signals for litter size (109 samples from 11 breeds; 14,574,050 SNPs), number of horns (146 samples from 15 breeds; 14,556,831 SNPs), and number of nipples (123 samples from 13 breeds; 14,415,949 SNPs) with the thresholds of −log 10 (P value) = 6, 6, 4, respectively (Fig. 5c,  Supplementary Fig. 22 and Supplementary Data 46). Furthermore, we detected 20, 56, and one of these respective GWAS signals to be overlapped with selective sweeps detected for the three traits (Supplementary Data 47), suggesting the significance of these genomic regions in shaping the traits.
Except for previously reported major candidate genes (e.g., BMPR1B, INHBB, and ESR1) 26 , annotation of the significant GWAS signals revealed that those for litter size were mapped to a number of novel genes, such as NOX4, IRF2, PDE11A, ZFAT, ZFP91, TENM1, BICC1, LRRTM3, CTNN3, SMYD3, KCNN3, and d Selective regions associated with tail configuration by XP-CLR using the SNP data with the threshold XP-CLR ≥ 8.26. Candidate genes overlapping with the regions, which are significantly selected by XP-CLR & ln(π ratio)/ln(2), XP-CLR & ln(π ratio)/ln(2) & iHS, and XP-CLR & ln(π ratio)/ln(2) & HKA are marked by gray, orange, and blue colors, respectively. Below this plot, genes near the peaks are indicated by green boxes. The pie charts represent the spectrum of allele frequencies at the non-synonymous loci of PDGFD in populations of different tail configurations. The type of variant allele is indicated in blue, while the reference allele in pink. e Genotype patterns for the promoter region of PDGFD among 11 fat-tailed/rumped, 11 thin-tailed sheep, and Asiatic mouflon. f Structures and expression levels of four isoforms of PDGFD. Expression levels are shown in varying shades of yellow color. g, i Expression pattern of control gene β-actin and target gene PDGFD in tail fat examined by RT-PCR g and western blot analysis i. h, j The relative expressions of PDGFD in tail fat by real-time PCR (qPCR) h and western blot analysis j. k Adipogenesis signaling pathway 46 and the inhibitory function of PDGFD in differentiation of white adipocytes 45 by activating PDGFRβ signaling 44 . All experiments were repeated three times with similar results. Samples derived from the same experiment and the blots were processed in parallel. g-j Experiments were performed with the control sample (the thin-tailed sheep; MFW) and target samples (long fat-tailed sheep (HDW), fat-rumped sheep (ALS) and short fat-tailed sheep (SXW)). The data in h and j are presented as the mean ± SD, n = 3 biologically independent samples; groups with significant differences (*P < 0.05; **P < 0.01) were performed by two-tailed unpaired t-test. Source data are provided as a Source Data file.  Supplementary Fig. 22a), which play roles in various reproductive functions, including embryogenesis, uterine remodeling, follicular development and ovulation 27,28 . Focusing on 600 significant SNPs, we found 323 SNPs in intergenic regions, one in the downstream and 276 SNPs in the coding regions (only one SNP in the exon but 275 in the introns). For the number of horns, several GWAS peaks were situated in HOXD1, HOXD3, and RXFP2 (Fig. 5c) as well-characterized functional genes for the polycerate and polled phenotypes in sheep 29 . For the number of nipples, most of the significant GWAS signals were located in genes associated with breast cancer, including five genes (LRP1B, GRM3, MACROD2, SETBP1, and GPC3) reported previously 30 .
Finally, we examined a total of 3,558 significant association signals with previously reported QTLs for reproduction, numbers of horns and nipples, and mapped 121 signals in five QTLs responsible for reproductive seasonality and total lambs born (permutation test, P < 0.001), 87 signals in one QTL associated with horns (permutation test, P < 0.001), and 31 signals in four QTLs related to teat placement, udder depth, udder shape and udder attachment (Supplementary Table 5 and Supplementary Data 49). We also implemented GWAS analyses using the CNV data for litter size, numbers of horns and nipples (Fig. 5b,  Supplementary Fig. 23 and Supplementary Tables 14 and 15). We detected a total of 11 significant association signals and associated functional genes for the three traits (litter size: SMARCA1 and APP; number of horns: RXFP2; and number of nipples: GPC5) (Supplementary Notes). Several significant CNVs found by GWAS were in common with CNVs identified by selective sweep analysis. Two of these common CNVs were associated with reproductive traits while three CNVs with horn related traits (Supplementary Table 16).

Discussion
In this study, we re-sequenced the whole genomes of 248 wild, landrace, and improved sheep with an emphasis on local breeds with distinct phenotypes that have not been studied previously at the genomic level. Our exploration of the sheep variome and selective sweeps focused specifically on domestication and selective breed formation.
Deciphering the genetic basis of animal domestication is an active research area. The availability of whole-genome sequences provides an opportunity to study this at the gene mutation level. Such genomic studies have been implemented in other domestic animals 33,34 and recently also in sheep 1,2 . These studies were based on genome sequences of low-to-medium coverages (8.4-17.2× in ref. 2 ; 12-14× in ref. 1 ). We have complemented this work using a hierarchically structured breed panel and highdepth whole-genome sequencing (mean depth of 25.7×).
We observed a lower level of genomic diversity in domestic breeds than their wild ancestors. This suggested that a substantial proportion of the genomic variation has been lost during and after domestication, whereas the genomic diversity in landraces has been largely retained in improved breeds. It should be noted that the nearly identical nucleotide diversity between landraces and improved breeds was also consistent with the observation that very strong positive selective pressure on modern breeds has only been in practice over the last~200 generations 9 . Similar patterns of changes in genomic diversity through domestication have been observed in yak 35   domestication. Results of functional enrichment analysis of the 261 candidate genes were in agreement with previous evidence that metabolic process and olfactory transduction were involved as primary functional categories in sheep domestication 2 . Furthermore, functions of the domestication-related genes and significant differences in non-synonymous SNP allele frequencies between Asiatic mouflon and domestic sheep provided additional evidence for signatures of selection on these genes and associated traits (e.g., reproduction, immunity, fat, photoreceptor, and olfaction) brought about by domestication. Nevertheless, by annotation of the 260 fixed SNPs (derived allele frequency ≥ 0.95) in the 209 putative selective sweeps, we detected 206 SNPs located in the non-coding regions but only two SNPs in the exon regions including a single non-synonymous mutation. This finding suggested domestication as a quantitative trait (e.g., litter size) to be affected mostly by mutations in non-coding regions 36,37 , whereas very few mutations were non-synonymous. In addition to the functional non-synonymous SNPs, we observed differentiated frequencies in SVs between Asiatic mouflon and domestic sheep (Supplementary Tables 9). Thus, SVs in functional genes could also account for the changes in phenotypic traits during domestication. In particular, our results indicated that not only SNPs but also CNVs associated with adipogenesis (BCO2) 19 and proteostasis (LOC101112255) 38 have been under selection during domestication. Interestingly, as SVs in GTF2I have been found to be linked to Williams-Beuren syndrome in humans 39 and hypersocial behavior (i.e., a feature of the domestication syndrome 40 ) in domestic dogs 41 , the specific SV present in GTF2I may account for some behavioral differences between Asiatic mouflon and domestic sheep gained or lost during domestication. This work may inform ongoing and future analysis of ancient DNA in order to pinpoint more accurately the origin and time of sheep domestication and associated impact at the genomic level. We noted that Asiatic mouflon have several subspecies (e.g., Ovis orientalis gmelini, O. o. ophion, and O. o. laristanica) 42 and are distributed in various geographic areas such as Iran, Turkey, Azerbaijan, and Cyprus, therefore a comprehensive sampling of all of them would be necessary in future investigations. It would also be interesting to compare genetic mechanisms responsible for specific domestication traits, and identify general patterns across different species of domestic animals 1,6 .

K C N N 3 C D 9 6 Z F A T P D E 1 1 A S M Y D 3 Z F P 9 1 N O X 4 B IC C 1 C T N N A 3 L R R T M 3 IR F 2 T E N M 1 R X F P 2 P H G D H H M G C S 2 R O B O 2 G L IS 3 L R P 1 B K D M 3 A F S H R G R M 3 C S N 1 S 1 M A C R O D 2 S E T B P 1 G P C 3
In addition, we detected a few selective signatures associated with phenotypic and production traits during breeding and improvement. Functions of the candidate genes implied potential roles of human-induced changes in growth rate (SPAG8), reproduction (FAM184B), photoreceptor development (PDE6B), and tail configuration (PDGFD) during the development of specific breeds (Supplementary Notes). In particular, a strong selective signature located near the PAPPA2 gene implied intense artificial selections for body fat or types of tail (e.g., fat and thin tails) and production traits in sheep towards unique breeds. It is worth noting that the most significantly enriched pathway for the 391 selected genes identified from F ST analysis (Supplementary Data 26) was cytokine-cytokine receptor interaction (ACVR2B, TNFSF4, CCL25, etc.) (Supplementary Data 28). Previous studies have revealed that, instead of artificial selection, this important pathway for immune function can also be impacted by long-term natural selection 21,22 .
We investigated various phenotypic traits using an integrated analysis of whole-genome selection scans and GWAS. Although some of the candidate genes have been functionally associated with various traits in previous investigations 3 , we used highcoverage sequencing data to identify a number of novel candidate genes. In particular, we mapped novel candidate genes associated with popular traits that were studied previously (e.g., litter size, ear size, and coat color) or unique traits that were rarely studied in sheep and other livestock (tail configuration and numbers of horns and nipples). Moreover, our findings will help to narrow down the functional sub-regions within the QTLs and pinpoint the causal genes associated with these breeding-related traits. In particular, we identified several previously reported and a few novel genes to be associated with the tail configuration in sheep (Supplementary Data 38). Previous studies revealed that aminoacid changes in PDGFD (cysteine/tyrosine) have an important role in fat metabolism and adipogenesis in humans 43 . We envisage that non-synonymous variants with deviating allele frequencies between wild and domestic sheep and large effects on specific phenotypes in domestic sheep might be useful, for instance, as targets in CRIPSR/Cas experiments. Transcriptome analysis showed significant differential expressions of PDGFD transcripts among populations with different tail configurations (Supplementary Data 40). This may be ascribed to the distinct genotype pattern in the promoter region of PDGFD (Fig. 4e). PDGFD functions by causing dimerization and further activating PDGF receptor PDGFRβ 44 . Early studies showed that PDGFRβ signaling has an essential role in inhibiting differentiation of white adipocytes by regulating the expression of PPARγ2 and C/ EBPα 45 , which were identified as the key transcriptional regulators of adipogenesis 46 (Fig. 4k). Therefore, our results provide in-depth insights into the genomic architecture and molecular mechanism for tail configuration in sheep at the genotype, variant allele, transcript, and protein levels.
All our efforts have resulted in a unique data resource in terms of the sheep variome and selective sweeps with different categories of genetic markers, allele distributions in different breeds, and associations with phenotypes with different degrees of experimental validation. This will underpin more-accurate identification of causative gene variants in the near future and facilitate novel breeding strategies, like marker-assisted or genomic selection and genome editing targeting favorable traits towards a cost-effective and environmentally friendly sheep industry.

Methods
Sample collection, DNA extraction, and sequencing. Blood samples were collected from a total of 248 individuals comprising 232 domestic sheep (O. aries) and 16 wild sheep (Asiatic mouflon O. orientalis). The domestic sheep samples represent 36 landraces (172 individuals) and six improved breeds (60 individuals) with different geographic origins from Asia, Europe, Africa, and the Middle East (Fig. 1a). More specifically, the domestic samples represent various geographic origins, morphological characteristics, and production traits (Supplementary Data 1). Breed origins of the domestic sheep samples included populations from geographic areas underrepresented in earlier work (China, Afghanistan, Iran, Iraq, Azerbaijan, South Africa, Ethiopia, Burkina Faso, Niger, Nigeria, Chad, and Cameroon) as well as Germany, Spain, England, Finland, France, Scotland, Sweden, and the Netherlands (Fig. 1a and Supplementary Data 1). All the domestic sheep were typical of the breeds and unrelated according to pedigree records or herdsman's information (Fig. 1a and Supplementary Data 1). The Asiatic mouflon were collected from captivity in Iran, which is within the putative geographic center of sheep domestication. To minimize potential bias as a result of overrepresentation of local effects (e.g., inbreeding), individuals from different locations were sampled. A full description of the samples is detailed in Supplementary Data 1. Genomic DNA was extracted following the standard phenol-chloroform extraction procedure. For genome sequencing, at least 0.5 μg of genomic DNA from each sample was used to construct a library with an insert size of~350 bp. Paired-end sequencing libraries were constructed according to the manufacturer's instructions (Illumina Inc., San Diego, CA, USA) and sequenced on the Illumina HiSeq X Ten Sequencer (Illumina Inc.).
Sequence read mapping. We obtained~82.86 Gb of raw sequences for each sample, giving an average depth of 25.7× coverage for clean reads (17.2-37.0×) (Supplementary Data 2). The 150-bp paired-end reads were mapped onto the sheep reference genome Oar v.4.0 (https://www.ncbi.nlm.nih.gov/assembly/ GCF_000298735.2) with the Burrows-Wheeler Aligner v.0.7.8 (ref. 47 ) using the default parameters. Mapping results were then converted into the BAM format and sorted with SAMtools v.1.3.1 (ref. 48 ). Duplicate reads were removed using SAMtools. If multiple read pairs had identical external coordinates, only the pair with the highest mapping quality was retained.
SNP calling, validation, and annotation. After mapping, we performed SNP calling separately for the two sets of samples (see below) using the Bayesian approach implemented in SAMtools and Genome Analysis Toolkit (GATK) v.3.7 (ref. 49 ), with all individuals in each set simultaneously. One set included the 228 wild/domestic sheep with at least five samples per breed/population, which was used in all analyses, whereas the other set consisted of the 20 domestic sheep from the Middle East and Africa with one individual per breed, which was used to explore population structure and demographic history of domestic sheep in the Old World and to identify selective signatures associated with domestication and improvement (e.g., global F ST analysis). Only SNPs detected by both methods were kept for further analyses. The detailed processes were as follows: (i) For the GATK, the UnifiedGenotyper parameters -stand_emit_conf and -stand_call_conf were both set as 30. The same aligned BAM files were used in SNP calling through the SAMtools mpileup package; and (ii) For filtering using the command parameters -mis 0.1-maf 0.05 -qd 2 -fs 60 -mq 40 -dp_min 6 -dp_max 120 -DP_min 100 -DP_max 30000 -gq 20 -MQRankSum -12.5 -ReadPosRankSum -8.0, the common sites were first identified by the GATK and SAMtools using the SelectVariants package, and then SNPs with missing rates ≥0.1 and minor allele frequencies (MAF) <0.05 from the three groups (Asiatic mouflon, landraces, and improved breeds) were filtered out from further analysis. For Asiatic mouflon and each breed of domestic sheep, we estimated the site frequency spectrum (Supplementary Fig. 24) based on individual genotype likelihoods assuming Hardy-Weinberg equilibrium using the ANGSD v.0.915 (ref. 50 ) with the parameters -dosaf 1 -fold 1 -maxIter 100.
To validate the SNPs detected, we first compared the identified set with the O. aries dbSNP v.151 at the National Center for Biotechnology Information (NCBI; http://www.ncbi.nlm.nih.gov/SNP). Next, we compared the genotypes of the called SNPs with those on the Ovine Infinium HD BeadChip array (~600 K SNPs) (Illumina, San Diego, CA) for 223 samples with available chip data (Supplementary Methods). In addition, to assess the performances of the GATK and SAMtools variant calling methods, we employed a false-positive measure by determining the rate of the monomorphic reference loci on the Ovine Infinium HD BeadChip array that were erroneously called as variant loci by the variant calling methods. We calculated the false-positive rates as the number of false heterozygous SNPs divided by the total number of homozygous reference loci 51 .
SNPs were annotated using the ANNOVAR v.2013-06-21 (ref. 52 ) based on the sheep reference genome Oar v.4.0 and then categorized as variations in exonic regions, splicing sites, intronic regions, upstream and downstream regions, and intergenic regions. Those in exons were further classified into synonymous or nonsynonymous SNPs.
Identification of Indels, CNVs, and SVs. Similar to SNP calling, the calling of INDELs was conducted using SAMtools with minimum depth ≥4 and GQ >20, and only INDELs <100 bp were retained. CNVs were detected using both CNVnator v.0.3.2 (ref. 53 ) and DELLY v.0.7.9 (ref. 54 ). For CNVnator, the analyses were performed on the BAM files with a bin size of 100 bp and with the length >200 bp (Supplementary Methods). For DELLY, the analyses were conducted with default parameters, and deletions and duplications were considered to be CNVs (Supplementary Methods). Only the CNV calls with >50% of their lengths being overlapped between the two approaches were retained in the final set of CNVs. The CNVs that overlapped with gaps or genomic repeats were removed, and the remaining CNVs were segregated into short CNV bins (≥100 bp) across the genomes among the 248 individuals for subsequent analyses 25 .
Population genetics analysis. After filtering, we generated a set of SNPs for the following analyses. First, an individual-based neighbor-joining (NJ) tree was constructed for all the samples based on the nucleotide p-distance matrix using the TreeBeST v.1.9.2 (ref. 57 58 ). Furthermore, population structure was assessed using the default setting in the ADMIXTURE v.1.23 (ref. 59 ). The number of assumed genetic clusters K ranged from 2 to 7. To construct a NJ tree for the four subgroups of domestic sheep (i.e., African, East Asian, Middle Eastern, and European groups; see "Results"), 1-2 individuals from each landrace and two individuals from each sampling site of Asiatic mouflon were selected, totaling 50 individuals (i.e., 10 for each group; Supplementary Table 17). SNPs for the 50 individuals were extracted from the dataset of landraces and Asiatic mouflon. To mitigate the possible effect of LD, we implemented LD pruning using the parameter-indep-pairwise (50 5 0.4) in PLINK v.1.07 (ref. 60 ). To eliminate the potential influence of selective SNPs, we only retained the SNPs located 150 kb away from genes and without missing genotypes. Eventually, a final set of 59,943 SNPs for the 50 individuals were kept for the construction of NJ tree. Reynolds genetic distances (ref. 61 ) among the five groups were calculated using Arlequin v3.5.2.2 (ref. 62 Table 18). A NJ tree was then constructed based on the Reynolds genetic distances with 1000 bootstraps using PHYLIP v.3.695 (ref. 63 ) and visualized using FigTree v.1.4.3. The parameter r 2 (ref. 64 ) for LD was calculated for pairwise SNPs within each chromosome using PLINK v.1.07 60 with the parameters (-ld-window-r 2 0 -ld-window 99999 -ld-window-kb 500). The average r 2 values were calculated for each length of distance and the whole-genome LD was averaged across all chromosomes. The LD decay plot was depicted against the length of distance using the R script (http:// www.r-project.org). Nucleotide diversity (π) and global F ST were calculated using the vcftools v.0.1.14 (ref. 65 ). The F ST values between populations were estimated using the ARLSUMSTAT implemented in the Arlequin v3.5.2.2 (ref. 62 ), with a sliding window of 50 kb. The genomic SNP data of variant call format (VCF) were converted into the Arlequin format (arp) using the VCF2Arlequin python script 62 . Statistical significance (P values) of the F ST values were tested through 100,000 Markov chains following 10,000 burn-in steps. The average F ST and associated P values over all sliding windows were regarded as the values at the wholegenome level.

) (Supplementary
Estimates of effective population size. We used the PSMC 66 method to estimate changes in effective population size (N e ) over the last one million years. The PSMC analysis was implemented in each of the 248 samples. The parameters were set as follows: -N30 -t15 -r5 -p '4 + 25*2 + 4 + 6', with the filtering criteria of read depth for each SNP as six at the individual level. An average mutation rate (μ) of 2.5 × 10 −8 per base per generation and a generation time (g) of 3 years 67 were used for the analysis. We also inferred recent Ne using the SNeP v.1.0 (ref. 68 ) with default settings. SNPs with missing data and a MAF smaller than 5% were excluded from the analysis. The different SNP marker distance bins for r 2 analysis were used to obtain different estimates of N e at t = 1/2c generations ago.
Detection of selective sweeps. For SNPs, we performed tests for selective sweeps during domestication and breeding using two approaches based on the SNPs with less than 10% missing data: the XP-CLR approach implemented in the XP-CLR v.1.0 (ref. 69 ), and by the comparison of π ratios calculated using the vcftools v.0.1.14 (ref. 65 ). To detect genomic regions under selection during domestication, we calculated the ln(π-O. orientalis /π-Landrace )/ln (2). Also, we estimated the ln (π-Control /π-Target )/ln(2) between populations of domestic sheep with contrasting phenotypes for a specific target trait. The specific populations involved in comparisons between wild and domestic sheep and pairwise comparisons between domestic populations for detecting the signals associated with particular traits are shown in Supplementary Table 12. Values of π were calculated with a 50 kb sliding window and a 25 kb sliding step. For the XP-CLR approach, a 0.5 cM sliding window with a spacing of 2 kb across the whole genomes were used for scanning, and 200 SNPs were assayed in each window with the parameters -w1 0.005 200 2,000 chrN -p0 0.95. To assess the statistical significance of the XP-CLR value for each window, we first estimated the proportion of SNPs with extreme XP-CLR values (i.e., top 1%) in the sliding windows, and then calculated the P values from the empirical distribution of the proportion scores obtained with these windows. In each comparison, the genomic regions in the top 1% XP-CLR values and ln(π ratio)/ln(2) values across the whole-genome were considered to be the selective sweeps.
Moreover, we estimated the iHS across the genomes of Asiatic mouflon and different groups of domestic sheep populations using the Selscan v.1.2.0 (ref. 70 ) after filtering all missing data, with 50 kb sliding windows and 25 kb stepwise, a recombination rate of 1 cM Mb −1 (ref. 9 ) and default parameters -max-extend 1,000,000 -max-gap 200,000 -cutoff 0.05 (Supplementary Methods). We computed the proportions of SNPs with normalized |iHS | >2 in non-overlapping windows, and identified those windows within the top 5% empirical cutoff (i.e., above the 95th percentile of genome-wide distribution) 71 in the tested group as the signals of positive selection. We also employed the HKA test 72 to identify the selective signals associated with domestication and specific traits using Asiatic mouflon as an outgroup after filtering all missing data (Supplementary Methods). We calculated the χ 2 statistic in 50 kb sliding windows and shift of 25 kb across the genome to find potential selective signals deviating from genome-wide neutral expectations. Two loci were analyzed each time, one was the 50 kb window taken from the tested genome and the other was the virtual neutral 50 kb window in terms of the average value of nucleotide statistics in the whole genome. After application of the HKA test for each sliding window, the χ 2 statistic used to measure the goodness-of-fit was obtained and subsequently used to identify the selective signals.
For CNVs, we calculated a statistic V ST , an analog to F ST . V ST estimates population differentiation based on the quantitative intensity data and varies from 0 to 1 (ref. 25 ). The statistic V ST of each CNV region was calculated to detect the selective signals between different comparisons 25  GWAS. Association analyses of litter size, numbers of horns and nipples were performed using the MLM in the GEMMA v.0.96 (ref. 73 ) based on a panel of 109, 146, and 123 samples collected from 11, 15, and 13 breeds, respectively (Supplementary Fig. 21). The effect of population stratification was corrected by adjusting the first three principal components (PCs) as derived from the whole-genome SNPs, and the proportion of variance explained by the markers was calculated using TASSEL v.5.0 (ref. 74 ). To avoid potential false positives in multiple comparisons, the whole-genome significance threshold was adjusted via the Bonferroni test 75 . For SNPs, we set the thresholds as −log 10 (P value) = 6, 6, and 4 for litter size, numbers of horns and nipples, respectively. For CNVs, we set the thresholds as −log 10 (0.05/total CNVs) = 5.28, 5.29 for litter size and number of horns, respectively, but −log 10 (P value) = 4 for number of nipples using GEMMA v.0.96 based on the genotypes of CNVs selected by the DELLY and CNVnator. In addition, the quantile-quantile (Q-Q) plots of the MLM for individual traits were implemented in R Bioconductor.
Permutation test for QTL overlaps. We performed permutation test to check if the overlaps between selective sweeps/GWAS hits and QTL regions were significantly different from those expected at random. To this end, we used BEDTools v2.26.0 shuffle to generate simulated data sets by randomly selecting genomic regions of equal number and size to the observed selective sweeps/GWAS hits in the sheep genome, and we replicated this process 10,000 times 76 . We compared the number of overlaps between the observed selective sweeps/GWAS hits and the QTL regions with the distribution of overlap statistics between the simulated selective sweeps/GWAS hits data sets and the QTL regions, and calculated the statistical significance of P values (i.e., the probability that a higher number of overlaps would be observed by chance).
RNA-Seq analysis. We collected tail adipose tissues from thin-tailed, short fattailed, long fat-tailed, and fat-rumped sheep for RNA-Seq analysis (Supplementary Table 13). Each tail type included three independent samples from different individuals as biological replicates. We used the Trizol RNA Reagent (Takara, Dalian, China) to extract total RNA from the tissues and measured the concentration and integrity of the RNA with the Agilent 2100 RNA 6000 Nano Kit (Agilent Technologies, Waldbronn, Germany). Subsequently, the libraries of mRNAs were constructed using the NEBNext® Ultra TM RNA Library Prep Kit (New England Biolabs, Ipswich, MA, USA) and sequenced using the Illumina HiSeq X Ten System to generate 150-bp paired-end reads. We used the HISAT v2.1.0, StringTie v2.0, and Ballgown package in R version 3.5.3 (ref. 77 ) to map the paired-end reads to the sheep reference genome, assemble the reads, and estimate the gene expression levels, respectively. The number of reads matched to an expressed gene was standardized as fragments per kilobase of exon per million mapped fragments (FPKM) values. We employed the stattest function in the Ballgown package 77 to search for transcripts that were differentially expressed between the thin-tailed breed (Chinese Merino sheep) and fat-tailed/fat-rumped breeds (Small-tailed Han sheep, Large-tailed Han sheep, and Altay sheep), following correction for any differences in expression owing to population variables. This allowed us to get the confounder-adjusted fold changes between the two tested groups. The genes that exhibited |log 2 (fold change)| ≥ 2 and adjusted P ≤ 0.05 in the comparisons between fat-tailed/fat-rumped and thin-tailed individuals were considered as differentially expressed genes.
Gene expression and western blot analyses of PDGFD gene. We examined the gene expression level of PDGFD in the adipose tissues from the thin-tailed, short fat-tailed, long fat-tailed, and fat-rumped sheep through reverse transcription PCR (RT-PCR) and qPCR. The three biological replicates of adipose samples from different individuals for each tail type were used. The total RNA was extracted using the Trizol RNA Reagent (Takara, Dalian, China) and were treated with RNase-free DNase I to remove DNA using the RapidOut DNA Removal Kit (Thermo Fisher Scientific, Waltham, MA, USA). The first-strand cDNA was synthesized using the RevertAid First Strand cDNA Synthesis Kit (Thermo Fisher Scientific, Waltham, MA, USA). Following the manufacturer's instruction, 500 ng of RNA was reverse transcribed as the template for RT-PCR in 40 μl volume (including 20 μl RNA, 2 μl Random Hexamer Primer (100 μM), 8 μl 5× Reaction Buffer (including 250 mM Tris-HCl (pH 8.3), 250 mM KCl, 20 mM MgCl 2 , 50 mM DTT), 2 μl RiboLock RNase Inhibitor (20 U μl −1 ), 4 μl 10 mM dNTP Mix, 2 μl RevertAid RT (200 U μl −1 ) and 2 μl nuclease-free water) and a thermocycling condition at 25°C for 5 min, 42°C for 60 min, and 70°C for 5 min. Subsequently, the qPCR with SYBR Green (Promega, Madison, WI, USA) was performed on the QuantStudio TM 6 Flex Real-Time PCR System (Life Technologies, Carlsbad, CA, USA) using the first-strand cDNA and the primers designed based on the 5′-and -3′ end sequences of PDGFD gene (PDGFD F: 5′-GCGGATGCTCTGGACAAA and PDGFD R: AAGGAGGCAGCGTGGAAA-3′). The qPCR reactions and conditions were set as those described above. Each qPCR was run three times for one sample as technical replicates. Based on the qPCR results, the expression level was calculated using the 2 −ΔΔCt method 78 and normalized according to the internal control β-actin gene. The primers for β-actin were β-actin F (5′-CCAACCGT-GAGAAGATGACC) and β-actin R (CCCGAGGCGTACAGGGACAG-3′).
Protein was extracted from the tail adipose tissue using the Total Protein Extraction Kit (Huaxingbio, Beijing, China). The protein extract was mixed with an equal amount of sample buffer and then separated on 10% sodium dodecyl sulfate-polyacrylamide gel electrophoresis (SDS-PAGE) gels (60 μg per lane). The SDS-PAGE-separated proteins were electrophoretically transferred to a polyvinylidene fluoride (PDVF) membrane and then incubated for 3 hours at room temperature in blocking buffer (5% BSA in PBS- Tween 20). Immunodetection was carried out with the Rabbit Anti-beta Actin antibody (ab8227, Abcam, dilution 1:1,000), Anti-SCDGFB/PDGFD antibody (ab181845, Abcam, dilution 1:1,000) and Goat Anti-Rabbit IgG H&L (ab205718, Abcam, dilution 1:10,000). The blot signals were imaged using Tanon 6100 Chemiluminescent Imaging System and quantified using ImageJ (NIH) software. Photoshop CS6 were used to crop images from unprocessed images.
Phenotyping. Individual phenotype for traits, such as coat color, classes of fiber fineness, ear size, numbers of nipples and horns, and tail configurations were recorded for all the breeds whenever possible during sampling. Five different coat colors or color patterns, including white, white body with black head, black, brown, and gray were recorded for the animals sampled. The number of nipples ranged from 2 (normal) to 5 (selected). The horn phenotypes varied from polled to horned animals with 2-5 horns. The wool was graded into three classes (coarse, fine, and super fine) according to the British Wool Grading System (http://www.eytest.com/ ey31f1.html). The tails were categorized into short fat-tailed, long fat-tailed, thintailed, and fat-rumped types according to the shape of tails of the animals as well as the recorded information for the breeds. Reproductive traits included number of litter per year, litter size in each birth, and seasonal or non-seasonal estrus extracted from the breeding records.
Ethics statement. All animal work was conducted according to a permit (no. IOZ13015) approved by the Committee for Animal Experiments of the Institute of Zoology, Chinese Academy of Sciences (CAS), China. For domestic sheep, animal sampling was also approved by local authorities where the samples were taken. For Asiatic mouflon, we collected peripheral blood samples from 14 captive Asiatic mouflon after receiving authorization for research from the Department of Environmental Protection in Iran (no. 93/34089). For other two Asiatic mouflon samples from Shahr-e Kord, Iran, sampling procedure was also approved by the governorate of Chaharmahal and Bakhtiari of Iran (no. 97.32.43.33165).