Whole-Genome Resequencing of Worldwide Wild and Domestic Sheep Elucidates Genetic Diversity, Introgression, and Agronomically Important Loci Feng-Hua Lv,†,1 Yin-Hong Cao,†,2,3 Guang-Jian Liu,†,4 Ling-Yun Luo,1 Ran Lu,1 Ming-Jun Liu,5 Wen-Rong Li,5 Ping Zhou,6,7 Xin-Hua Wang,6,7 Min Shen,6,7 Lei Gao,6,7 Jing-Quan Yang,6,7 Hua Yang,6,7 Yong-Lin Yang,6,7 Chang-Bin Liu,6,7 Peng-Cheng Wan,6,7 Yun-Sheng Zhang,6,7 Wen-Hui Pi,6,7 Yan-Ling Ren,8 Zhi-Qiang Shen,8 Feng Wang,9 Yu-Tao Wang,10 Jin-Quan Li,11 Hosein Salehian-Dehkordi,2,3 Eer Hehua,12 Yong-Gang Liu,13 Jian-Fei Chen,1 Jian-Kui Wang,1 Xue-Mei Deng,1 Ali Esmailizadeh ,14 Mostafa Dehghani-Qanatqestani,14 Hadi Charati,14 Maryam Nosrati,15 Ondrej Stepanek,16 Hossam E. Rushdi ,17 Ingrid Olsaker,18 Ino Curik,19 Neena A. Gorkhali,20 Samuel R. Paiva,21 Alexandre R. Caetano,21 Elena Ciani,22 Marcel Amills,23,24 Christina Weimann,25 Georg Erhardt,25 Agraw Amane,26,27 Joram M. Mwacharo,28,29 Jian-Lin Han,30,31 Olivier Hanotte,27,32 Kathiravan Periasamy,33 Anna M. Johansson,34 J on H. Hallsson,35 Juha Kantanen,36 David W. Coltman,37 38,39 Michael W. Bruford , Johannes A. Lenstra,40 and Meng-Hua Li *,1 1 College of Animal Science and Technology, China Agricultural University, Beijing, China 2 CAS Key Laboratory of Animal Ecology and Conservation Biology, Institute of Zoology, Chinese Academy of Sciences (CAS), Beijing, China 3 College of Life Sciences, University of Chinese Academy of Sciences (UCAS), Beijing, China 4 Novogene Bioinformatics Institute, Beijing, China 5 Animal Biotechnological Research Center, Xinjiang Academy of Animal Science, Urumqi, China 6 Institute of Animal Husbandry and Veterinary Medicine, Xinjiang Academy of Agricultural and Reclamation Sciences, Shihezi, China 7 State Key Laboratory of Sheep Genetic Improvement and Healthy Production, Xinjiang Academy of Agricultural and Reclamation Sciences, Shihezi, China 8 Shandong Binzhou Academy of Animal Science and Veterinary Medicine, Binzhou, China 9 Institute of Sheep and Goat Science, Nanjing Agricultural University, Nanjing, China 10 College of Life and Geographic Sciences, Kashi University, Kashi, China 11 College of Animal Science, Inner Mongolia Agricultural University, Hohhot, China 12 Grass-Feeding Livestock Engineering Technology Research Center, Ningxia Academy of Agriculture and Forestry Sciences, Yinchuan, Article China 13 College of Animal Science and Technology, Yunnan Agricultural University, Kunming, China 14 Department of Animal Science, Faculty of Agriculture, Shahid Bahonar University of Kerman, Kerman, Iran 15 Department of Agriculture, Payame Noor University, Tehran, Iran 16 Department of Virology, State Veterinary Institute Jihlava, Jihlava, Czech Republic 17 Department of Animal Production, Faculty of Agriculture, Cairo University, 12613 Giza, Egypt 18 Department of Preclinical Sciences and Pathology, Faculty of Veterinary Medicine, Norwegian University of Life Sciences, Ås, Norway 19 Department of Animal Science, Faculty of Agriculture, University of Zagreb, Zagreb, Croatia 20 Animal Breeding Division, National Animal Science Institute, Nepal Agriculture Research Council (NARC), Kathmandu, Nepal 21 Embrapa Recursos Geneticos e Biotecnologia, Parque Estaç~ao Biologica, PqEB, Brasılia, DF, Brazil 22 Dipartimento di Bioscienze, Biotecnologie e Biofarmaceutica, Universita degli Studi di Bari Aldo 24 Moro, Bari, Italy 23 Department of Animal Genetics, Center for Research in Agricultural Genomics (CRAG), CSIC-IRTA-UAB-UB, Campus de la Universitat Aut onoma de Barcelona, Bellaterra, Spain 24 Department of Animal Sciences, Universitat Autonoma de Barcelona, Bellaterra, Spain 25 Department of Animal Breeding and Genetics, Justus-Liebig-University Giessen, Giessen, Germany 26 Department of Microbial, Cellular and Molecular Biology, Addis Ababa University, Addis Ababa, Ethiopia 27 LiveGene Program, International Livestock Research Institute, Addis Ababa, Ethiopia 28 Small Ruminant Genomics, International Centre for Agricultural Research in the Dry Areas (ICARDA), Addis Ababa, Ethiopia ß The Author(s) 2021. Published by Oxford University Press on behalf of the Society for Molecular Biology and Evolution. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons. org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited. Open Access Mol. Biol. Evol. 39(2):msab353 doi:10.1093/molbev/msab353 Advance Access publication December 10, 2021 1 Lv et al. . doi:10.1093/molbev/msab353 MBE 29 CTLGH and SRUC, The Roslin Institute Building, Easter Bush Campus, Edinburgh, Scotland 30 CAAS-ILRI Joint Laboratory on Livestock and Forage Genetic Resources, Institute of Animal Science, Chinese Academy of Agricultural Sciences (CAAS), Beijing, China 31 Livestock Genetics Program, International Livestock Research Institute (ILRI), Nairobi, Kenya 32 School of Life Sciences, University of Nottingham, University Park, Nottingham, United Kingdom 33 Animal Production and Health Laboratory, Joint FAO/IAEA Division, International Atomic Energy Agency (IAEA), Vienna, Austria 34 Department of Animal Breeding and Genetics, Faculty of Veterinary Medicine and Animal Science, Swedish University of Agricultural Sciences, Uppsala, Sweden 35 Faculty of Natural Resources and Environmental Sciences, Agricultural University of Iceland, Borgarnes, Iceland 36 Production Systems, Natural Resources Institute Finland (Luke), Jokioinen, Finland 37 Department of Biological Sciences, University of Alberta, Edmonton, Alberta, Canada 38 School of Biosciences, Cardiff University, Cathays Park, Cardiff, Wales, United Kingdom 39 Sustainable Places Research Institute, Cardiff University, Wales, United Kingdom 40 Faculty of Veterinary Medicine, Utrecht University, Utrecht, the Netherlands † These authors contributed equally to this work. *Corresponding author: E-mail: menghua.li@cau.edu.cn. Associate editor: Katja Nowick Abstract Domestic sheep and their wild relatives harbor substantial genetic variants that can form the backbone of molecular breeding, but their genome landscapes remain understudied. Here, we present a comprehensive genome resource for wild ovine species, landraces and improved breeds of domestic sheep, comprising high-coverage (16.10) whole genomes of 810 samples from 7 wild species and 158 diverse domestic populations. We detected, in total, 121.2 million single nucleotide polymorphisms, 61 million of which are novel. Some display significant (P < 0.001) differences in frequency between wild and domestic species, or are private to continent-wide or individual sheep populations. Retained or introgressed wild gene variants in domestic populations have contributed to local adaptation, such as the variation in the HBB associated with plateau adaptation. We identified novel and previously reported targets of selection on mor- phological and agronomic traits such as stature, horn, tail configuration, and wool fineness. We explored the genetic basis of wool fineness and unveiled a novel mutation (chr25: T7,068,586C) in the 30 -UTR of IRF2BP2 as plausible causal variant for fleece fiber diameter. We reconstructed prehistorical migrations from the Near Eastern domestication center to South-and-Southeast Asia and found two main waves of migrations across the Eurasian Steppe and the Iranian Plateau in the Early and Late Bronze Ages. Our findings refine our understanding of genome variation as shaped by continental migrations, introgression, adaptation, and selection of sheep. Key words: whole-genome sequences, genetic diversity, adaptive introgression, genetic selection, agronomic traits, migration. Introduction 2017; Hu et al. 2019; Ciani et al. 2020; Cao et al. 2021). Thus, Since their domestication, sheep (Ovis aries) have been pro- sheep have acquired a worldwide distribution through adap- viding meat, wool, skin, milk, and other products to humans. tation to a diverse range of environments and genetic im- They are essential for welfare of hundreds of millions of peo- provement under different production systems, developing ple living in rural communities of some Asian, African, and into plenty of unique breeds (Lv et al. 2015; Cao et al. 2021). South American countries (Nkedianye et al. 2019). Since their Hence, the extensive variations in landraces and improved domestication from Asiatic mouflon 10,000 years ago in the breeds underpin the genomic variations, adaptive character- Near East (Zeder 2008), sheep expanded to other parts of Asia istics, and important agronomic traits of domestic sheep (Lv et al. 2015; Zhao et al. 2017) and to Europe (Ryder 1984; (Yang et al. 2016; Li et al. 2020). Tapio et al. 2006; Ciani et al. 2020), Africa (Bradford and In this study, we sequenced genomes with an average Fitzhugh 1984; Muigai and Hanotte 2013), and the New depth of 16.7 for one Asiatic mouflon and a diverse set World (Crispim et al. 2012; Campos et al. 2020; Revelo et al. of 165 samples from 82 domestic breeds. These include 2020; Thorne et al. 2021). An early and major event was the African, South Asian, Southeast Asian, Central Asian, East post-domestication transition of hair sheep to wool sheep Asian, South American, European, and the Middle Eastern (Chessa et al. 2009). In addition, introgression from wild rel- breeds which were underrepresented in early work (Alberto atives into domestic populations has occurred (Barbato et al. et al. 2018; Naval-Sanchez et al. 2018; Deng et al. 2020; Li et al. 2 Genetic Diversity, Introgression, and Agronomically Important Loci of Worldwide Sheep . doi:10.1093/molbev/msab353 MBE 2020), but are essential to characterize the genetic variation identified varies from 9.82 to 12.9 million per domestic indi- and the introgression from wild into domestic sheep, and to vidual, 9.9 million (European mouflon) to 18.46 million explore for the first time the history of South-and-Southeast (Bighorn) per wild individual (supplementary fig. S1B, Asian sheep. By combining these novel genomic resources Supplementary Material online), and 10.4 to 13.0 million with previously published whole-genome sequences (WGS; per domestic breed, which did not depend on the population >10), we achieved a total of 810 samples from 158 domes- sample size (supplementary fig. S1C, Supplementary Material tic breeds and seven wild ovine species (fig. 1A, supplemen- online). tary table S1, Supplementary Material online). Using this We identified 61.0 million SNPs (50.2%) that have not comprehensive data set, we explored patterns of genome- been reported previously in the dbSNP (ftp://ftp.ncbi.nih.gov/ wide sequence variation among species/breeds, characterized snp/organisms/archive/sheep_9940/, last accessed July 10, regional and population-specific variants, inferred demo- 2021) or the European Variation Archive (EVA, including graphic expansions, discovered retained or introgressed var- accessions PRJEB33111, PRJEB23437, PRJEB6025, PRJEB15642, iants from wild sheep, and identified genes and variants and PRJEB14685) (fig. 1C), which could be due to prior un- associated with agronomically important traits in domestic derrepresentation of the sheep breeds studied here. We ob- sheep. served 4.51 million (3.7% of the total SNPs) SNPs shared among all the eight ovine species, 5.64 million (European mouflon vs. snow sheep) to 34.6 million (domestic sheep vs. Results Asiatic mouflon) between pairwise species, 50.7 million Whole-Genome Sequencing and Variation only in wild species, and 40.96 million SNPs shared between Individual genomes of one Asiatic mouflon and 165 samples landraces and improved breeds of domestic sheep (fig. 1D). from 82 domestic sheep breeds across the world were se- We identified 0.024 million (Brazilian Creole) to 1.74 mil- quenced to an average depth of 16.7 and average genome lion (local Moroccan sheep) breed-specific SNPs (supplemen- coverage of 96.40% to the Oar_rambouillet_v1.0 reference tary fig. S2, Supplementary Material online). The number of genome. These data were combined with 644 available SNPs present in each major geographic region, including genomes (fig. 1A and supplementary table S1, region-specific SNPs, varies from 27.8 million in the Supplementary Material online). Of the 82 domestic breeds, American breeds to 49.4 million in the Central-and-East whole genomes of 24 breeds from regions such as Africa (e.g., Asian breeds (supplementary fig. S3 and table S3, Menz), Central-and-East Asia (e.g., Minxian Black Fur, Gala, Supplementary Material online), which could be partially as- Azerbaijian Mountain Merino, Kazakh Arkhar-Merino, and cribed to the differences in individual WGS variability and Kazakh Edilbai), South-and-Southeast Asia (e.g., Indonesia number of samples included in the study. Thin-tailed, Kajli, Vembur, Mecheri, Madras Red, Pattanam, Kilakarsal, and Baruwal), Europe (e.g., Russian Edilbai, Lleyn, Variant Accuracy Hebridean, Ryelands, Swedish Finewool, Icelandic Leader, On average, 9,954,398 of the SNPs (82.05%) identified in do- Icelandic, and Oxford Down), and America (e.g., Morada mestic sheep and 8,252,754–13,976,424 of the SNPs detected Nova) were sequenced for the first time and extended our in the wild species (60.26–81.47%) were confirmed in the coverage of the sheep genetic resources (supplementary table sheep dbSNP database v.151 (supplementary table S4, S1, Supplementary Material online). Supplementary Material online). The proportions of SNPs In total, we analyzed 31.5 Tb raw data from 810 samples, confirmed were comparable with those reported in a previ- comprising 41 billion 100-bp and 150-bp paired-end reads ous study (Li et al. 2020). In addition, we inspected 133 se- with average genome depth of 16.10 per individual (or lected SNPs in candidate functional genes below from 2 to 16 39.8 Gb raw data per sample). The reads were mapped, individuals of ten populations obtained by Sanger sequencing resulting in varied alignment rate, genome coverage and se- approach, giving an overall validation rate of 98.58% (supple- quence depth among individuals, and across the whole ge- mentary table S5, Supplementary Material online). Overall, nome (fig. 1B and supplementary table S1, Supplementary the results provided confidence on the accuracy of variant Material online). After filtering, we identified 137.7 million calling. genome-wide genetic variants, including 121.2 million sin- gle nucleotide polymorphisms (SNPs) (12.1 million/individ- Genome Variant Map ual), 16.4 million small insertions and deletions (INDELs, The average SNP densities were 43.67/kb, 33.15/kb, and 19.20/ <50 bp, 15.5 million/individual), and 74,173 structural var- kb in autosomes, X, and Y-chromosomes, respectively, illus- iants (SVs, 51–886,322 bp, 9,190/individual; 58,256 deletions, trating the notion that the nucleotide diversity was lower on 10,890 translocations, 3,501 duplications, and 1,524 inver- the nonrecombining sex chromosomes (Sayres 2018). From sions) across all the samples (fig. 1B and table 1). The number all SNPs, 44.21% were located in or near genes, including of SNP variants per individual showed significant (Mann– upstream (6.71%) and downstream (2.66%) of open reading Whitney, P < 0.001) difference between species (supplemen- frame, introns (32.84%) and untranslated regions (UTRs, tary table S2, Supplementary Material online). Of the set of 0.88%). In total, 1,042,367 SNPs were located in protein- 121.2 million high-quality SNPs, 77.9% had a minor allele coding exons, causing 492,609 nonsynonymous and 549,128 frequency (MAF) < 5% (supplementary fig. S1A, synonymous mutations, with 45,679 start codon and 8,803 Supplementary Material online). The number of SNP stop codon changes, showing their potentially large impact 3 Lv et al. . doi:10.1093/molbev/msab353 MBE A B C D FIG. 1. Geographic distributions and genetic variants of wild and domestic sheep. (A) Sampling locations of 158 domestic sheep populations, and geographic distribution ranges and sampling locations of the seven wild sheep species. The brown and black dots indicate the sampling locations of domestic and wild sheep. The geographic distribution range of each wild sheep species is shown on the map with different colors as follows: Ovis musimon (yellow), O. orientalis (blue), O. vignei (green), O. ammon (orange), O. nivicola (purple), O. dalli (pink), and O canadensis (red). All European mouflon populations in wild parks of Europe have been imported recently from Coasica and Sardinia. (B) Diagram of SNPs found by resequencing of 734 individuals. Circles represent from outermost to innermost, 27 chromosomes (chr. 1–26, X) denoted by different colors (a), average sequence depth for 810 individuals (b), candidate genes selection based on the top 1% values of global FST (c), SNP abundance bars in improved populations (d) and landraces (e), nucleotide diversity (p) abundance bars in improved populations (f) and landraces (g), and INDEL abundance bars in improved populations (h) and landraces (i). (C) Venn diagrams for novel variants detected in wild and domestic sheep. (D) Venn diagrams for variants shared among wild sheep species and improved populations and landraces of domestic sheep. on function of the relevant genes (supplementary table S6, (e.g., BMPR1B and MYPN; supplementary table S7, Supplementary Material online). Supplementary Material online). Gene ontology (GO) and From the 26,342 SVs discovered only in domestic sheep, Kyoto Encyclopedia of Genes and Genomes (KEGG) annota- 11,391 (43.24%) were predicted to have functional conse- tions of the SVs indicated these genes to be involved in quences, and 6,816 (25.88%) are in 3,592 functional genes axonogenesis and axon guidance (supplementary fig. S4 and 4 Genetic Diversity, Introgression, and Agronomically Important Loci of Worldwide Sheep . doi:10.1093/molbev/msab353 MBE Table 1. Summary of Whole-Genome Variations Identified in Wild and Domestic Sheep. Variations Wild Sheep Domestic Sheep Total SNPs 90,771,286 70,626,497 121,253,260 INDELs 12,304,629 10,694,323 16,409,609 Structural variants Deletions 42,048 40,885 58,256 Translocations 4,823 9,315 10,890 Duplications 1,875 2,742 3,501 Inversions 1,035 1,064 1,524 Insertions 0 2 2 table S8, Supplementary Material online). We found 19,875 (Oar v.4.0 or Oar_rambouillet_v1.0) used in the SNP calling, SVs present only in the seven wild sheep species. GO and etc. KEGG enrichment analyses for 8,215 genes associated with Linkage disequilibrium (LD) analysis indicated that the these unique SVs in wild sheep uncovered their essential roles physical distance between SNPs measured as half of its max- in MAPK signaling pathway, cell morphogenesis, and GTPase imal value occurred at 40.0 kb (r2 ¼ 0.35) for domestic sheep activator activity (supplementary fig. S5 and table S9, and at 29.1–37.7 kb (r2 ¼ 0.31–0.40) for wild sheep except the Supplementary Material online). inbred European mouflons (fig. 2D). For landraces, LD values ranged from 14 to 296 kb (but >1,000 kb for the Cameroon Genomic Diversity population) and were higher than for improved breeds (14.6– Nucleotide diversities (p) in the eight Ovine species were 67 kb) (supplementary fig. S7A and B, Supplementary estimated at individual level after the correction for sample Material online). The LD decay in domestic sheep was slower size. Asiatic mouflon and urial have the highest nucleotide than in cattle (2–10 kb) (Mei et al. 2018), goat (1–2.5 kb) diversity, whereas the lowest was observed in bighorn sheep (Zheng et al. 2020), pig ( 8–20 kb) (Li et al. 2013), and dog (fig. 2A). Asiatic mouflon and urial, which genetically are not (1–3 kb) (Wang et al. 2015). completely separated (Demirci et al. 2013; Deng et al. 2020) Global FST value for domestic sheep (FST ¼ 0.082, (see also Results), are potential resources to expand the ge- P < 0.001) is lower than those across other domestic species netic resources of domestic improvement. The diversity of such as goats (FST ¼ 0.112) (Kim et al. 2019) and pigs (FST ¼ domestic sheep varies considerably with p values ranging 0.115) (Mu~ noz et al. 2019), which may partially reflect a rel- from p ¼ 3.1  10-3 (Bashibai breed) to 1.7  10-3 (Sv€ardsjö atively frequent genetic exchange during the development of breed; fig. 2B and supplementary table S10, Supplementary modern breeds (Kijas et al. 2012). FST value for breeds from Material online). There are no large differences among breeds the same region ranges from 0.050 (the Middle East) to 0.193 from different geographic regions (fig. 2C). Breeds with a low (America; supplementary fig. S8, Supplementary Material diversity show a high coverage by runs of homozygosity online). (ROH; supplementary fig. S6 and table S10, Supplementary Material online). Sheep from central and northern Europe Effective Population Size with a long history of breed formation and selective breeding Pairwise sequentially Markovian coalescent (PSMC) analysis have a relatively low diversity compared to most other breeds. revealed that all wild and domestic sheep species experienced There are no systematic differences between landraces and a dramatic contraction in Ne 100–400 thousand years ago productive breeds: several fine-wool breeds and breeds from (ka). Ancestors of domestic sheep (i.e., Asiatic mouflon) marginal regions such as the Balkan, the Caucasus, Libya, underwent a second decline in Ne during 10,000–30,000 Kazakhstan, and Tibet and Xinjiang of China have retained years BP, coinciding with the glacial periods (supplementary high diversity, whereas small population sizes have decreased fig. S9, Supplementary Material online). SMCþþ analysis the diversity in breeds from both categories (fig. 2B and sup- showed a continual decline of all the domestic sheep popu- plementary table S10, Supplementary Material online). In lations during 8,000–10,000 years BP, which probably general, the values of nucleotide diversity calculated here reflects the bottleneck effect of domestication (fig. 2E). are in the range of earlier estimates in domestic sheep Their subsequent population expansions show the expansion (p ¼ 1.9–2.5  10-3, Yang et al. 2016; p ¼ 1.6  10-3, Naval- of agricultural in sedentary societies (Gignoux et al. 2011). Sanchez et al. 2018; p ¼ 2.15–2.68  10-3, Alberto et al. 2018; Estimates of recent Ne by SNeP for the wild relatives and p ¼ 2.44–2.84  10-3, Pan et al. 2018; p ¼ 1.2–42  10-3, Hu domestic sheep breeds 50 generations ago were 38–54 et al. 2019; p ¼ 1.05–1.18  10-3, Li et al. 2020; p ¼ 3.2  10-3, and 42–53, respectively. Relatively smaller Ne was observed in Chen et al. 2021) (supplementary table S11, Supplementary South-and-Southeast Asian and American domestic popula- Material online). For the same populations such as Merino, tions, which may be due to strong artificial selection and/or Tan, Lop, and Tibetan sheep, most showed similar estimates genetic isolation (Li et al. 2020). in this and early studies (supplementary table S12, Supplementary Material online), but the discrepancy in Phylogenetic Relationships and Population Structure some populations (e.g., Iran sheep) can be due to the differ- Phylogenetic relationships among the wild and domestic ence in sequencing depth, sample size (e.g., the number sheep populations were inferred from whole-genome SNPs of animals sequenced) and the sheep reference genomes for 714 unrelated individuals (supplementary fig. S10 and 5 Lv et al. . doi:10.1093/molbev/msab353 MBE A B 0.008 50 0.004 π Latitude π 0 0.0028 0.000 0.0024 As mo p n n l li p or eep p ria ga ee lo flo ee ee 0.0020 uf U Ar ou sh sh sh sh m -50 tic ow rn n an tic ho es Sn gh ia pe in om Bi Th -100 0 100 200 ro D Eu Longitude C D E Predomestication (10,000-12,000 BP) Domestic sheep Argali Domestication episode (8,000-10,000 BP) European mouflon Snow sheep Asiatic mouflon Thinhorn sheep 105 0.008 Urial 0.8 Bighorn sheep π 0.6 Ne 104 0.004 African populations r2 American populations Central-and-East Asian 0.4 populations European populations the Middle Eastern populations 103 0.000 South-and-Southeast Asian populations 0.2 1 2 3 4 5 6 0 100 200 300 400 500 10 3 10 4 10 5 1. Central-and-East Asian populations 4. African populations Distance(Kb) Years (generation time = 3 years, mutation rate=1.51×10-8) 2. South-and-Southeast Asian populations 5. European populations 3. the Middle Eastern populations 6. American populations FIG. 2. Genetic variability of wild and domestic sheep. (A) Average nucleotide diversity (p) of eight wild and domestic sheep species. (B) Nucleotide diversity (p) of 158 domestic populations across the world. (C) Average nucleotide diversity (p) of six geographic groups of domestic sheep populations. (D) Decay of LD for eight wild and domestic sheep species, with one line per species. (E) Effective population sizes (Ne) for six domestic sheep populations inferred using SMCþþ, with one line per population. A generation time of 3 years and mutation rate of 1.51  108 per site per generation were used to convert coalescent scaling to calendar time. table S13, Supplementary Material online). The phylogenetic and-West African and East-and-South African breeds (Muigai tree (fig. 3A) shows a close relationship between the European and Hanotte 2013). mouflon and domestic sheep, followed by Asiatic mouflon. The geographic pattern is further supported by principal This supports the notion that the modern European mouflon component analysis (PCA) and model-based clustering via is the feral descendent of early domestic sheep rather than a the Admixture. A PCA plot of 653 unrelated domestic sheep genuine species of wild sheep (fig. 3A) (Ciani et al. 2020). shows largely separates the six geographic groups, most of Domestic sheep were split into six geographically structured which could be connected to continental origins (fig. 3B and clades (Central-and-East Asian, South-and-Southeast Asian, supplementary fig. S11, Supplementary Material online). The the Middle Eastern, African, European, and American popu- inbred Ouessant and Gotland breeds have the extreme posi- lations) with obvious explanations of apparent exceptions: tions. European and African populations are largely separated the East European Baidarak sheep (blue) cluster within from Asian populations with overlaps such as Merino-derived Central-and-East Asian sheep; thin-tailed Tibetan and individuals (fig. 3B). Southwest Chinese sheep (red) are close to Indian breeds Genetic clustering analysis using sNMF (Frichot et al. 2014) (yellow); Azerbaijan sheep (red) cluster within the Middle recapitulates the same patterns observed in the phylogenetic East sheep; and Asian fine-wool breeds (red) cluster are close tree (fig. 3A) and PCA (fig. 3B). At the optimal number K ¼ 11 to European Merino sheep (fig. 3A). The phylogenetic tree with the smallest cross validation error (supplementary fig. using the wild relatives as the outgroup revealed clear S12, Supplementary Material online), we observed five genetic between-population genetic differentiation in domestic clusters of the seven wild sheep species (Asiatic mouflon and sheep as most of the samples from the same breeds clustered urial are in the same genetic clusters; bighorn and thinhorn together. The Middle Eastern breeds were placed near the sheep are in one genetic cluster) and the six geographically root within domestic sheep, being the descendant of the first distributed genomic components of domestic sheep (fig. 3C). domesticated sheep there (Zeder 2008). The European Interestingly, sNMF analysis revealed the genetic components branch confirms separate subclusters for British and Nordic of wild relatives in domestic breeds, for example, the genetic breeds, respectively (supplementary fig. S11, Supplementary ancestry of Asiatic mouflon in Awassi and Bangladeshi sheep Material online) (Kijas et al. 2012; Ciani et al. 2020), and the breeds, argali in Bashibai sheep, and European mouflon in African branch indicates a genetic division between North- Cameroon and Sv€ardsjö sheep (fig. 3C). Additionally, we 6 Genetic Diversity, Introgression, and Agronomically Important Loci of Worldwide Sheep . doi:10.1093/molbev/msab353 MBE A 4 /OULE/LOCA North African d-South African 5 RD/LOCA East-an BENI/DM D/TIMA/LOCA West African DMAN/SARD/LOCA /LBS/RAH AN/SARD BOUJ/DMAN/SA S P/FSS/GT DMAN/LOCA DMAN/LOCA DMAN/LOCA DMAN/LOCA TIMA/LOCA DMAN/SAR 3 OSS/BAR SAH/MOS AWD/DPS WAD/MB NQA/RDA DMAN MOS/DJI DMAN DMAN LOCA LOCA LOCA CAM UDA ILS/OS MAK/M A YAN GOT SFS ME MS HEL FIN SOM SVS AFS BOG MEN SHE NRS AFH EFR NDZ /DA BENT HEB WS/W OLS DRS S AWA IROA KRS OUE GALA RYE AW OH ME OXD BOZ LIE MAZ TEX ZWA GHE ROS SWA /N /AZ SHA S AFH 2 HAS WH GTS GSS AB BAL KAR HAM SFK ALT /WBS DE G J KAS 1 N C/Z IB TH GC ITT SUM M ZR KU ZLX S GA BA SU AU TCS A/K L RH /GU DQ S SO R/B S BC /ME TO/C BB n AJ/ KAC C VA M S PR S R LE A MK flo GE KA 1 WG PIS L A R S MS E/K G DA R RIP ou JLS NE 6 XIS M MR S PA F ME S AR CA S S S/N IL LA S NE S S m CH C R U E VE 1 SAOM K DE SB MN L ic M KA AJ2 MX KA R B G B LS ED iat Z/A LS S/H SIS R ZS BE C P As LO S /G DL G BA 1 MA KH BA R /S DA AY O/T WZ J AB B M A S /S AJ HU Q TA HA /C S N TH S SS S W XW D S W HD 1. Central-and-East Asian populations 70.0 2. South-and-Southeast Asian populations 3. the Middle Eastern populations 4. African populations 5. European populations 6. American populations B Central-and-East Asian populations 0.10 South-and-Southeast Asian populations the Middle Eastern populations African populations European populations American populations PC2 (17.2%) 0.05 Sn ow sh ee p Bighorn sheep Thinhorn sheep 0.00 Argali 0.04 rial u flon/U Ouessant c mo Asiati lo n ouf -0.05 m an ope Eur Gotland -0.05 0.00 0.05 PC1 (22.6%) C North African East-and-South African 1.0 Ancestry 0 .0 West African 1 2 3 4 5 67 8 1. Central-and-East Asian populations 4. African populations 7. European mouflon 10. Snow sheep 2. South-and-Southeast Asian populations 5. European populations 8. Asiatic mouflon/Urial 11. Thinhorn sheep/Bighorn sheep 3. the Middle Eastern populations 6. American populations 9. Argali D E TEU 9.25 ka Ghost TAF (9.14-9.33) 8.74 ka (8.66-8.83) 3 TME 7.75 ka 2 10 (7.68-7.83) TCEA 7.35 ka NAME 1 14 11 (7.27-7.43) 7 12 5 13 TSSA2 5.93 ka 1 The Mediterranean routes (5.85-6.01) 2.12 9 (2.08-2.16) 2 The Danubian route 9 6 4 4.95 3 The northern Europe route (4.86-5.04) 1.68 4 The ancient sea trade route to the Indian subcontinent (1.67-1.69) 8 5 Thin-tailed sheep routes in Africa 1.19 (1.18-1.20) 6 Fat-tailed sheep routes in Africa aE=0.30 TSSA1 4.05 ka (0.293-0.316) 7 The first route to America ~15th century 0.302 2.42 2.80 1.83 (3.93-4.18) (0.296-0.307) (2.37-2.46) (2.78-2.83) (1.47-2.20) 8 The second route to America ~ 16th to 19th century 9 The recent route to America 8.87 0.268 3.23 2.63 (8.70-9.05) (0.263-0.274) (3.19-3.27) (2.35-2.91) 10 The Eastern Eurasia route Present Asiatic mouflon 1 2 3 4 5 11 The Qinghai-Tibetan Plateau route The major post-domestication migratory events of sheep All arrows are approximate 12 The Yunnan-Kweichow Plateau route The approximate major migratory events of sheep 13 The South-and-Southeast Asian route in South-and-Southeast Asia inferred in this study 1. Central-and-East Asian populations 4. African populations 14 The admixture event from Middle Eastern sheep The approximate domestication center into South-and-Southeast Asian sheep The geographic distribution of Asiatic mouflon 2. South-and-Southeast Asian populations 5. European populations 3. the Middle Eastern populations FIG. 3. Population genetic structure of wild and domestic sheep and the demographic history. (A) Neighbor-joining (NJ) trees of wild and domestic sheep based on whole-genome SNPs using the identity-by-state (IBS) genetic distances implemented in PLINK v1.90 (Purcell et al. 2007). BAJ1 and BAJ2 are from the population BAJ, and DEG1 and DEG2 are from the population DEG. (B) PCA of domestic sheep populations. (C) Population genetic structure of wild and domestic sheep inferred from the sNMF analyses (K ¼ 11) using whole-genome SNP data. (D) Demographic history reconstruction of domestic sheep. The best-supported demographic model of South-and-Southeast Asian population: South-and-Southeast sheep diverged from Central-and-East Asian population 5.93 ka and was later admixed by the Middle Eastern population 2.42 ka. The Middle Eastern population contributed 30% of their gene pool to the South-and-Southeast Asian population. Dates of events (ka) are indicated on the 7 Lv et al. . doi:10.1093/molbev/msab353 MBE noted clear evidence of genetic heterozygosity among First, sequences from Asiatic and European mouflons that European, African, and Asian populations. This indicated sub- were inherited directly by the first domestic sheep may stantial historical genetic flows among elite genetic stocks have been retained in specific breeds, whereas other parts from different continents in the development of local breeds of the genome adapted to the domestic status and breeding (Galal et al. 2008; Ciani et al. 2020). objectives, for instance during the transition from hair to One objective of this study was to fill the gaps of the global wool sheep (Chessa et al. 2009). Second, introgression may genomic coverage of population-scale WGS data by including have taken place during later contacts between wild and samples from the so far neglected South-and-Southeast domestic sheep. Here we investigated the pattern and con- Asian, African, and South American sheep breeds. South- sequences of the retained or introgressed sequences (RIS) and-Southeast Asian breeds cluster with breeds from using a combination of D and fd statistics. Pakistan and India. South American breeds share ancestry We used ABBA-BABA statistics to detect ancestry shared by with sheep from both Western Africa and Europe as reported wild donors to domestic breeds. In agreement with Ciani et al. previously (Kijas et al. 2012). (2020), we found substantial influence (negative D values) of European mouflon on several breeds with the higher jZj values Genetic Origins of South-and-Southeast Asian Sheep in North European sheep such as Welsh Mountain, Icelandic In order to explore the origin of South-and-Southeast Asian Leader, and Shetland (fig. 4A and supplementary table S17, sheep (Singh et al. 2013; Lv et al. 2015), we tested three alter- Supplementary Material online). Retained or introgressed an- native models (supplementary fig. S13, Supplementary Material cestry from Asiatic mouflon and the related urial are much online). We computed the expected multidimensional site fre- lower (D value is closer to zero), although Z values indicate quency spectra (SFS) under specific models and compared it statistical significance for a few breeds. Asiatic mouflon and with the observed multidimensional SFS by a composite like- urial appear to have influenced the same breeds, confirming lihood method. The best supported SFS with the highest aver- that these wild species have become mutually admixed (fig. 3C) age value of estimated log10(likelihood) (supplementary fig. S14 (Demirci et al. 2013; Deng et al. 2020; Chen et al. 2021). and table S14, Supplementary Material online) suggested that Interestingly, the genomes of Bashibai sheep carry a sizable South-and-Southeast Asian sheep descended from Central- argali component. In addition to several other potential explan- and-East Asian population 5.85 to 6.01 ka (95% CIs). The op- ations such as horizontal gene transfer, hybridization, and com- timal model also showed a later admixture from the Middle mon ancestry (Wang et al. 2021), the observation could be very Eastern sheep 3.93–4.18 ka, which contributed approximately likely the result of introgression (fig. 4A and supplementary 29–32% of the gene pool (fig. 3D and supplementary table S15, table S17, Supplementary Material online). Supplementary Material online). To further locate the RIS genomic regions, we computed The best model further indicated a divergence of the the fd values for 100-kb windows and 50-kb step across the European, African, and the Middle Eastern sheep from their genomes of relevant domestic populations. After the P-values common ancestor 9.16–9.33, 8.66–8.83, and 7.68–7.83 ka, re- adjusted (false discovery rate, FDR), we detected 258 (fd > spectively (fig. 3D and supplementary table S15, 0.1714)  434 (fd > 0.1585) significant introgressed blocks Supplementary Material online). This was in agreement from argali, 283 (fd > 0.2271)  335 (fd > 0.2363) significant with the initial expansions of ancestral sheep as evidenced introgressed tracts from urail, 311 (fd > 0.2374)  344 (fd > by the archaeological findings and early genetic studies (Zeder 0.5322) significant introgressed blocks from Asiatic mouflon, 2008; Deng et al. 2020). Asymmetric gene flow occurred from 296 (fd > 0.5117)  323 (fd > 0.3157) significant introgressed the Central-and-East Asia to the South-and-Southeast Asia tracts from European mouflon, respectively (supplementary (8.70–9.05 migrants per generation) and from the South-and- fig. S15 and table S18, Supplementary Material online). The Southeast Asia to the Middle East (4.86–5.04 migrants per proportion of genome introgression (PGI) of domestic sheep generation) (fig. 3D and E and supplementary table S16, breeds from different donor wild species on the basis of the fd Supplementary Material online). statistics gives different patterns (fig. 4B–E and supplementary table S19, Supplementary Material online). The proportion of Retained or Introgressed Sequences from Wild Sheep RIS segments, which reflects the actual genome size of do- in the Domestic Genomes mestic sheep retained or introgressed from wild relatives was The sharing of ancestry by wild and domestic sheep, indicated estimated to be 7.64–13.40%, 5.94–9.61%, 4.31–6.39%, and by the Admixture analysis, may be the result of two scenarios. 3.67–5.03% of their whole genomes originating from left and right, and estimated migration rates (migrants per generation) are scaled by the effective population size were given above/below the corresponding arrows. The 95% confidence intervals (CIs) are shown in parentheses. (E) Major sheep migrations across the world. (1) The Mediterranean routes (Ryder 1984), (2) the Danubian route (Ryder 1984), (3) the northern Europe route (Tapio et al. 2006), (4) the ancient sea trade route to the Indian subcontinent (Singh et al. 2013), (5) thin-tailed sheep routes in Africa (Muigai and Hanotte 2013), (6) fat-tailed sheep routes in Africa (Muigai and Hanotte 2013), (7) the first route to America 15th century (Dunmire 2013), (8) the second route to America 16th to 19th century (Spangler et al. 2017), (9) the recent route to America (Ryder 1984), (10) the Eastern Eurasia route (Lv et al. 2015), (11) the Qinghai- Tibet Plateau route (Zhao et al. 2017; Hu et al. 2019), (12) the Yunnan-Kweichow Plateau route (Zhao et al. 2017; Hu et al. 2019), (13) the South- and-Southeastern Asia route (this study), and (14) the admixture event from Middle Eastern sheep into South-and-Southeastern Asian sheep (this study). 8 Genetic Diversity, Introgression, and Agronomically Important Loci of Worldwide Sheep . doi:10.1093/molbev/msab353 MBE European mouflon Asiatic mouflon Urial Argali A 0 D(Menz, X; Wild sheep, Bighorn sheep) D-value -0.09 0 British North European South-and-East European Asian and African LOCA Z-value IROA GAR SWA BOUJ NDZ SUM RDA MER GUR NQA MER BGE GUR RDA SUM NQA WDS BGE RED KRS MNR KAED HUS SSS THO DAR BAG CHA BMN GAR ITT BOZ PRAM BAJ KAZ BAS MOH GUR TCS SUM SOM MAK BSI BER SIS ZLX VAL CAM MKS BENI YAN BAR ROM BCS DAS UDA LBS WAD NQA DPS OULE RDA AWD MBS DAME DEG ARME AZME DMAN DQS SVS BENT RHS KATO SAL SARD UNKN CHU TIMA Milk HEL CAUC LEC ZAB MOS RYE WBS SUMA GCN Meat BBS ALT PIS DJI SAH EFR ZWA LOCA FSS BSB SOL SFS GOT OLS WHSF ZRJ MFW NRS DRS NWS KUIB OUE OSP XIS FIN CAS HEB SFK LIE TEX SWA ICS MSF MER SHE OXD RIP HAM GTS ROS ILS WMS -30 B C PGI from European mouflon PGI from Asiatic mouflon 0.0764 - 0.0914 0.0594 - 0.0638 0.0915 - 0.1014 0.0639 - 0.0694 0.1015 - 0.1128 0.0695 - 0.0773 0.1129 - 0.1227 0.0774 - 0.0892 0.1228 - 0.1340 0.0893 - 0.0961 0 0 D E PGI from Urial 0.0431 0.0432 - 0.0477 PGI from Argali 0.0478 - 0.0486 0.0367 0.0487 - 0.0553 0.0368 - 0.0503 0.0554 - 0.0638 0 0 F 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 AKT3 HBB MRP4 Argali Urial 0 0.2 0.4 0 0.2 0.4 0.6 European mouflon Asiatic mouflon 0 0.2 0.4 0.6 0.8 0 0.2 0.4 0.6 NPFFR1 16 17 18 19 20 21 22 23 24 25 26 G 140 European mouflon 120 Asiatic mouflon Urial Length of accumulated Argali RIS fragments (Mb) 100 80 60 40 20 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 Chromosomes FIG. 4. Genetic introgression of wild sheep species in domestic populations. (A) D statistic tests in the form (Menz sheep, X, wild sheep, Bighorn sheep) where wild sheep indicates European mouflon, Asiatic mouflon, urial, or argali, respectively. X indicate a candidate domestic population. Introgressions from European mouflon into British, North European, Southwest-and-East European, and Asian, and African populations are indicated by green, gold, purple, and black color, respectively. (B, C, D, E) The PGI across the whole genome of each wild sheep species in worldwide domestic populations. The levels of PGI from European mouflon, Asiatic mouflon, urial, and argali are mapped based on the results of fd statistics: fd (Menz sheep, X, European mouflon, Bighorn sheep), fd (Menz sheep, X, Asiatic mouflon, Bighorn sheep), fd (Menz sheep, X, Urial, Bighorn sheep), and fd (Menz sheep, X, Argali, Bighorn sheep), separately. X indicates a candidate domestic population. PGI can be estimated with the equation: PGI =(RfdiGi)/G, where fdi and Gi refer to the fd value and the window size in base pairs for the ith window, and G is the genome size in base pairs (Zhou et al. 2020). (F) RIS of wild sheep species in domestic sheep on each chromosome. Colored heatmaps of blue, red, purple, and gray show the abundance of genomic regions containing RIS from the four wild sheep species based on fd statistics (Argali, Asiatic mouflon, Urial, and European mouflon). The genes with the most significant signals of introgression are highlighted. (G) Length of accumulated RIS from the four wild sheep species (Argali, Asiatic mouflon, Urial, and European mouflon) in domestic sheep. 9 Lv et al. . doi:10.1093/molbev/msab353 MBE European mouflon, Asiatic mouflon, urial, and argali, respec- 2017), a generation time of 3 years for domestic sheep (Zhao tively (supplementary table S19, Supplementary Material on- et al. 2017), a divergence time of 2.36 Ma between argali and line). Counting the RIS tracts of union set for every wild donor domestic sheep (Yang et al. 2016), 1.26 Ma between urial into all domestic recipients yielded 566, 1437, 1718, and 4808 and domestic sheep (Rezaei et al. 2010), 11 ka between RIS tracts from argali, urial, Asiatic mouflon, and European Asiatic mouflon and domestic sheep (Chessa et al. 2009), mouflon, respectively (fig. 4F and supplementary table S20, and 5–6 ka between European mouflon and domestic sheep Supplementary Material online). The accumulated genome (Vigne 1992). This gives a length of tracts in domestic sheep sizes of RIS segments on each chromosome in domestic shared with argali of 42.4 bp, with urial of 79.4 bp, with Asiatic genomes were shown in figure 4G. The number and size of mouflons of 9,091 bp, and with European mouflon of 16.7– detected introgression windows show different patterns 20 kb. The probability of a length of 85.4 kb (i.e., the observed among different chromosomes and different donor wild spe- introgressed region from argali, urial, Asiatic mouflon, and cies (fig. 4F and G and supplementary table S20, European mouflon) is zero, which ruled out that the intro- Supplementary Material online). gressed region identified in Changthangi sheep were due to Of the tracts shared with the wild relatives, 79 common random ILS without selection. blocks were shared with all four wild species and span 125 functional genes (supplementary tables S21 and S22, Novel Signatures of Selection Supplementary Material online). In the GO analysis of the The distribution of genetic differentiation statistic FST values retained or introgressed genes, we identified a significant (ad- for 50-kb windows with 25-kb step across the genomes of justed P < 0.05) overrepresentation of genes involved in orga- domestic sheep, an indicator of selection, is summarized in nonitrogen compound (GO:0010243), hemoglobin complex figure 6A and supplementary table S25, Supplementary (GO:0005833), olfactory receptor activity (GO:0004984), and Material online. Selective sweeps contained 559 genes, of oxygen transporter activity (GO:0005344) (supplementary fig. which 259 (46%) were novel, and 300 were identified by sev- S16 and table S23, Supplementary Material online). The eral previous investigations on the basis of WGS (supplemen- KEGG analysis identified a major and significant enrichment tary table S25, Supplementary Material online) (Yang et al. of genes associated with insulin signaling pathway and 2016; Alberto et al. 2018; Li et al. 2020). Functional annotation African trypanosomiasis (oas05143) (supplementary fig. S16, of putatively selected genes revealed that they were predom- Supplementary Material online). In particular, the significantly inantly associated with morphological and agronomic traits, enriched GO terms and KEGG pathways contained the oxy- including genes such as KRT71 and FGF7 (skin development gen transporter activity, patterning of blood vessels, and olfac- function), Olr226, Olr51B2, and Olr2D3-like (olfactory func- tory receptor activity genes such as HBB, VEGFC, and Olr51L1. tion), and HBB, HBB-like, HBBC, and HBE2 (oxygen transporter In order to illustrate the potential of this approach, we focus function, see above). In addition, the previously characterized on the Changthangi sheep in the high-altitude region of Ladakh genes with relevant traits such as RXFP2 (see also fig. 5B), ca. 4,000 m above sea level. Genome scans of fd values of PDGFD (Dong et al. 2020), and TMEM154 (Heaton et al. Changthangi versus Asiatic mouflon were shown in figure 5A. 2012; Hu et al. 2019; Li et al. 2020) have elevated FST values One of the introgression signals is the HBB region on chromo- (fig. 6A). some 15 (chr15: 51.9–52.0 Mb). Various analyses such as the We found significant overlap of the FST outliers with population branch statistics (PBS, Yi et al. 2010; fig. 5B and known QTLs in sheep (permutation test, P < 0.001), most supplementary table S24, Supplementary Material online), dxy, notably with milk-, meat-, and disease-resistance-related FST distances, and haplotype pattern (supplementary fig. S17, QTLs (supplementary table S26, Supplementary Material on- Supplementary Material online) suggested an anomalous phy- line). Within the protein-coding exons in the selective regions, logeny as the result of wild introgression of variants of HBB and/ we identified 8,753 nonsynonymous SNP mutations in 559 or neighboring genes that confer hypoxia adaptation in the putatively selected genes. The allele frequencies of nonsynon- Changthangi sheep (Hu et al. 2019). In the GO analysis of the ymous SNPs and their genotypes pattern in ten genes (e.g., introgressed genes from wild species into Changtangi sheep, we EVC, GDF6, RXFP2, TMEM154) showed differences among identified a significant (P < 0.05) overrepresentation of genes populations with varied phenotypes such as body size (e.g., involved in oxygen transporter activity (GO:0005344), oxygen normal/dwarf), horn status (e.g., horned/polled), and disease binding (GO:0019825), pigmentation (GO:0043473), cellular re- resistance (e.g., pneumonia susceptibility/resistance) (supple- sponse to UV-B (GO:0071493), and energy homeostasis mentary table S27, Supplementary Material online). (GO:0097009) (fig. 5C). The KEGG analysis identified a major A genome-wide reduction of diversity (ROD) as alterna- and significant enrichment of genes associated with olfactory tive selection signature indicated 134 genomic regions with transduction (oas04740) and PI3K-Akt signaling pathway top 1% of both ROD values between landraces and im- (oas04151) (fig. 5C). The two highest peaks in the PBS scan proved breeds of domestic sheep (fig. 6B and supplemen- (fig. 5B) cover IGF2BP2, encoding an important protein in cell tary table S28, Supplementary Material online). We found metabolism and development, and RXFP2, a gene involved in that ROD chr15: 49,000,001–50,000,000 were in concor- horn development (Li et al. 2018; Aldersey et al. 2020). Both also dance with high nonsynonymous to synonymous (N/S) have elevated FST values (fig. 6A). ratio (N/S ratio ¼ 1.49; supplementary table S29, We calculated the probability of incomplete lineage sort- Supplementary Material online). Additionally, 59 SVs under ing (ILS) using a recombination rate of 1.5 cM/Mb (Petit et al. selection during sheep improvement including 1,601 10 Genetic Diversity, Introgression, and Agronomically Important Loci of Worldwide Sheep . doi:10.1093/molbev/msab353 MBE FIG. 5. Signals of introgression and selection in the Changthangi sheep. (A) Manhattan plot showing the introgression signal from Asiatic mouflon to Changthangi sheep. The locations of the HBB, HBE2, HBBC, HBE1, RAB6A, and MRPL48 genes are shown by the highest fd value. (B) Signals of selection in the Changthangi sheep detected by the PBS (Yi et al. 2010). PBS values in windows of 100 kb and names of genes associated with the highest peaks are shown. (C) GO and KEGG pathway enrichment analysis based on genes across significant introgressed regions from wild species to Changthangi sheep. The dot size shows the gene count enriched in the pathway, and the dot color shows adjusted significance value of the enriched pathways. candidate genes (supplementary table S30, Supplementary P < 0.001) differences in a few nonsynonymous SNP variants Material online), several of which are related to regulation and genotype patterns between the populations with differ- of ovulation rate (e.g., BMP15 and CYP17) (Estienne et al. ent statures (supplementary fig. S18, Supplementary Material 2017). online). Several of the genes (e.g., IGF2BP1, SMOX family mem- The composite likelihood ratio (CLR) statistic scores for 10- ber, and FGFR family member) related to stature in sheep kb windows in Ouessant sheep with an extreme dwarf phe- have been associated with body size in cattle and dogs notype was shown in figure 6C. Outlier windows with CLR > (Rimbault et al. 2013; Hayward et al. 2016; Bouwman et al. 1.31 (P < 0.001) were considered to be under artificial selec- 2018). tion and cover 1,702 genes (supplementary table S31, Supplementary Material online). These genes were mainly Selective Signatures Associated with Wool Traits involved in bone resorption (e.g., PTGER4), bone mineraliza- Genome-wide XP-CLR selection tests of domestic breeds pos- tion (e.g., BMP6), and dwarfism (e.g., GDF6 and PTDSS1). In sessing different fleece fiber variations (e.g., hair, coarse wool, gene GDF6, we also detected significant (Mann–Whitney, medium wool, and fine wool; supplementary table S32, 11 Lv et al. . doi:10.1093/molbev/msab353 MBE 0.8 A 0.6 RXFP2 FOXO4 SSBP2 Global FST ABCC4 HOXB13 LOC101111733 MSRB3 0.4 LOC114118424 IRF2BP2 EVC KIF1C NFX1 SOCS2 NF1 TMEM154 CSMD1 HOXC13 PDGFD IGF2BP2 ASTN2 OSTF1 HOXA1 CDH8 KIF9 0.2 0.0 1 2 3 4 5 6 7 8 9 10 12 14 16 18 20 23 26 113,776,250 113,776,275 113,776,300 113,776,325 30,967,900 30,967,950 30,968,000 30,968,050 5,776,900 5,777,000 EVC CDS RXFP2 CDS TMEM154 CDS 113,776,294-bp(T/G;Tyr/Asp) 30,967,896−bp(G/A,Glu/Lys) 5,776,842−bp(A/G;Lys/Glu) B 0.3 MRPL17 ATPAF1 TPP1 TMEM275 TRIM3 POMGNT1 SPIRE1 CCKBR CHMP1B MAST2 BMP2 ACTR1B CAMK4 HEBP2 RXFP2 GPCPD1 KNCN ZNF804A ANKRD39 SMIM28 HSPH1 KIF5B TMEM151A 0.2 NAB1 EXOC6B TMEM176B FBXL5 PERP SLC25A10 TSSK1B KIF15 IRF2BP2 NEMP2 ROD CNOT8 STBD1 OSR2 DYRK3 MC1R DGAT2 MCTP2 MDGA1 FGF8 SMIM18 0.1 0.0 1 2 3 4 5 6 7 8 9 10 12 14 16 18 20 22 25 C 14 89,112,389-bp(C/A;Thr/Lys) GDF6 CDS OGFOD2 89,112,200 89,112,400 12 10 CHRM2 CLR score TMED7 8 PRSS16 SOX6 GDF6 ATP2A1 6 PTDSS1 FGF10 PRR14 SAMD5 ESR1 IGF2BP1 PTGER4 BMP6 TGFBR3 PDGFC SMOX PTPN5 4 INPP4B BMP5 XPR1 EPHA2 2 0 1 2 3 4 5 6 7 8 9 10 12 14 16 18 20 22 24 26 FIG. 6. Genome-wide selective signals across all the domestic sheep populations based on global FST and ROD, and associated with the stature (i.e., body size) phenotype. (A) Genome-wide distribution of global FST across all the 158 domestic sheep populations. (B) Candidate selective genomic regions based on ROD analysis. (C) Manhattan plot for the dwarf (i.e., small body size) phenotype showing the CLR score distribution in the Ouessant sheep population. Genes within regions with significant outlier scores are highlighted. Supplementary Material online) were shown in all the six PTCH1) genes. We also identified 414 selective SVs associated pairwise comparisons (fig. 7A and supplementary figs. S19– with wool traits in the selective sweeps of SVs (supplementary S23 and tables S33–S38, Supplementary Material online). tables S39–S44, Supplementary Material online). Annotation From 2,759 genes associated with 5,995 selective signals, 86 of the SVs indicated their colocalizations with genes or QTLs are scored in at least two of the six comparisons and with associated with mean fiber diameter (e.g., FST gene) (Ma et al. functions involved in fiber and wool development, including 2017). Positional gene enrichment analysis of the gene sets novel (e.g., HR, ASTN2, TWIST2, IRF4, INHBA, and FGF7) and showed significant (Padj. < 0.05) overrepresentation of genes previously reported (e.g., MCIR, FOXQ1, TP63, VEGFA, and related to fiber diameters in these regions (De Preter et al. 12 Genetic Diversity, Introgression, and Agronomically Important Loci of Worldwide Sheep . doi:10.1093/molbev/msab353 MBE A C B D E F FIG. 7. Genomic evolution of hair/coarse-wool/medium-wool/fine-wool sheep. (A) Whole-genome selective signals between fine-wool and hairy populations of domestic sheep by the cross-population composite likelihood ratio (XP-CLR) test. (B) Six comparisons of p ratio (pHair sheep/pCoarse- wool sheep, pHair sheep/pMedium-wool sheep, pHair sheep/pFine-wool sheep, pCoarse-wool sheep/pMedium-wool sheep, pCoarse-wool sheep/pFine-wool sheep and pMedium- wool sheep/pFine-wool sheep) indicates strong selective signal in IRF2BP2 gene. (C) Allele frequency of the mutation site and frequencies of genotypes CC, TT, and CT in 30 -UTR of IRF2BP2 (chr25: T7,068,586C), and significance of the genotype difference is tested. (D) Haplotype pattern of IRF2BP2 among wild sheep, hairy, coarse-wool, medium-wool, and fine-wool populations of domestic sheep. (E) Sequence comparison among different species at mutations chr25: 7,068,586, and chr25: 7,070,630. (F) GO and KEGG pathway enrichment analyses, with the significant (P < 0.05) GO terms and pathways and associated genes shown. 13 Lv et al. . doi:10.1093/molbev/msab353 MBE 2008). All candidate genes under selection are directly or in- miR-20a-3p reduced reporter activity compared with a neg- directly associated with wool growth and development ative control mimic (fig. 8D). However, mutated 30 -UTR of (fig. 7F). IRF2BP2 abolished the reduction of luciferase by the oar-miR- We further dissected the genomic architecture of these 20a-3p-mediated (fig. 8D and E and supplementary fig. S26 strong candidate genes by calculating allele frequencies at and table S46, Supplementary Material online). Altogether, nonsynonymous SNPs in exons. We observed significant dif- these data supported that the WT IRF2BP2 30 -UTR was a ferences in allele frequencies (Mann–Whiteney test, direct and functional target of oar-miR-20a-3p, which could P < 0.001), and genotype or haplotype patterns of nonsynon- directly bind IRF2BP2 and potentially regulate IRF2BP2 expres- ymous SNPs in IRF2BP2 and MC1R genes between breeds sion in the skin tissues of coarse-wool sheep. displaying different wool phenotypes (supplementary table We also detected a strong and consistent selective signal in S27 and fig. S24, Supplementary Material online). the gene VEGFA (vascular endothelial growth factor A) in four comparison tests (hair vs. fine wool, hair vs. medium wool, Functional Analysis of IRF2BP2 Variation coarse wool vs. fine wool, and coarse wool vs. medium wool, Remarkably, the signals overlapping with IRF2BP2 are outliers supplementary tables S26–S31, Supplementary Material on- in four of the six comparisons in the XP-CLR analysis and p line). Early evidence suggested that IRF2BP2 coactivates ratio analysis (fig. 7B and supplementary tables S26–S31, VEGFA expression, leading to the activation of secondary Supplementary Material online). We annotated the gene con- wool follicles (Ho et al. 2012; Li, Lu, et al. 2012; Li, Man, taining 85 SNPs and detected one intron and two exons in et al. 2012; Zhang et al. 2018). In the histological analysis, the gene. The 30 -UTR variant (chr25: 7,068,586) and the hematoxylin and eosin (H&E) staining showed a larger ratio intronic mutation (chr25: 7,070,630) showed significantly of diameters (DS/DP) but a smaller number (NS/NP) of sec- (P < 0.05) differentiated allele and genotype frequencies be- ondary to primary wool follicles in Wadi sheep (coarse wool) tween the groups of breeds with different fleece fiber diam- than in Chinese Merino sheep (fine wool) in the important eters, such as hair, coarse wool, medium wool, and fine wool developmental stage of wool follicles (fig. 8F), which conse- (fig. 7C and supplementary fig. S25, Supplementary Material quently accounted for the higher yields and finer wool of online). In domestic sheep, the ancestral T (chr25: 7,068,586) uniformly small-diameter fibers in fine-wool sheep and G alleles (chr25: 7,070,630) showed higher frequencies in (Galbraith 2010). Thus, oar-miR-20a-3p may act on the de- hair and coarse-wool populations, corresponding to the an- velopment of wool follicles in coarse-wool sheep through its cestral alleles, whereas the derived reference C (chr25: effect on IRF2BP2, which potentially inhibits VEGFA signaling 7,068,586) and A (chr25: 7,070,630) alleles were dominant in expression in wool follicles (fig. 8F). medium-wool and fine-wool populations (fig. 7C and D). Further characterization of alleles of these two variants in a Discussion large variety of species including wild sheep, wild yak, goat, In this study, we performed a whole-genome sequencing and horse showed both ancestral alleles in animals with hair analysis of 810 wild and domestic sheep. This is thus the or coarse-wool coat (fig. 7E). Thus, this observation indicated hitherto most comprehensive data set on population genetic that the two variants in IRF2BP2 may have advantageous structure of worldwide sheep using more than 2 times as effects on fleece fiber diameter. In our worldwide sheep panel, many individuals and WGS as in previous studies (Alberto these SNPs do not cosegregate with the EIF2S2 insertion up- et al. 2018; Naval-Sanchez et al. 2018; Deng et al. 2020; Li et al. stream of the IRF2BP2 gene previously found to be associated 2020). We explored the occurrence of sequences from wild with wool and fine-wool genotype (Demars et al. 2017). species in the domestic genomes and further identified useful The fragments per kilobase million (FPKM) value for nonsynonymous SNP mutations and SVs involved local ad- IRF2BP2 mRNA expression showed significant (P < 0.05) dif- aptation and agronomic traits. ference in skin tissues at different prenatal stages between coarse-wool and fine-wool sheep (fig. 8A and supplementary The Influence of Wild Species on Domestic Genomes table S45, Supplementary Material online). No expression was In general, a high level of diversity with novel, rare, and private observed in IRF2BP2 for coarse-wool sheep (fig. 8A). We iden- variants was observed in wild relatives compared with both tified binding sites of the highly conserved miRNAs oar-miR- improved and landrace breeds (fig. 1C and D). In domestic 496-3p, oar-miR-379-3p, oar-miR-411a-3p, and oar-miR-20a- sheep, we observed many region-specific variants and a clear 3p binding in the 30 -UTR of IRF2BP2, which suggested clustering of populations according to their geographic ori- IRF2BP2 as a target of these four miRNAs (supplementary gins (supplementary fig. S3, Supplementary Material online). table S46, Supplementary Material online). MiRNA expression Several valuable genetic variants for modern sheep appeared profiling in skin samples of three coarse-wool Hu sheep and to originate from wild sheep species. Genetic influence of wild nine Aohan fine-wool sheep (supplementary table S47, on domestic species as a mode of adaptation have already Supplementary Material online) revealed expression of only been described for chicken (Wang et al. 2020), pigs (Ai et al. oar-miR-20a-3p in coarse-wool sheep (fig. 8B). 2015), dogs (Miao et al. 2016), and cattle (Chen et al. 2018) Binding of oar-miR-20a-3p to the 30 -UTR of wild-type and is supported by ancient DNA studies (Frantz et al. 2020). (WT) IRF2BP2 (fig. 8C) was confirmed by luciferase assay in We observed that among the wild sheep species, the HEK293T cells. We introduced the WT IRF2BP2 30 -UTR into a European mouflon has the most widespread influence on luciferase expression vector and found that expression of oar- domestic sheep, although the influence of Asiatic mouflon 14 Genetic Diversity, Introgression, and Agronomically Important Loci of Worldwide Sheep . doi:10.1093/molbev/msab353 MBE A 100 100 100 B C Fine-wool sheep Coarse-wool sheep oar-miR-20a-3p 80 80 80 0.6 IRF2BP2 60 60 60 7,068,586 TPM FPKM 0.4 40 40 40 3’ -…UGAAAUUCACGAGUAUUACGUCA…-5’ 0.2 20 20 20 5’-…GUUGGUUACAUAAUACACA…-3’ 3’ UTR WT (coarse-wool sheep) 5’-…GUUGGUUACGUAAUACACA…-3’ 3’ UTR MU (fine-wool sheep) 0 0 0 0 r-m -3 3p 3p p p 2 ths 5 5 th p l s p ol d ee oo E1 0d th 5 m th -3 -3 m th E5 E7 Bir ar 20 oa -miR 0a- iR 79- E6 ee o on on 5 on 6 Bir r-m 11a 96 sh -w E9 ye sh e-w m -2 ne -4 -4 s iR 1 iR Fi ar r-m Co r oa oa oa D E F Wadi sheep (coarse-wool) primary follicle N.S. Coarse-wool * transverse luciferase activity 15 oar-miR-20a-3p secondary PF negative control (NC) VEGFA follicle Relative 10 oar-miR-20a-3p expression SF MU: fine-wool mRNA Coarse-wool sheep 40× 5 WT: coarse-wool mRNA oar-miR-20a-3p IRF2BP2 mRNA PF transverse 0 Fine-wool × hair SF IRF2BP2 NC U T U T 3’UTR W +W M +M IRF2BP2 follicle secondary p+ p+ NC mRNA follicle -3 -3 IRF2BP2 coactivates VEGFA 0a 0a VEGFA -2 -2 ① expression 40× iR iR r-m r-m Fine-wool sheep ② Merino sheep (fine-wool) oa oa FIG. 8. Histological analysis, mRNA, and miRNA expressions and genetic mechanisms for the gene IRF2BP2 and oar-miR-20a-3p involved in the fleece fiber variation in domestic sheep. (A) Expressions of IRF2BP2 in embryonic developmental stages and 2 years old of fine-wool (left), coarse- wool sheep (medium) and boxplot of FPKM distribution between fine-wool and coarse-wool sheep (right). (B) Expressions of miRNAs in coarse- wool sheep (Hu sheep) normalized by reads per million reads mapped to miRNAs (RPM). (C) Putative binding sites for oar-miR-20a-3p to the 30 - UTR of IRF2BP2. (D) Dual-luciferase reporter assays, *P < 0.05 (t-test); N.S., not significant. (E) Proposed genetic mechanisms for IRF2BP2 and oar- miR-20a-3p involved in the fleece fiber variation in domestic sheep. ‹Ho et al. (2012) and ›Zhang et al. (2018). (F) Histological transverse section of skin in Chinese Merino (fine wool, at embryonic 85 days) and Wadi sheep (coarse wool, at embryonic 90 days); PF, primary wool follicle; SF, secondary wool follicle. may have been underestimated by a substantial influence of discrepancy is likely explained by sensitive detection of short the urial on the current Asiatic mouflon population (this mouflon fragments by WGS, which are missed by 50K SNP study) (Demirci et al. 2013; Deng et al. 2020). The European BeadChip genotyping with an average spacing between SNPs mouflon is supposed to have emerged as a feral population of about 54 kb. Introgression is the most plausible explanation when the original hair-coat domestic sheep were replaced by for the argali influence on the Bashibai sheep. We observed wool sheep. During this process, mouflon DNA may have that domestic sheep may have acquired advantageous alleles been retained in domestic sheep, in particular in North of several immune and sensory genes via natural or managed European breeds, some of which have retained a wild appear- hybridization with the wild relatives (supplementary table ance (Chessa et al. 2009; Ciani et al. 2020; Deng et al. 2020). S22, Supplementary Material online). This may very well Several English breeds such as the Suffolk, Texel, Dorset, and have facilitated the worldwide range of domestic sheep the Longwools became transboundary mutton sheep and adapting to a wide range of environments (Cao et al. 2021), have been used for improving local breeds on different con- for instance to the extreme conditions of the Arctic climate tinents (Kijas et al. 2012; Deniskova et al. 2018; Cao et al. and the Himalaya Highland. 2021), thus spreading the ancestry of European mouflon. Previous studies have showed genetic introgression of wild Alternatively, mouflon influence may be the result of a relatives in domestic animals such as sheep, goat, pig, dog, later introgression in a sympatric domestic population cattle, yak, cattle, rabbit, horse, and chicken (Ai et al. 2015; (Barbato et al. 2017; Cao et al. 2021). Historical records Medugorac et al. 2017; Jones et al. 2018; Wu et al. 2018; Lawal from ancient Rome and the 18th century have described et al. 2020; Wang et al. 2020; Zheng et al. 2020; Cao et al. 2021; hybridization between wild and domestic sheep in Europe Yu et al. 2021; ). The proportion of wild introgression varies (Cetti 2000). Even nowadays, farmers in Sardinia typically al- from 0.05% to 27.3% depending on the populations of do- low sheep to graze in the wild, where European mouflon mestic animals and the relevant wild relatives (supplementary reside (Lorenzini et al. 2011). Our estimates of the proportion table S48, Supplementary Material online). These introgressed of the mouflon component in domestic breeds (7.64–13.40%) variants and segments were located in genes with functions are higher than the 1.0–4.1% estimated by Barbato et al. such as environmental adaptation (e.g., the response to hyp- (2017) and the 0.6–1.3% estimated by Cao et al. (2021), oxia), innate immunity, olfactory relevance, and morpholog- both of which were on the basis of 50K genotypes. This ical traits (e.g., coat color, the polled, and the different shapes 15 Lv et al. . doi:10.1093/molbev/msab353 MBE of horns). Introgressed variants in genes associated with sim- originated in the Near East (around 5,000 years BP; Ryder ilar functions such as the oxygen transportation and hemo- 1984), as supported by the Y-chromosomal genetic variation globin complex have been observed in this and early studies analyses (Deng et al. 2020). Alternatively, an early mtDNA (Hu et al. 2019). Nevertheless, genes with different functions study suggested that sheep could have been introduced to were also detected in these studies (supplementary table S48, Indian subcontinent from Near East, probably by ancient sea Supplementary Material online), which could be due to dif- trade route (Singh et al. 2013). Overall, our study comple- ferentiated natural and artificial selection forces in a diverse mented the origins and migrations of sheep in previously range of environments and varied production systems. understudied regions and presented a comprehensive picture of expansions of sheep across the world (fig. 3E). We noted Demographic History of South-and-Southeast Asian that the estimated time of the main demographic events Sheep were in the approximate periods of the recent paleogenomic Reconstruction of the history of the Asian sheep populations studies (Lazaridis et al. 2016; Frachetti et al. 2017; McColl et al. by demographic modeling suggested that Central-and- 2018; Narasimhan et al. 2019; Shinde et al. 2019; Skourtanioti Eastern and South-and-Southeast Asian sheep diverged et al. 2020), but the time difference observed here could be 5.85–6.01 ka with subsequent gene flow from the Middle due to a limited number of modern samples and the lack of East 3.93–4.01 ka. This was possibly mediated by gene flow ancient samples in Central, South-and-Southeast Asia. Also, along the Southern Silk Road route and the Tibet–Nepal and there were controversial estimates of time for the population Burma–India routes (Pe 1959; Shaha 1970; Lv et al. 2015). In demographic histories in these paleogenomic investigations. the Holocene and Neolithic time, the major transformation of Additional samples in the regions (e.g., Kazakhstan, farming and pastoralism in Eurasia were accompanied by Afghanistan, Pakistan, Vietnam, Thailand, and Malaysia) in- human migrations (Narasimhan et al. 2019) and cultural cluding particularly ancient DNA will be needed to elucidate transitions, such as the Late Neolithic to Early Bronze Age the genetic origins and prehistoric population interactions of material cultures in the Near East, the Kura-Araxes culture in sheep in South-and-Southeast Asia. Our results of sheep his- the Southern Caucasus, the Yamnaya culture across the tory in South-and-Southeast Asia can provide insights into Eurasian Steppe, and the Indus Valley Civilization the waves of human migration in the regions. (Narasimhan et al. 2019; Skourtanioti et al. 2020). Our findings For a more complete visualization of the major migrations are consistent with recent paleogenomic evidence for the of domestic sheep (fig. 3E), we combined the evidence of the dynamic population history on the Eurasian Steppe and the modeling for the Asian populations (fig. 3D) with 1) the to- genetic makeup of most present-day South-and-Southeast pology with separate positions for British and Nordic breeds Asian populations during the Early and Late Bronze ages, (fig. 3A and supplementary fig. S11, Supplementary Material which share ancestry of: 1) Anatolian and Iranian farmers online) (Tapio et al. 2006; Kijas et al. 2012; Ciani et al. 2020); 2) across the broad region of Iranian plateau 6,500–6,000 BP the clear separation of North-and-West African sheep and (Narasimhan et al. 2019; Raghavan et al. 2019; Shinde et al. East-and-South African sheep (fig. 3A and supplementary fig. 2019); and 2) North Eurasian and Central Asian Steppe pas- S11, Supplementary Material online); and 3) tentatively, the toralists along the Inner Asian Mountain Corridor, the ancient model-based clustering for the American populations (fig. 3C) “Silk Roads” across Asia (e.g., Kazakhstan, Turkmenistan, (Kijas et al. 2012), which is in agreement with earlier results Afghanistan, and Pakistan) 4,500–4,000 BP (Majumder based on mtDNA and SNP BeadChip (Spangler et al. 2017; 2010; Frachetti et al. 2017; Lipson et al. 2017; Narasimhan Paim et al. 2021). et al. 2019). Thus, our results illustrate the impact of the Middle Eastern farming to both the sheep husbandry in Molecular Basis of Adaptive and Agronomic Traits Iran and to the pastoralism on the Eurasian Steppe spreading We performed several alternative methods for genome-wide eastward into South-and-Southeast Asia (Lazaridis et al. searching of selection signatures in order to identify novel 2016). trait-associated genes and variants. With our data set, PBS, The demographic reconstruction and population genetic FST, and CLR gave clear signals above a uniform background, differentiation supported close genetic relationship between which is in contrast to rather noisy fd and ROD scans. In South-and-Southeastern Asian and Central-and-Eastern general, the various approaches each detect different panels Asian populations (fig. 3A–D). Our results supported a rela- of genes with only little overlap, indicating for all methods a tively recent divergence and a close genetic affinity between high frequency of false negatives. At the same time, all meth- South-and-Southeastern Asian and Chinese sheep popula- ods are prone to delivering false positives, further complicat- tions (including Tibetan and Southwestern Chinese popula- ing the interpretation (Utsunomiya et al. 2015; Saravanan tions). The observation could also be explained by the strong et al. 2020). On the other hand, our large WGS data set allows gene flow from Central-and-Eastern Asian populations to to identify candidate causative mutations, at least in coding South-and-Southeastern Asian populations indicated by the sequences but also in plausible miRNA binding sites, and to modeling analyses (fig. 3D), which could be through several test their correlation with the breed-dependent phenotypes. popular ancient land (the Southern Silk Road route and the In addition, alignment of haplotype patterns (figs. 5D and 7D) Tibet–Nepal and Burma–India routes) (Pe 1959; Shaha 1970; for selected regions appeared to be informative. Lv et al. 2015). The later Near Eastern infusion could be as- From the multitude of potentially informative selection sociated with the introduction of fat-tailed sheep, which signatures, we have highlighted HBB and neighboring genes 16 Genetic Diversity, Introgression, and Agronomically Important Loci of Worldwide Sheep . doi:10.1093/molbev/msab353 MBE (adaptation to hypoxia), EVC (development), RXFP2 (horn genes (e.g., EFNA5, PIK3R3, and VTN) were involved in the development), TMEM154 (disease resistance), GDF6 (growth), PI3K-Akt signaling pathway, which were functionally related PDGFD (tail configuration), and MC1R (coat color), all of to cell proliferation and survival (supplementary fig. S28, which are key genes for the most relevant phenotypes of Supplementary Material online). Functions of these genes domestic sheep. For fleece fiber variation, we found an asso- and network of these relevant pathways indicated that ciation of mutated allele chr25:7,068,586(C) in the 30 -UTR of IRF2BP2, VEGFA, and oar-miR-20a-3p are the potentially im- IRF2BP2 with the medium- or fine-wool phenotype. Although portant factors regulating the fleece fiber diameter in domes- we observed homozygous genotypes of ancestral and derived tic sheep. alleles at the SNP mutation (chr25: T7,068,586C) in the hair We conclude that the historical wild-relative genetic intro- and coarse-wool sheep populations, the ancestral allele gression played an important role in shaping the genetic di- chr25:7,068,586(T) showed much higher frequencies than versity and agronomic traits and introduced adaptive variants those in medium-wool and fine-wool sheep (fig. 7C). Thus, in domestic sheep. The large data set presented here is a we concluded the ancestral allele chr25:7,068,586 (T) is asso- useful resource by providing a genomic framework for explo- ciated with the phenotypes of hair and coarse wool in wild ration of the genomic basis underlying phenotypic traits and and domestic sheep. genetic improvement of sheep. Variation between hairy and woolly fleece has been previ- ously found to be associated with an insertion of an antisense Materials and Methods EIF2S2 retrogene in the 30 -UTR of IRF2BP2, using a mainly French panel of populations (Demars et al. 2017). However, Data Sets, Samples, and DNA Extraction this association could not be entirely reproduced in our The sequencing data sets consisted of WGS generated in this worldwide panel of wild and domestic sheep, particularly in study and published previously, totaling 810 individuals from the wild species (with the phenotype of hair), hair sheep, and 7 wild relatives (72 individuals) and 158 domestic sheep pop- coarse-wool sheep populations (supplementary fig. S27, ulations (738 individuals) around the world. The novel Supplementary Material online). We estimated LD between sequences included genomes of one Asiatic mouflon and the mutation (chr25: T7,068,586C) and the insertion 165 individuals from 82 domestic populations (fig. 1A and (asEIF2S2) by calculating the parameter of LD (r2) (Hill and supplementary table S1, Supplementary Material online), Robertson1968). We found strong linkage between them in which are from regions that so far have been understudied wild species (frequency of IRF2BP2wt ¼ 1, frequency of in previous investigations, such as India, Indonesia, Malaysia, T ¼ 0.96) and fine-wool (r2 ¼ 0.27) sheep populations, and Pakistan (supplementary table S1, Supplementary whereas lower (r2 < 0.1) or no linkage was observed in the Material online). hair, coarse-wool, and medium-wool sheep populations (sup- The publicly available data sets included WGS of 644 indi- plementary table S49, Supplementary Material online). We viduals (572 domestic sheep, 32 Asiatic mouflon, 6 bighorn noted that only seven fine-wool individuals were included in sheep, 6 thinhorn sheep, 9 urial, 8 argali, 8 snow sheep, and 3 the LD estimation, which might have inadequate statistical European mouflon), and they were derived from five sources: power and significance. The results showed cosegregation the NextGen Consortium, the International Sheep Genomics between the ancestral allele chr25:7,068,586T and the absence Consortium (ISGC), and three of our previous studies (Deng of EIF2S2 insertion in the 30 -UTR of IRF2BP2 in wild sheep et al. 2020; Li et al. 2020; Chen et al. 2021). Summary infor- species. Thus, the expression of IRF2BP2 was regulated by mation of these newly sequenced and publicly available sam- both the EIF2S2 insertion and the SNP (chr25: T7,068,586C) ples, including population names, codes, geographic origins, in the 30 -UTR, which is associated with the fleece fiber phe- sample sizes, and contributors, is detailed in supplementary notypes in wild and domestic sheep. Our results provided table S1, Supplementary Material online. Coordinates of geo- complementary evidence for the genetic mechanisms under- graphic origins of the populations were provided by the con- lying the fleece fiber variation in sheep. Different from Demars tributors or were assigned as the centroid of their known core et al. (2017), we did not find IRF2BP2asEIF2S2/IRF2BP2asEIF2S2 distribution areas. Blood samples, ear or skin tissues were homozygotes. The difference in genotype pattern between collected; ear and skin samples were preserved in 95% ethanol this and the previous study could be due to that 1) multiple at 80 C. approaches and strict thresholds have been adopted in obtaining the accurate genotypes; 2) the EIF2S2 insertion Resequencing and Variant Calling was genotyped in much less samples of medium- and fine- DNA of the samples sequenced in this study was extracted wool sheep here, which, however, showed a high frequency in using the DNeasy Blood & Tissue kit (QiaGen, Shanghai, medium- and fine-wool sheep (Demars et al. 2017). China), including an RNase A treatment. DNA integrity was Besides these, we detected a set of candidate genes asso- inspected on agarose gels and concentration was quantified ciated with fleece fiber variation in sheep. Among these and using a Qubit 2.0 Fluorometer (Life Sciences, CA). At least 1mg previously identified candidate genes, 7 genes (e.g., VEGFA, of genomic DNA was used to construct a sequencing library HRAS, and AKT) were located in the pathways of VEGF sig- according to the manufacturer’s specifications for the TruSeq naling, which plays a central role in regulating cell prolifera- Nano Sample Prep Kit (Illumina Inc., San Diego, CA). Briefly, tion, migration, and survival; 17 genes (e.g., CACNA2D3, FGF7, the DNA was fragmented, end polished, A-tailed, and ligated and MAP4K4) were found in MAPK signaling pathway, and 16 with the full-length adapter. Fragments of 400–500 bp were 17 Lv et al. . doi:10.1093/molbev/msab353 MBE selected, PCR amplified and purified using the AMPure XP Sanger sequencing. The primers used for PCR were designed system (Beckman Coulter, IN). The prepared libraries were with FastPCR v6.7 (Kalendar et al. 2017; supplementary table assessed on an Agilent2100 Bioanalyzer and quantified using S5, Supplementary Material online). The PCR reactions were real-time PCR. Paired-end libraries that had average insert carried out in 20 ll volume containing 10 ll of 2taq PCR sizes of approximately 350 bp were sequenced on the Master Mix (GeneBetter Biotech, Beijing, China), 0.8 ll Illumina HiSeq X Ten platform (Illumina.) by Berry (10 pmol/ml) for each forward and reverse primer (supple- Genomics Co. Ltd, (Beijing, China), yielding 150-bp reads mentary table S5, Supplementary Material online), 1ll DNA with a target depth of 20-fold coverage per genome. templates (30–100 ng/ml), and the remainder supplied with We obtained 39 Gb of raw sequences for per sample and ddH2O. The reactions were performed by an Eppendorf 6331 used the FastQC v0.11.9 (https://www.bioinformatics.babraham. Flexlid Mastercycler Nexus Cycler with conditions of an initial ac.uk/projects/fastqc/, last accessed April 22, 2019) for assessing denaturation at 95 C for 15 min, followed by 35 cycles at a per-base sequence quality. Further processing and removal of 95 C for 30 s, annealing at 60 C for 90 s and extension at low-quality bases and artifact sequences were implemented 72 C for 30 s, and then a final extension at 72 C for 7 min. using the Trimommatic v.0.36 (Bolger et al. 2014). The high- The PCR products were sequenced by the Sanger method. All quality 150-bp/100-bp paired-end reads were aligned to the the reads were assessed manually and genotypes of each site sheep reference genome Oar_rambouillet_v1.0. (https://www. were called by the sequencing peaks. Subsequently, we com- ncbi.nlm.nih.gov/assembly/GCF_002742125.1/, last accessed pared genotypes of each site identified by whole-genome May 12, 2019), using the Burrows–Wheeler aligner (BWA resequencing and obtained by the Sanger sequencing for mem) v.0.7.8 (Li and Durbin 2009) with default parameters. the same individuals. We then converted the mapping reads into bam files and sorted the files using the SAMtools. Duplicates were removed by the Structural Variant Calling MarkDuplicates module in GATK v 4.1.2.0 (McKenna et al. SVs were called from the sequences of 534 high-depth 2010). We calculated whole-genome sequencing coverage and samples (>15). We used three approaches of the depth of each sample using the SAMtools v.1.9. We only kept Manta v1.3.2 (Chen et al. 2016), Delly v0.8.3 (Rausch et al. samples with depth >10 to call short variations, and selected 2012), and smoove-0.2.6 (https://github.com/brentp/ samples with depth >15 to identify SVs. smoove, last accessed June 15, 2020) with default parame- SNP and indels were called from the bam files by the GATK ters to detect SVs such as deletions (DEL), inversion (INV), HaplotypeCaller module with the GATK best-practice recom- and duplications (DUP) with sizes of 50 bp–1Mb and trans- mendations (McKenna et al. 2010). Raw GVCFs with the sam- locations (TRA). We filtered all SVs based on the following ples called individually were merged using the CombineGVCFs criteria: 1) 50 bp  SV1 Mb; and 2) the quality filtering and called for SNPs using the GenotypeGVCFs. We then selected parameters suggested by two approaches (flag PASS). The the candidate SNPs and created the selected SNP data using the smoove-0.2.6 was implemented to combine the three SV GATK module SelectVariants. To avoid potential false-positive data sets generated by the above three software and per- calls, we implemented “VariantFiltering” of the GATK for the form genotyping. The merged SV genotypes were retained selected SNPs and INDELs using the best practice parameters following the criteria: 1) genotypes of SVs detected by all “QUAL < 30.0 | QD < 2.0 | MQ < 40.0 | FS > 60.0 | SOR > 3.0 the three approaches were identical for at least two | MQRankSum < 12.5 | ReadPosRankSum < 8.0” and “QD approaches; and 2) for SVs detected only by two out of < 2.0 | QUAL < 30.0 | FS > 200.0 | ReadPosRankSum < the three approaches, their genotypes identified by 20.0,” respectively. We then filtered out nonbiallelic SNPs. smoove-0.2.6, or Manta v1.3.2 were kept. After the quality screening, all the identified SNPs were further annotated using SnpEff v4.3t (Cingolani et al. 2012) Genetic Diversity based on the gene annotations of the sheep reference ge- We used the filtered set of 121.2 million high-quality SNPs for nome Oar_rambouillet_v1.0. Locations for SNPs in various estimation of genetic diversity, population differentiation and genic and intergenic regions as well as synonymous or non- LD decay metrics including expected heterozygosity (He), nu- synonymous SNPs in exonic regions were annotated. cleotide diversity (p), runs of homozygous (ROH), and LD Additionally, we enriched for SNPs specific for species, popu- decay. We summarized the numbers of homozygous and lation, or geographic groups of populations using the SnpEff heterozygous SNPs against the reference genome. Values of v4.3t. We annotated each of the significant species-, popula- p and He were estimated for nonoverlapping 1 Mb regions tion- or region-specific enriched SNPs, in particular the non- along the genome of each individual using the vcftools v0.1.17 synonymous SNPs, using SnpEff v4.3t. (Danecek et al. 2011). We used PLINK v1.90p (Purcell et al. 2007) (–homozyg-density 50 –homozyg-window-het 1 – SNP Validation homozyg-window-kb 200 –homozyg-window-snp 50) to es- To check the confidence of SNPs called, we compared the timate the number (NROH) and size (SROH) of ROH and indi- SNPs identified with the O. aries dbSNP v.151 (http://www. vidual autozygosity (FROH) (McQuillan et al. 2008). ncbi.milm.nih.gov/SNP, last accessed July 10, 2021). In addi- Overall genetic differentiation across domestic popula- tion, we validated the SNPs in specific genes by the Sanger tions measured by Weir and Cockerham’s estimator of FST method of DNA sequencing. A total of 133 SNPs from 16 (Weir and Cockerham 1984) was calculated using the vcftools individuals of ten populations were genotyped by PCR and v0.1.17 with 100 kb bins. FST values between populations or 18 Genetic Diversity, Introgression, and Agronomically Important Loci of Worldwide Sheep . doi:10.1093/molbev/msab353 MBE groups of populations were estimated with 220,350 unlinked and-East Asian, the Middle Eastern, South-and-Southeast SNPs using the gcta_1.93.2beta (Yang et al. 2011). The statis- Asian, and American populations). Bi-allelic SNPs of the cor- tical significance of FST estimates was tested using the ap- responding individual were extracted from the raw vcf files. proach t test implemented in the R program (R Core Team We further masked low-complexity regions of the reference 2020). genome Oar_rambouillet_v1.0 following the method of Li We examined the patterns of LD decay within each species Heng (http://lh3lh3.users.sourceforge.net/snpable.shtml, last or population by random selection of three individuals in accessed October 20, 2020). The mutation rate was set at order to avoid biases by difference in sample sizes and using 1.51 10-8 mutations per nucleotide site per generation the vcftools v0.1.17 to extract individual data for LD analysis. (Chen et al. 2019). All other parameters were set to their Pairwise LD estimates were measured as parameter r2 with a default values while the same mutation rate (Chen et al. maximum distance of 300 kb using the PopLDdecay v3.41 2019) and generation time (Zhao et al. 2017; Alberto et al. (Zhang et al. 2019). 2018) were applied to the final representation of recent de- mographic events of Ne. Population Genetics Analysis Historical effective population size (Ne) was also estimated To remove the potential impact of close relatives on estima- using the program SNeP v.1.11 (Barbato et al. 2015) with tion, we calculated the pair-wise KING-robust kinship estima- default settings. The program SNeP estimates Ne based on tor using the KING v2.2.5 (Manichaikul et al. 2010). We the relationship between r2, Ne, and c (recombination rate) excluded one of each sample pair showing close relatedness following the equation E (r2) ¼ 1/(1 þ 4 Nec) (Sved 1971). with kinship coefficient > 0.0884 (e.g., duplicates/monozy- Only the SNPs with nonmissing genotypes and MAF > gotic twins, 1st-degree and 2nd-degree relatives) (Manichaikul 0.05 were included in this estimation. We used the different et al. 2010), which resulted in 714 unrelated individuals for SNP marker distance bins in the LD (r2) analysis to estimate genetic diversity analyses (supplementary table S13, Ne at t ¼ 1/2c generation ago. Supplementary Material online). We further filtered the SNP data set with a MAF < 0.05, Hardy–Weinberg equilib- Migration and Divergence Modeling rium <0.001 and a proportion of missing genotypes >10%. Data Set We implemented LD pruning with the PLINK option “– Demographic history of sheep postdomestication was recon- indep-pairwise 50 5 0.2.” After the filtering and LD pruning, structed using the fastsimcoal v2.6 (Excoffier et al. 2013). 24,363,608 SNPs were retained for population genetics Forty-eight high-depth (15.1–22. 6) domestic samples analyses. from Africa (n ¼ 10), Europe (n ¼ 10), Central-and-East Asia First, we calculated the identity-by-state (IBS) genetic dis- (n ¼ 10), South-and-Southeast Asia (n ¼ 10), and the Middle tance matrix between the individuals using the PLINK v1.90 East (n ¼ 8) were selected as the representative samples of (–distance 1-ibs) and visualized by an unrooted neighbor- different geographic and genetic origins based on their pop- joining (NJ) phylogenetic tree using the SplitsTree5 v5.1.7- ulation history and results of the population genetics analyses beta (Huson and Bryant 2006). PCA was performed using (supplementary fig. S29 and table S50, Supplementary the smartpca function built in the EIGENSOFT v.6.0.1 Material online). (Patterson et al. 2006) with default settings. Furthermore, For these statistical modeling analyses, we first selected model-based clustering was performed using the sNMF v1.2 6,477,008 high-quality autosomal SNPs that: 1) are located with the number of cluster K from 1 to 20. The optimal K was within intergenic regions and 2) are outside CpG islands, de- determined by the minimal value of the cross-entropy crite- fined according to the UCSC annotation (https://hgdown- rion of sNMF. load.soe.ucsc.edu/hubs/GCF/002/742/125/GCF_002742125.1/ bbi/, last accessed October 27, 2020), 3) do not have missing Demographic Inference genotypes, 4) have a depth >10 in all samples, and 5) have We used the PSMC model (Li and Durbin 2011) to recon- an ancestral state inferred from the goat and cattle reference struct changes in historical effective population size (Ne) genomes. We then generated the multidimensional SFS file through time. Analyses within species or populations were using the python script vcf2sfs.py (https://github.com/mar- based on the three individuals with the highest mean se- queda/SFS-scripts, last accessed November 12, 2020). In order quencing coverage per population. The psmc function built to calculate confidence intervals for the parameters, we in the PSMC package was used to estimate Ne with the resampled 2,647  1-Mb blocks on the autosomes with sim- parameters “–N25 –t15 –r5 –p ‘4 þ 25*2 þ 4 þ 6.” We ilar cumulative genomic length as in the observed data used a mutation rate of 1.51  10-8 mutations per nucleotide (2,234 Mb) and generated 100 multidimensional SFS data site per generation (Chen et al. 2019) and a generation time of sets. 3 years (Zhao et al. 2017; Alberto et al. 2018) to convert the resulting scaled values into years and individuals. Statistical Modeling We also inferred the dynamics of recent historical Ne dur- We investigated the migration and divergence of the South- ing the past 1,000 to 100,000 years using the SMCþþ and-Southeast Asian populations and estimated the time of v1.14.0.dev0 (Terhorst et al. 2017). We selected 7–10 high- their main expansions and gene flows among populations in depth unrelated individuals in each of the six geographically different continents. On the basis of the model inferred by Lv differentiated genetic groups (African, European, Central- et al. (2015) and Zhao et al. (2017), three alternative models 19 Lv et al. . doi:10.1093/molbev/msab353 MBE without gene flow among populations from the different AdmixTools v7.0.1 (Patterson et al. 2012). We used Menz main geographic regions were tested for their origins and sheep, a domestic population which did not show any signal demographic history of South-and-Southeast Asian popula- of wild introgression in the admixture analysis, as the refer- tions (supplementary fig. S13, Supplementary Material on- ence population (W), individual populations of domestic line): scenario I: descent of the ancestral Middle Eastern sheep as the target population (X), and each wild species as breeds; scenario II: descent of Central-and-East Asian breeds; the donor (Y). We assumed that outgroup (Z) was fixed for and scenario III: originated from early admixture between the ancestral alleles when the target population was snow Central-and-East Asian and the Middle Eastern breeds. We sheep, bighorn, or thinhorn sheep, whereas bighorn sheep assumed a ghost population as the common ancestor of was used as the outgroup (Z) when Asiatic mouflon, modern domestic sheep that was present in the Middle European mouflon, urial, or argali was treated as the donor East after domestication, but before their expansions to other population. A value of jZj > 3 was considered to support an regions. influence of the donor wild species and tested domestic In the optimal model, asymmetric gene flow between populations. modern populations from neighboring geographical regions were added, such as between African and European breeds, Localization and Quantification of Introgression between African and the Middle Eastern breeds, between the In order to quantify the introgressed proportions of wild Middle Eastern and South-and-Southeast Asian breeds, be- relatives in the genomes of domestic sheep and map the tween the Middle Eastern and Central-and-East Asian pop- introgressed loci at the chromosome-wide level, we further ulations, between South-and-Southeast Asian and Central- conducted fd analysis (Martin et al. 2015) based on the results and-East Asian populations, but not between the ghost and of D statistics. Statistics fd, a modified version of f estimator modern populations (supplementary fig. S30, Supplementary (Green et al. 2010), is prone to extreme values in the genomic Material online). regions of low diversity, providing an accurate and robust method for detecting and quantifying the introgressed loci Maximum-Likelihood Parameter Estimation (Martin et al. 2015). The domestic populations with jZ scoresj We estimated the maximum likelihood values for scenarios I, > 3 were selected as potential recipient populations. For each II, and III, independently, after performing 1,000,000 simula- individual in the recipient populations, we used fd to quantify tions, 65 conditional maximization cycles, and 100 optimiza- the proportion of introgression as well as to identify the tion runs starting from random initial conditions. We used introgressed regions. For the windows with fd < 0, the fd approaches described in Choin et al. (2021) and Malaspinas statistic doesn’t have biological meaning. Therefore, the fd et al. (2016) to optimize the fit between the expected and values were converted to 0. We then calculated P value of observed SFS (supplementary fig. S31, Supplementary each fd by Z transform, and further corrected the P values Material online). We calculated the maximum likelihood of with the Benjamini-Hochberg’s FDR method (Benjamini and each model with all entries of SFS for the first 25 cycles (-l 25) Hochberg 1995). and then optimized the maximum likelihood of each model Considering the distribution of introgressed tract for the remaining 40 cycles (-L 65). lengths and to minimize nonindependences among adja- cent windows and ensure sufficient resolutions, we used Nonparametric Bootstrapping Analysis ABBABABAwindows.py to calculate the fd statistic using For the optimal model, we estimated confidence intervals for sliding windows of 100 kb with a step of 50 kb across the the parameters with nonparametric bootstrapping approach genomes (Martin et al. 2016). We set the minimum SNPs in using the 100 multidimensional SFS data sets. We re- each window of 500 with the option “-w 100000 -m 500 -s estimated parameters using the fastsimcoal v2.6 with same 20000.” settings as for the original observed data set. We performed We used the popgenWindows.py to calculate p, FST, and 20 replicate runs per bootstrap data set starting from random dxy as a measure of genetic variation, genetic differentiation, initial conditions. The 95% CIs for each parameter were and absolute divergence, with the same sliding windows obtained using the parameter estimates of all bootstrapping parameters used in the fd analysis, for the following two replicates with the R boot package. All analyses were imple- groups of comparisons: donor versus recipient population mented using the fastsimcoal v2.6 with the same mutation (-p donor -p recipient) and donor versus reference popula- rate (Chen et al. 2019) and generation time (Zhao et al. 2017; tion (-p donor -p reference) (Martin et al. 2015). We con- Alberto et al. 2018). verted the meaningless or noisy fd value (D < 0, or D > 0, but fd > 1) to zero because theoretically only a positive fd statistic Interspecies Introgression Analysis value was considered as evidence of introgression from pop- D Statistics ulation P3 into P2 under a four-taxon topology ((P1, P2), P3, We tested for the presence of interspecific genetic introgres- O) (Martin et al. 2015). We estimated the fd statistic value sion from wild relatives into domestic populations using using Menz sheep as P1, individual domestic populations as different statistical analyses such as D statistics (Patterson P2, individual wild species as P3, and ancestral alleles (when P3 et al. 2012) and fd statistics (Martin et al. 2015). We calculated was snow sheep, bighorn sheep, or thinhorn sheep) or big- the D statistics for the four-taxon model (W, X, Y, Z) horn sheep (when P3 was Asiatic mouflon, European mou- using the function qpDstat implemented in the flon, urial or argali) as O in the analysis. 20 Genetic Diversity, Introgression, and Agronomically Important Loci of Worldwide Sheep . doi:10.1093/molbev/msab353 MBE Incomplete Lineage Sorting selection on phenotypes such as dwarf (e.g., Ouessant sheep), We inferred that a genomic tract in the gene HBB was intro- we calculated CLR scores using the SweeD (Pavlidis et al. gressed from one of the wild sheep species into Changthangi 2013). We used 10-kb nonoverlapping sliding windows in sheep. To verify that the target genomic tract is the result of the tested Ouessant population. Empirical P values were cal- introgression, we calculated the probability of the alternative culated for CLR windows and the top 5% of windows were scenario of ILS. The expected length of a shared ancestral considered as candidate signatures of selection. sequence is L ¼ 1/(rt), where r is the recombination rate per generation per base pair and t is the branch length be- Selection Related to High-Altitude Adaptation tween argali/Asiatic mouflon/European mouflon/urial and We conducted PBS analysis (Yi et al. 2010) to identify the domestic sheep since their divergence. The probability of a selection signatures associated with high altitude adaptation. length of > m is 1  GammaCDF (m, shape ¼ 2, r ¼ 1/L), Because we detected the introgression signals in HBB and where GammaCDF is the Gamma distribution and m is the HBBC genes from argali to Changthangi sheep in Ladakh length of introgressed tracts (Huerta-Sanchez et al. 2014). (ca. 4,000 m above sea level), we used Changthangi sheep (CHA) as the representative high-altitude population, East Detection of Selective Signatures Friesian sheep (EFR, living around the sea level) as the low- Selective Signatures in Domestic Sheep during Improvement altitude population, and Hu sheep (HUS, living around the sea To identify the genomic signatures of selection in domestic level) as the reference low-altitude population to perform the sheep, we calculated the average SNP FST values (Weir and PBS analysis. We used R scripts from the github (https:// Cockerham 1984) using the VCFtools v1.17. All genomic github.com/genevol-usp/EvoGen_course, last accessed regions in 50-kb sliding windows were scanned across genomes January 5, 2021) to estimate pairwise population FST values of the 158 populations using 10-kb steps. Regions with top 1% between three comparisons of Changthangi versus East of the average FST distribution were identified as selection Friesian sheep, Changthangi versus Hu sheep and East signatures. Furthermore, we examined the changes in genetic Friesian versus Hu sheep. PBS value was calculated following diversity (as measured in p) due to artificial selection. The ROD the equation: PBS ¼ ((log(1  FST.EFR.CHA)) þ (log(1  values in 10-kb nonoverlapping windows were calculated be- FST.CHA.HUS))  (log(1  FST.EFR.HUS)))/2 (Yi et al. 2010). tween landraces and improved populations using the formula: The PBS values with 99th empirical thresholds were identified ROD ¼ 1  pimproved populations/plandraces (Wu et al. 2020). as outlier signals under selection. Gene Annotation, GO, and KEGG Analyses Selection Associated with Fleece Fiber Diameter We annotated genes associated with the candidate SNPs, We tested for the potential selective sweeps in the pairwise CNVs, or genomic regions using the sheep reference assembly comparisons between groups of hair, coarse-wool, medium- Oar_rambouillet_v1.0. GO term enrichment and KEGG path- wool, and fine-wool populations (supplementary table S32, way analyses of a candidate gene set were carried out by a Supplementary Material online). A XP-CLR test was imple- statistical overrepresentation test with the default setting in mented to scan the genomes for selective sweeps using the the Database for Annotation, Visualization, and Integrated program XP-CLR v1.0 (Chen et al. 2010). The SNPs with less Discovery (DAVID) database v6.8 (Huang et al. 2009). than 10% missing data were allowed (“–max-missing 0.9”) in Categories with at least six genes and the threshold of ad- the analysis. We calculated XP-CLR scores using the grid points justed P value < 0.05 after the Bonferroni correction were spaced by 2 kb, a maximum of 200 SNPs in a window of 0.5 considered as significantly enriched terms and pathways. cM, and the down-weighting contributions of highly corre- lated SNPs (r2 > 0.95) with the parameters (-w1 0.005 200 Histological Analysis 2000 2 -p0 0.95). We considered genomic regions with the top To examine the numbers and the ratio of secondary to pri- 1% region-wise XP-CLR scores as candidate selective sweeps. mary wool follicles in the important development stage of Moreover, we conducted a selection test based on the SVs by embryonic skin (embryo 85–90 days), fixed dorsal skin sam- calculating the statistic FST (Bertolotti et al. 2020). For each SV, ples for coarse-wool (e.g., Wadi sheep) and fine-wool (e.g., we calculated the FST values between the above-mentioned Chinese Merino) sheep were dehydrated with gradient alco- pairwise comparisons of groups of populations with different hol, processed in paraffin and cut into histological longitudi- fleece fiber diameters using the VCFtools v0.1.17. To identify nal sections of 5lm. The dorsal skin sections were processed the significance of each FST value of individual SVs, we gener- by dewaxing, stained with H&E, and photographed for the ated 200 randomly sampled FST values, and obtained P values identification of the secondary and primary wool follicles. of per SV by calculating the proportion of FST values obtained in the 200 random distributions higher than the FST estimates RNA-seq and miRNA-seq Analysis in the observed distribution (Bertolotti et al. 2020). We se- RNA-seq Analysis lected the SVs with the P < 0.01 as the candidate selective SVs. RNA-seq data for samples from eight fine-wool and medium- wool populations (e.g., Aohan fine wool, Chinese Merino, Selection for Particular Phenotypes Churra, Erdos fine -wool, Gansu Alpine fine wool, Rambouillet, To identify alleles that are either close to fixation or already Super fine wool Merino, and Texel sheep) and seven coarse-wool reached fixation in particular populations due to past populations (e.g., Bashibai, Hu, Kazakh, Min Xian Black Fur, 21 Lv et al. . doi:10.1093/molbev/msab353 MBE Small-tailed Han, Tan, and Tibetan sheep) were retrieved from 2000. WT/mutated IRF2BP2 30 -UTR and miRNAs mimics public database (supplementary table S45, Supplementary were generated according to the mRNA and miRNA sequen- Material online). RNA-seq reads were quality filtered with the ces (supplementary table S46, Supplementary Material FastQC v0.11.9 (https://www.bioinformatics.babraham.ac.uk/ online). projects/fastqc/, last accessed April 22, 2019). Low-quality ends HEK293T cells were transfected with Lipofectamin 2000 in of the reads were removed using Cutadapt (Martin 2011). Reads 24-well plates. The psiCHECK-2 vector containing WT or mu- from each sample were aligned to the Oar rambouillet v1.0 tated (MU) 30 -UTR of IRF2BP2 vector was cotransfected into reference genome (https://www.ncbi.nlm.nih.gov/assembly/ HEK293T with mimic or negative control of miRNA (i.e., WT GCF_002742125.1, last accessed May 12, 2019) using the þ miRNA mimic, WT þ negative control, MU þ miRNA TopHat (Trapnell et al. 2012) with default options. Transcript mimic, and MU þ negative control). One day after the trans- quantification, normalization, and assembly were implemented fection, the cells were lysed, and then the Renilla and Firefly using the Cufflinks (Trapnell et al. 2012). The FPKM of IRF2BP2 substrates were added using the dual-luciferase reporter kit gene (chr25:7,067,974–7,071,785) was then calculated to quantify (Promega, WI) following the manufacturer’s recommenda- its expression level. tions. The Renlilla luciferase signals were normalized to firefly luminescence. These experiments were repeated three times. Prediction of the Target miRNA for IRF2BP2 A two-sided Student’s t-test was performed to evaluate sta- tistical significance of the differences between the miRNA Three different programs Targetscan (Agarwal et al. 2015), negative control and miRNA mimics. miRBase (Kozomara et al. 2019), and miRDeep2 (Mackowiak 2011) were used to jointly predict the target Supplementary Material miRNA based on the mRNA sequence of the IRF2BP2 gene (fine-wool sheep: 50 -. . .GUUGGUUACGUAAUACACA. . .-30 Supplementary data are available at Molecular Biology and and coarse-wool sheep: 50 -. . .GUUGGUUACAUAAUACACA Evolution online. 0 . . .-3 ). The program Targetscan was used to analyze the evo- Acknowledgments lutionary conservation of the miRNAs bound to the 30 -UTR of IRF2BP2, and the binding sites of IRF2BP2 and correspond- This study was financially supported by grants from the ing miRNAs in different species (e.g., sheep, cattle, and goat). National Natural Science Foundation of China (Nos. We used the miRBase database to obtain miRNA sequences 32061133010, 31825024, 31661143014, and 31972527), the of different species, evaluated site conservation, and selected National Key Research and Development Program-Key miRNAs whose binding sites were relatively conserved in Projects of International Innovation Cooperation between sheep and other species. The target identified by all the three Governments (2017YFE0117900), the External Cooperation programs above was considered to be the potential target Program of Chinese Academy of Sciences miRNA for IRF2BP2. (152111KYSB20150010), and the Second Tibetan Plateau Scientific Expedition and Research Program (STEP) (No. 2019QZKK0501), and the Special Funds of the State Key miRNA-seq Analysis Laboratory of Sheep Genetic Improvement and Healthy Small RNA raw data were obtained from public database and Production (Nos. 2018CA001, and 2019CA009, filtered through the Cutadapt (Martin 2011) by excluding 2020CA001). We express our thanks to the owners of the adaptor sequences, low-quality reads, and sequences sheep for donating samples (see supplementary table S1, <18 nt or >30 nt (supplementary table S47, Supplementary Supplementary Material online). For the Ethiopian popula- Material online). Reads were mapped to the Oar rambouillet tions, the sampling and DNA extraction were supported by v1.0 reference genome and aligned with mature miRNA the CGIAR Research Program on Livestock (CRP Livestock) sequences in the miRBase database to identify known and the Chinese government contribution to the CGIAR miRNAs. The expression levels of miRNAs in coarse-wool through ICARDA. We also thank Juan Deng, Ya-Xi Xu, Jia- and fine-wool sheep breeds were estimated by using the Nan Jing, Amadou Traore (Institut de l’Environnement et sRNABench (Barturen et al. 2014), which normalized reads de Recherches Agricoles [INERA], Ouagadougou, Burkina count number of each miRNA reads per million (RPM) using Faso), and Masroor Ellahi Babar (Virtual University, Lahore, the following formula: RPM ¼ (miRNA reads number/total Pakistan) for their help during sample collection. For mapped reads per library)  1,000,000. Iranian samples, data collection was supported by Iran National Science Foundation (INSF Grant No. 98028814). Cell Lines and Dual Luciferase Reporter Assay The Chinese government contribution to CAAS-ILRI Joint HEK293T cells were used to validate the IRF2BP2 target. Cells Laboratory on Livestock and Forage Genetic Resources in were cultured in the Dulbecco0 s modified eagle medium plus Beijing is appreciated. 10% fetal bovine serum and penicillin/streptomycin. The cells were propagated at 37 C, 90% of air humidity and with 5% of Author Contributions CO2. All reagents used for cell culture were provided by the M.-H. L. designed and supervised the study. F.-H.L. and Y.-H.C. Invitrogen/Gibco (Beijing, China). The cells were seeded into performed the genome data analyses. Y.-H. C. did the labo- 24-cell plates and transfected 24 h later with Lipofectamine ratory work. M.-H. L. wrote the manuscript with main 22 Genetic Diversity, Introgression, and Agronomically Important Loci of Worldwide Sheep . doi:10.1093/molbev/msab353 MBE contributions from F.-H. L., Y.-H.C., J.-L. H., and J.A.L. All Chen X, Schulz-Trieglaff O, Shaw R, Barnes B, Schlesinger F, Kallberg M, authors have prepared the samples or provided help during Cox AJ, Kruglyak S, Saunders CT. 2016. Manta: rapid detection of structural variants and indels for germline and cancer sequencing the sample collection. All authors reviewed and approved the applications. Bioinformatics 32(8):1220–1222. final version of the manuscript. Chen ZH, Xu YX, Xie XL, Wang DF, Aguilar-Gomez D, Liu GJ, Li X, Esmailizadeh A, Rezaei V, Kantanen J, et al. 2021. Whole-genome sequence analysis unveils different origins of European and Asiatic mouflon and domestication-related genes in sheep. Commun. Biol. Data Availability 4:1307. The whole-genome sequences reported in this study are avail- Chessa B, Pereira F, Arnaud F, Amorim A, Goyache F, Mainland I, Kao RR, able upon request for research purpose. Pemberton JM, Beraldi D, Stear MJ, et al. 2009. Revealing the history of sheep domestication using retrovirus integrations. Science 324(5926):532–536. References Choin J, Mendoza-Revilla J, Arauna LR, Cuadros-Espinoza S, Cassar O, Larena M, Ko AMS, Harmant C, Laurent R, Verdu P, et al. 2021. Agarwal V, Bell GW, Nam JW, Bartel DP. 2015. Predicting effective Genomic insights into population history and biological adaptation microRNA target sites in mammalian mRNAs. eLife 4:e05005. in Oceania. Nature 92:583–589. Ai H, Fang X, Yang B, Huang Z, Chen H, Mao L, Zhang F, Zhang L, Cui L, Ciani E, Mastrangelo S, Da Silva A, Marroni F, Ferencakovic M, Ajmone- He W, et al. 2015. Adaptation and possible ancient interspecies Marsan P, Baird H, Barbato M, Colli L, Delvento C, et al. 2020. On the introgression in pigs identified by whole-genome sequencing. Nat origin of European sheep as revealed by the diversity of the Balkan Genet. 47(3):217–225. breeds and by optimizing population-genetic analysis tools. Genet Alberto FJ, Boyer F, Orozco-terWengel P, Streeter I, Servin B, de Sel Evol. 52(1):25. Villemereuil P, Benjelloun B, Librado P, Biscarini F, Colli L, et al. Cingolani P, Platts A, Wang le L, Coon M, Nguyen T, Wang L, Land SJ, Lu 2018. Convergent genomic signatures of domestication in sheep X, Ruden DM. 2012. A program for annotating and predicting the and goats. Nat Commun. 9(1):813. effects of single nucleotide polymorphisms, SnpEff: SNPs in the ge- Aldersey JE, Sonstegard TS, Williams JL, Bottema CDK. 2020. nome of Drosophila melanogaster strain w1118; iso-2; iso-3. Fly Understanding the effects of the bovine POLLED variants. Anim (Austin). 6(2):80–92. Genet. 51(2):166–176. Crispim BA, Matos MC, Seno LO, Grisolia AB. 2012. Molecular markers Barbato M, Hailer F, Orozco-terWengel P, Kijas J, Mereu P, Cabras P, for genetic diversity and phylogeny research of Brazilian sheep Mazza R, Pirastru M, Bruford MW. 2017. Genomic signatures of breeds. Afr J Biotechnol. 11(90):15617–15625. adaptive introgression from European mouflon into domestic sheep. Danecek P, Auton A, Abecasis G, Albers CA, Banks E, DePristo MA, Sci Rep. 7(1):7623. Handsaker RE, Lunter G, Marth GT, Sherry ST, et al. 2011. The variant Barbato M, Orozco-terWengel P, Tapio M, Bruford MW. 2015. SNeP: a call format and VCFtools. Bioinformatics 27(15):2156–2158. tool to estimate trends in recent effective population size trajectories De Preter K, Barriot R, Speleman F, Vandesompele J, Moreau Y. 2008. using genome-wide SNP data. Front Genet. 6:109. Positional gene enrichment analysis of gene sets for high-resolution Benjamini Y, Hochberg Y. 1995. Controlling the false discovery rate: a identification of overrepresented chromosomal regions. Nucleic practical and powerful approach to multiple testing. J R Stat Soc Acids Res. 36(7):e43 Series B 57(1):289–300. Demars J, Cano M, Drouilhet L, Plisson-Petit F, Bardou P, Fabre S, Servin Bertolotti AC, Layer RM, Gundappa MK, Gallagher MD, Pehlivanoglu E, B, Sarry J, Woloszyn F, Mulsant P, et al. 2017. Genome-wide identi- Nome T, Robledo D, Kent MP, Røsæg LL, Holen MM, et al. 2020. The fication of the mutation underlying fleece variation and discriminat- structural variation landscape in 492 Atlantic salmon genomes. Nat ing ancestral hairy species from modern woolly sheep. Mol Biol Evol. Commun. 11(1):5176. 34(7):1722–1729. Bolger AM, Lohse M, Usadel B. 2014. Trimmomatic: a flexible trimmer for Demirci S, Koban Baştanlar E, Dagtaş ND, Pişkin E, Engin A, Ozer F, Illumina sequence data. Bioinformatics 30(15):2114–2120. Yüncü E, Dogan SA, Togan I. 2013. Mitochondrial DNA diversity Bouwman AC, Daetwyler HD, Chamberlain AJ, Ponce CH, Sargolzaei M, of modern, ancient and wild sheep (Ovis gmelinii anatolica) from Schenkel FS, Sahana G, Govignon-Gion A, Boitard S, Dolezal M, et al. Turkey: new insights on the evolutionary history of sheep. PLoS One. 2018. Meta-analysis of genome-wide association studies for cattle 8(12):e81952. stature identifies common genes that regulate body size in mam- Deng J, Xie X-L, Wang D-F, Zhao C, Lv F-H, Li X, Yang J, Yu J-L, Shen M, mals. Nat Genet. 50(3):362–367. Gao L, et al. 2020. Paternal origins and migratory episodes of domes- Bradford GE, Fitzhugh HA. 1984. Hair sheep of Western Africa and the tic sheep. Curr Biol. 30(20):4085–4095. Americas. In: A genetic resource for the tropics. Vol. 2. 1st ed. Boca Deniskova TE, Dotsev AV, Selionova MI, Kunz E, Medugorac I, Reyer H, Raton (FL): Cambridge University Press. p. 163–170. Wimmers K, Barbato M, Traspov AA, Brem G, et al. 2018. Population Campos E, Cuellar J, Salvador O, Garcıa-Trejo EA, Pereira F. 2020. The structure and genetic diversity of 25 Russian sheep breeds based on genetic diversity and phylogeography of Mexican domestic sheep. whole-genome genotyping. Genet Sel Evol. 50(1):1–16. Small Ruminant Res. 187:106109. Dong K, Yang M, Han J, Ma Q, Han J, Song Z, Luosang C, Gorkhali NA, Cao YH, Xu SS, Shen M, Chen ZH, Gao L, Lv FH, Xie XL, Wang XH, Yang Yang B, He X, et al. 2020. Genomic analysis of worldwide sheep H, Liu CB, et al. 2021. Historical introgression from wild relatives breeds reveals PDGFD as a major target of fat-tail selection in sheep. enhanced climatic adaptation and resistance to pneumonia in BMC Genomics. 21(1):1–12. sheep. Mol Biol Evol. 38(3):838–855. Dunmire WW. 2013. New Mexico’s Spanish Livestock Heritage: Four Cetti F. 2000. Storia Naturale di Sardegna. Nuoro, Italy: Ilisso. Centuries of Animals, Land, and People. Albuqueque: UNM Press. Chen H, Patterson N, Reich D. 2010. Population differentiation as a test Estienne A, Lahoz B, Jarrier P, Bodin L, Folch J, Alabart JL, Fabre S, Monniaux for selective sweeps. Genome Res. 20(3):393–402. D. 2017. BMP15 regulates the inhibin/activin system independently of Chen L, Qiu Q, Jiang Y, Wang K, Lin Z, Li Z, Bibi F, Yang Y, Wang J, Nie ovulation rate control in sheep. Reproduction 153(4):395–404. W, et al. 2019. Large-scale ruminant genome sequencing Excoffier L, Dupanloup I, Huerta-Sanchez E, Sousa VC, Foll M. 2013. provides insights into their evolution and distinct traits. Science Robust demographic inference from genomic and SNP data. PLoS 364(6446):6446. Genet. 9(10):e1003905. Chen N, Cai Y, Chen Q, Li R, Wang K, Huang Y, Hu S, Huang S, Zhang H, Frachetti MD, Smith CE, Traub CM, Williams T. 2017. Nomadic ecology Zheng Z, et al. 2018. Whole-genome resequencing reveals world- shaped the highland geography of Asia’s Silk Roads. Nature wide ancestry and adaptive introgression events of domesticated 543(7644):193–198. cattle in East Asia. Nat Commun. 9:1–13. 23 Lv et al. . doi:10.1093/molbev/msab353 MBE Frantz LAF, Bradley DG, Larson G, Orlando L. 2020. Animal domestica- Li M, Tian S, Jin L, Zhou G, Li Y, Zhang Y, Wang T, Yeung CKL, Chen L, Ma tion in the era of ancient genomics. Nat Rev Genet. 21(8):449–460. J, et al. 2013. Genomic analyses identify distinct patterns of selection Frichot E, Mathieu F, Trouillon T, Bouchard G, Francois O. 2014. Fast and in domesticated pigs and Tibetan wild boars. Nat Genet. efficient estimation of individual ancestry coefficients. Genetics 45(12):1431–1438. 196(4):973–983. Li M, Wu X, Guo X, Bao P, Ding X, Chu M, Liang C, Yan P. 2018. Galal S, Gürsoy O, Shaat I. 2008. Awassi sheep as a genetic resource and Comparative iTRAQ proteomics revealed proteins associated with efforts for their genetic improvement—A review. Small Rumin Res. horn development in yak. Proteome Sci. 16(1):1–11. 79(2–3):99–108. Li W, Lu ZF, Man XY, Li CM, Zhou J, Chen JQ, Yang XH, Wu XJ, Cai SQ, Galbraith H. 2010. Fundamental hair follicle biology and fine fibre pro- Zheng M. 2012. VEGF upregulates VEGF receptor-2 on human outer duction in animals. Animal 4(9):1490–1509. root sheath cells and stimulates proliferation through ERK pathway. Gignoux CR, Henn BM, Mountain JL. 2011. Rapid, global demographic Mol Biol Rep. 39(9):8687–8694. expansions after the origins of agriculture. Proc Natl Acad Sci U S A. Li W, Man XY, Li CM, Chen JQ, Zhou J, Cai SQ, Lu ZF, Zheng M. 108(15):6044–6049. 2012. VEGF induces proliferation of human hair follicle dermal pa- Green RE, Krause J, Briggs AW, Maricic T, Stenzel U, Kircher M, Patterson pilla cells through VEGFR-2-mediated activation of ERK. Exp Cell Res. N, Li H, Zhai W, Fritz MHY, et al. 2010. A draft sequence of the 318(14):1633–1640. neandertal genome. Science 328(5979):710–722. Li X, Yang J, Shen M, Xie X-L, Liu G-J, Xu Y-X, Lv F-H, Yang H, Yang Y-L, Hayward JJ, Castelhano MG, Oliveira KC, Corey E, Balkman C, Baxter TL, Liu C-B, et al. 2020. Whole-genome resequencing of wild and do- Casal ML, Center SA, Fang M, Garrison SJ, et al. 2016. Complex mestic sheep identifies genes associated with morphological and disease and phenotype mapping in the domestic dog. Nat agronomic traits. Nat Commun. 11(1):2815. Commun. 7(1):1–11. Lipson M, Szecsenyi-Nagy A, Mallick S, Posa A, Stegmar B, Keerl V, Heaton MP, Clawson ML, Chitko-Mckown CG, Leymaster KA, Smith Rohland N, Stewardson K, Ferry M, Michel M, et al. 2017. Parallel TPL, Harhay GP, White SN, Herrmann-Hoesing LM, Mousel MR, palaeogenomic transects reveal complex genetic history of early Lewis GS, et al. 2012. Reduced lentivirus susceptibility in sheep European farmers. Nature 551(7680):368–372. with TMEM154 mutations. PLoS Genet. 8(1):e1002467. Lorenzini R, Cabras P, Fanelli R, Carboni GL. 2011. Wildlife molecular Hill WG, Robertson A. 1968. Linkage disequilibrium in finite populations. forensics: identification of the Sardinian mouflon using STR profiling Theor Appl Genet. 38(6):226–231. and the Bayesian assignment test. Forensic Sci Int Genet. 5(4):345–349. Ho TV, Yap N, Almontashiri N, Stewart AF. 2012. Embryonic lethality of Lv FH, Peng WF, Yang J, Zhao YX, Li WR, Liu MJ, Ma YH, Zhao QJ, Yang mice deficient for IRF2BP2: an essential ischemia-inducible transcrip- GL, Wang F, et al. 2015. Mitogenomic meta-analysis identifies two tional coactivator of VEGFA. Circulation 126:A19581. phases of migration in the history of Eastern Eurasian sheep. Mol Biol Hu XJ, Yang J, Xie XL, Lv FH, Cao YH, Li WR, Liu MJ, Wang YT, Li JQ, Liu Evol. 32(10):2515–2533. YG, et al. 2019. The genome landscape of Tibetan sheep reveals Ma GW, Chu YK, Zhang WJ, Qin FY, Xu SS, Yang H, Rong EG, Du ZQ, adaptive introgression from argali and the history of early human Wang SZ, Li H, et al. 2017. Polymorphisms of FST gene and their settlements on the Qinghai-Tibetan Plateau. Mol Biol Evol. association with wool quality traits in Chinese Merino sheep. PLoS 36(2):283–303. One 12(4):e0174868. Huang DW, Sherman BT, Lempicki RA. 2009. Bioinformatics enrichment Mackowiak SD. 2011. Identification of novel and known miRNAs in tools: paths toward the comprehensive functional analysis of large deep-sequencing data with miRDeep2. Curr Protoc Bioinform. gene lists. Nucleic Acids Res. 37(1):1–13. 36:12–10. Huerta-Sanchez E, Jin X, Asan BZ, Peter BM, Vinckenbosch N, Liang Y, Yi Majumder PP. 2010. The human genetic history of South Asia. Curr Biol. X, He M, Somel M, et al. 2014. Altitude adaptation in Tibetans 20(4):R184–R187. caused by introgression of Denisovan-like DNA. Nature 512:194–197. Malaspinas AS, Westaway MC, Muller C, Sousa VC, Lao O, Alves I, Huson DH, Bryant D. 2006. Application of phylogenetic networks in Bergström A, Athanasiadis G, Cheng JY, Crawford JE, et al. 2016. A evolutionary studies. Mol Biol Evol. 23(2):254–267. genomic history of aboriginal Australia. Nature 538(7624):207–214. Jones MR, Mills LS, Alves PC, Callahan CM, Alves JM, Lafferty DJR, Jiggins Manichaikul A, Mychaleckyj JC, Rich SS, Daly K, Sale M, Chen WM. 2010. FM, Jensen JD, Melo-Ferreira J, Good JM. 2018. Adaptive introgres- Robust relationship inference in genome-wide association studies. sion underlies polymorphic seasonal camouflage in snowshoe hares. Bioinformatics 26(22):2867–2873. Science 360(6395):1355–1358. Martin M. 2011. Cutadapt removes adapter sequences from high- Kalendar R, Khassenov B, Ramankulov Y, Samuilova O, Ivanov KI, 2017. throughput sequencing reads. EMBnetjournal 17:10. FastPCR: an in silico tool for fast primer and probe design and ad- Martin SH, Davey JW, Jiggins CD. 2015. Evaluating the use of ABBA– vanced sequence analysis. Genomics 109(3–4):312–319.1. BABA statistics to locate introgressed loci. Mol Biol Evol. Kijas JW, Lenstra JA, Hayes B, Boitard S, Neto LRP, San Cristobal M, Servin 32(1):244–257. B, McCulloch R, Whan V, Gietzen K, et al. 2012. Genome-wide Martin SH, Möst M, Palmer WJ, Salazar C, McMillan WO, Jiggins FM, analysis of the world’s sheep breeds reveals high levels of historic Jiggins CD. 2016. Natural selection and genetic diversity in the but- mixture and strong recent selection. PLoS Biol. 10(2):e1001258. terfly heliconius melpomene. Genetics 203(1):525–541. Kim JY, Jeong S, Kim KH, Lim WJ, Lee HY, Kim N. 2019. Discovery of McColl H, Racimo F, Vinner L, Demeter F, Gakuhari T, Moreno-Mayar JV, genomic characteristics and selection signatures in Korean indigenous van Driem G, Gram Wilken U, Seguin-Orlando A, de la Fuente goats through comparison of 10 goat breeds. Front Genet. 10:699. Castro C, et al. 2018. The prehistoric peopling of Southeast Asia. Kozomara A, Birgaoanu M, Griffiths-Jones S. 2019. miRBase: from microRNA Science 361(6397):88–92. sequences to function. Nucleic Acids Res. 47(D1):D155–D162. McKenna A, Hanna M, Banks E, Sivachenko A, Cibulskis K, Kernytsky A, Lawal RA, Martin SH, Vanmechelen K, Vereijken A, Silva P, Al-Atiyat RM, Garimella K, Altshuler D, Gabriel S, Daly M, et al. 2010. The genome Aljumaah RS, Mwacharo JM, Wu DD, Zhang YP, et al. 2020. The wild analysis toolkit: a mapreduce framework for analyzing next- species genome ancestry of domestic chickens. BMC Biol. 18(1):13. generation DNA sequencing data. Genome Res. 20(9):1297–1303. Lazaridis I, Nadel D, Rollefson G, Merrett DC, Rohland N, Mallick S, McQuillan R, Leutenegger AL, Abdel-Rahman R, Franklin CS, Pericic M, Fernandes D, Novak M, Gamarra B, Sirak K, et al. 2016. Genomic Barac-Lauc L, Smolej-Narancic N, Janicijevic B, Polasek O, Tenesa A, insights into the origin of farming in the ancient Near East. Nature et al. 2008. Runs of homozygosity in European populations. Am J 536(7617):419–424. Hum Genet. 83(3):359–372. Li H, Durbin R. 2009. Fast and accurate short read alignment with Medugorac I, Graf A, Grohs C, Rothammer S, Zagdsuren Y, Gladyr E, Burrows–Wheeler transform. Bioinformatics 25(14):1754–1760. Zinovieva N, Barbieri J, Seichter D, Russ I, et al. 2017. Whole-genome Li H, Durbin R. 2011. Inference of human population history from indi- analysis of introgressive hybridization and characterization of the vidual whole-genome sequences. Nature 475(7357):493–496. bovine legacy of Mongolian yaks. Nat Genet. 49(3):470–475. 24 Genetic Diversity, Introgression, and Agronomically Important Loci of Worldwide Sheep . doi:10.1093/molbev/msab353 MBE Mei C, Wang H, Liao Q, Wang L, Cheng G, Wang H, Zhao C, Zhao S, Song Rimbault M, Beale HC, Schoenebeck JJ, Hoopes BC, Allen JJ, Kilroy-Glynn J, Guang X, et al. 2018. Genetic architecture and selection of Chinese P, Wayne RK, Sutter NB, Ostrander EA. 2013. Derived variants at six cattle revealed by whole genome resequencing. Mol Biol Evol. genes explain nearly half of size reduction in dog breeds. Genome Res. 35(3):688–699. 23(12):1985–1995. Miao B, Wang Z, Li Y. 2016. Genomic analysis reveals hypoxia adaptation Ryder ML. 1984. Sheep. In: Mason IL, editor. Evolution of domesticated in the Tibetan Mastiff by introgression of the gray wolf from the animals. London/New York: Longman. p. 63–85. Tibetan Plateau. Mol Biol Evol. 34:734–743. Saravanan KA, Panigrahi M, Kumar H, Bhushan B, Dutt T, Mishra BP. Muigai AT, Hanotte O. 2013. The origin of African sheep: archaeological 2020. Selection signatures in livestock genome: A review of concepts, and genetic perspectives. Afr Archaeol Rev. 30(1):39–50. approaches and applications. Livest Sci. 241:104257. Mu~ noz M, Bozzi R, Garcıa-Casco J, Nun~ez Y, Ribani A, Franci O, Garcıa F, Sayres MAW. 2018. Genetic diversity on the sex chromosomes. Genome Skrlep M, Schiavo G, Bovo S, et al. 2019. Genomic diversity, linkage Biol Evol. 10:1064–1078. disequilibrium and selection signatures in European local pig breeds Shaha R. 1970. Nepal, Tibet and China. J Nepal Council World Aff. 3:13–82. assessed with a high density SNP chip. Sci Rep. 9:1–14. Shinde V, Narasimhan VM, Rohland N, Mallick S, Mah M, Lipson M, Narasimhan VM, Patterson N, Moorjani P, Rohland N, Bernardos R, Nakatsuka N, Adamski N, Broomandkhoshbacht N, Ferry M, et al. Mallick S, Lazaridis I, Nakatsuka N, Olalde I, Lipson M, et al. 2019. 2019. An ancient harappan genome lacks ancestry from steppe The formation of human populations in South and Central Asia. pastoralists or Iranian farmers. Cell 179(3):729–735. Science 365(6457):eaat7487. Singh S, Kumar S Jr, Kolte AP, Kumar S. 2013. Extensive variation and Naval-Sanchez M, Nguyen Q, McWilliam S, Porto-Neto LR, Tellam R, sub-structuring in lineage A mtDNA in Indian sheep: genetic evi- Vuocolo T, Reverter A, Perez-Enciso M, Brauning R, Clarke S, et al. dence for domestication of sheep in India. PLoS One 8(11):e77858. 2018. Sheep genome functional annotation reveals proximal regula- Skourtanioti E, Erdal YS, Frangipane M, Balossi Restelli F, Yener KA, tory elements contributed to the evolution of modern breeds. Nat € Pinnock F, Matthiae P, Ozbal R, Schoop UD, Guliyev F, et al. 2020. Commun. 9(1):859. Genomic history of Neolithic to Bronze Age Anatolia, Northern Nkedianye DK, Ogutu JO, Said MY, Kifugo S, de Leeuw J, Van Gardingen Levant, and Southern Caucasus. Cell 181(5):1158–1175. P, Reid RS. 2019. Livestock-wealth inequalities and uptake of crop Spangler GL, Rosen BD, Ilori MB, Hanotte O, Kim E-S, Sonstegard TS, cultivation among the Maasai of Kenya and Tanzania. World Dev Burke JM, Morgan JLM, Notter DR, Van Tassell CP. 2017. Whole Perspect. 14:100106. genome structural analysis of Caribbean hair sheep reveals quanti- Paim TP, Paiva SR, Toledo NM, Yamaghishi MB, Carneiro PLS, Faco O, tative link to West African ancestry. PLoS One 12(6):e0179021. Ara ujo AM, Azevedo HC, Caetano AR, Braga RM, et al. 2021. Origin Sved JA. 1971. Linkage disequilibrium and homozygosity of chromosome and population structure of Brazilian hair sheep breeds. Anim Genet. segments in finite populations. Theor Popul Biol. 2(2):125–141. 52(4):492–504. Tapio M, Marzanov N, Ozerov M, Cinkulov M, Gonzarenko G, Kiselyova Pan Z, Li S, Liu Q, Wang Z, Zhou Z, Di R, Miao B, Hu W, Wang X, Hu X, et T, Murawski M, Viinalass H, Kantanen J. 2006. Sheep mitochondrial al. 2018. Whole-genome sequences of 89 Chinese sheep suggest role DNA variation in European, Caucasian, and Central Asian areas. Mol of RXFP2 in the development of unique horn phenotype as response Biol Evol. 23(9):1776–1783. to semi-feralization. Gigascience 7(4):giy019. Terhorst J, Kamm JA, Song YS. 2017. Robust and scalable inference of Patterson N, Price AL, Reich D. 2006. Population structure and eigena- population history from hundreds of unphased whole genomes. Nat nalysis. PLoS Genet. 2(12):e190. Genet. 49(2):303–309. Patterson NJ, Moorjani P, Luo Y, Mallick S, Rohland N, Zhan Y, Thorne JW, Murdoch BM, Freking BA, Redden RR, Murphy TW, Taylor Genschoreck T, Webster T, Reich D. 2012. Ancient admixture in JB, Blackburn HD. 2021. Evolution of the sheep industry and genetic human history. Genetics 192(3):1065–1093. research in the United States: opportunities for convergence in the Pavlidis P, Z ivkovic D, Stamatakis A, Alachiotis N. 2013. SweeD: twenty-first century. Anim Genet. 52(4):395–408. likelihood-based detection of selective sweeps in thousands of Trapnell C, Roberts A, Goff L, Pertea G, Kim D, Kelley DR, Pimentel H, genomes. Mol Biol Evol. 30(9):2224–2234. Salzberg SL, Rinn JL, Pachter L. 2012. Differential gene and transcript Pe H. 1959. G. H. Luce and Pe Maung Tin (ed.): Inscriptions of Burma. expression analysis of RNA-seq experiments with TopHat and Portfolio IV. Down to 702 B.E. (1340 A.D.).—Portfolio V. 703–726 B.E. Cufflinks. Nat Protoc. 7(3):562–578. (1341–1364 A.D.). (University of Rangoon Oriental Studies Publication Utsunomiya YT, Perez O’Brien AM, Sonstegard TS, Sölkner J, Garcia JF. No. 5, No. 6.) 63 pp., plates 346–462; 38 pp., plates 463–609. Oxford: 2015. Genomic data as the “hitchhiker’s guide” to cattle adaptation: University Press, 1956. Bull Sch Orient Afr Stud. 22:177. tracking the milestones of past selection in the bovine genome. Front Petit M, Astruc JM, Sarry J, Drouilhet L, Fabre S, Moreno CR, Servin B. Genet. 6:36. 2017. Variation in recombination rate and its genetic determinism in Vigne JD. 1992. Zooarchaeology and the biogeographical history of the sheep populations. Genetics 207(2):767–784. mammals of Corsica and Sardinia since the Last Ice-Age. Mamm Rev. Purcell S, Neale B, Todd-Brown K, Thomas L, Ferreira MAR, Bender D, 22(2):87–96. Maller J, Sklar P, de Bakker PIW, Daly MJ, et al. 2007. PLINK: a tool set Wang GD, Zhai W, Yang HC, Wang L, Zhong L, Liu YH, Fan RX, Yin TT, for whole-genome association and population-based linkage analy- Zhu CL, Poyarkov AD. 2015. Out of southern East Asia: the natural ses. Am J Hum Genet. 81(3):559–575. history of domestic dogs across the world. Cell Res. 47:217–225. R Core Team. 2020. R: a language and environment for statistical com- Wang MS, Wang S, Li Y, Jhala Y, Thakur M, Otecko NO, Si J-F, Chen H-M, puting. R Foundation for Statistical Computing, Vienna, Austria. Shapiro B, Nielsen R, et al. 2020. Ancient hybridization with an un- Raghavan M, Schroeder H, Malaspinas AS. 2019. An ancient genome known population facilitated high altitude adaptation of canids. Mol from the Indus Valley civilization. Cell 179(3):586–588. Biol Evol. 37(9):2616–2629. Rausch T, Zichner T, Schlattl A, Stutz AM, Benes V, Korbel JO. 2012. Wang Y, Cao Z, Ogilvie HA, Nakhleh L. 2021. Phylogenomic assessment DELLY: structural variant discovery by integrated paired-end and of the role of hybridization and introgression in trait evolution. PLoS split-read analysis. Bioinformatics 28(18):i333–i339. Genet. 17(8):e1009701. Revelo HA, Lopez-Alvarez D, Landi V, Rizzo L, Alvarez LA. 2020. Weir BS, Cockerham CC. 1984. Estimating F-statistics for the analysis of Mitochondrial DNA variations in Colombian Creole sheep confirm population-structure. Evolution 38(6):1358–1370. an Iberian origin and shed light on the dynamics of introduction Wu DD, Ding XD, Wang S, Wojcik JM, Zhang Y, Tokarska M, Li Y, Wang events of African genotypes. Animals 10(9):1594. MS, Faruque O, Nielsen R, et al. 2018. Pervasive introgression facil- Rezaei HR, Naderi S, Chintauan-Marquier IC, Taberlet P, Virk AT, itated domestication and adaptation in the Bos species complex. Nat Naghash HR, Rioux D, Kaboli M, Pompanon F. 2010. Evolution Ecol Evol. 2(7):1139–1145. and taxonomy of the wild species of the genus Ovis (Mammalia, Wu J, Wang L, Fu J, Chen J, Wei S, Zhang S, Zhang J, Tang Y, Chen M, Zhu Artiodactyla, Bovidae). Mol Phylogenet Evol. 54(2):315–326. J, et al. 2020. Resequencing of 683 common bean genotypes 25 Lv et al. . doi:10.1093/molbev/msab353 MBE identifies yield component trait associations across a north–south Zhang C, Dong SS, Xu JY, He WM, Yang TL. 2019. PopLDdecay: a fast and cline. Nat Genet. 52(1):118–125. effective tool for linkage disequilibrium decay analysis based on var- Yang J, Lee SH, Goddard ME, Visscher PM. 2011. GCTA: a tool for iant call format files. Bioinformatics 35(10):1786–1788. genome-wide complex trait analysis. Am J Hum Genet. 88:76–82. Zhang JJ, Wang DG, Zhou XB, Gao Y, He XL, Chen YL, Zhang EP. Yang J, Li WR, Lv FH, He SG, Tian SL, Peng WF, Sun YW, Zhao YX, Tu XL, 2018. Effect of VEGF on secondary hair follicle outer root Zhang M, et al. 2016. Whole-genome sequencing of native sheep sheath cells of cashmere goat in vitro. Acta Vet Zootech Sin. provides insights into rapid adaptations to extreme environments. 49:1124–1133. Mol Biol Evol. 33(10):2576–2592. Zhao YX, Yang J, Lv FH, Hu XJ, Xie XL, Zhang M, Li WR, Liu MJ, Wang YT, Yi X, Liang Y, Huerta-Sanchez E, Jin X, Cuo ZXP, Pool JE, Xu X, Jiang H, Li JQ, et al. 2017. Genomic reconstruction of the history Vinckenbosch N, Korneliussen TS, et al. 2010. Sequencing of 50 hu- of native sheep reveals the peopling patterns of nomads and the man exomes reveals adaptation to high altitude. Science 329:75–78. expansion of early pastoralism in East Asia. Mol Biol Evol. Yu H, Xing YT, Meng H, He B, Li WJ, Qi XZ, Zhao JY, Zhuang Y, Xu X, 34(9):2380–2395. Yamaguchi N, et al. 2021. Genomic evidence for the Chinese moun- Zheng Z, Wang X, Li M, Li Y, Yang Z, Wang X, Pan X, Gong M, Zhang Y, tain cat as a wildcat conspecific (Felis silvestris bieti) and its intro- Guo Y, et al. 2020. The origin of domestication genes in goats. Sci gression to domestic cats. Sci Adv. 7(26):eabg0221. Adv. 6(21):eaaz5216. Zeder MA. 2008. Domestication and early agriculture in the Zhou Y, Zhao X, Li Y, Xu J, Bi A, Kang L, Xu D, Chen H, Wang Y, Wang YG, Mediterranean Basin: origins, diffusion, and impact. Proc Natl Acad et al. 2020. Triticum population sequencing provides insights into Sci U S A. 105(33):11597–11604. wheat adaptation. Nat Genet. 52(12):1412–1422. 26