Mitogenomic meta-analysis identifies two phases of migration in the history of Eastern Eurasian sheep

Article Fast Track Mitogenomic Meta-Analysis Identifies Two Phases of Migration in the History of Eastern Eurasian Sheep Feng-Hua Lv, y,1 Wei-Feng Peng, y,1,2 Ji Yang, y,1 Yong-Xin Zhao, 1,2 Wen-Rong Li, 3 Ming-Jun Liu, 3 Yue-Hui Ma, 4 Qian-Jun Zhao, 4 Guang-Li Yang, 1,5 Feng Wang, 6 Jin-Quan Li, 7 Yong-Gang Liu, 8 Zhi-Qiang Shen, 9 Sheng-Guo Zhao, 10 EEr Hehua, 11 Neena A. Gorkhali, 4,12 S. M. Farhad Vahidi, 13 Muhammad Muladno, 14 Arifa N. Naqvi, 15 Jonna Tabell, 16 Terhi Iso-Touru, 16 Michael W. Bruford, 17 Juha Kantanen, 16,18 Jian-Lin Han,* ,4,19 Meng-Hua Li* ,1 1 CAS Key Laboratory of Animal Ecology and Conservation Biology, Institute of Zoology, Chinese Academy of Sciences (CAS), Beijing, China 2 University of Chinese Academy of Sciences (UCAS), Beijing, China 3 Animal Biotechnology Research Institute, Xinjiang Academy of Animal Science, Urumqi, China 4 CAAS-ILRI Joint Laboratory on Livestock and Forage Genetic Resources, Institute of Animal Science, Chinese Academy of Agricultural Sciences (CAAS), Beijing, China 5 College of Life Sciences, Shangqiu Normal University, Shangqiu, China 6 Institute of Sheep and Goat Science, Nanjing Agricultural University, Nanjing, China 7 College of Animal Science, Inner Mongolia Agricultural University, Hohhot, China 8 College of Animal Science and Technology, Yunnan Agricultural University, Kunming, China 9 Shandong Binzhou Academy of Animal Science and Veterinary Medicine, Binzhou, China 10 College of Animal Science and Technology, Gansu Agricultural University, Lanzhou, China 11 Grass-Feeding Livestock Engineering Technology Research Center, Ningxia Academy of Agriculture and Forestry Sciences, Yinchuan, China 12 Animal Breeding Division, National Animal Science Institute, Nepal Agriculture Research Council, Kathmandu, Nepal 13 Agricultural Biotechnology Research Institute of Iran-North Branch (ABRII), Rasht, Iran 14 Department of Animal Technology and Production Science, Bogor Agricultural University, Darmaga Campus, Bogor, Indonesia 15 Faculty of Life Sciences, Karakoram International University, Gilgit, Baltistan, Pakistan 16 Green Technology, Natural Resources Institute Finland (LUKE), Jokioinen, Finland 17 School of Biosciences and Sustainable Places Research Institute, Cardiff University, Cardiff, United Kingdom 18 Department of Biology, University of Eastern Finland, Kuopio, Finland 19 International Livestock Research Institute (ILRI), Nairobi, Kenya y These authors contributed equally to this work. *Corresponding author: E-mail: [email protected]; m[email protected]. Associate editor: David Irwin Abstract Despite much attention, history of sheep (Ovis aries) evolution, including its dating, demographic trajectory and geographic spread, remains controversial. To address these questions, we generated 45 complete and 875 partial mitogenomic sequences, and performed a meta-analysis of these and published ovine mitochondrial DNA sequences (n= 3,229) across Eurasia. We inferred that O. orientalis and O. musimon share the most recent female ancestor with O. aries at approximately 0.790 Ma (95% CI: 0.637–0.934 Ma) during the Middle Pleistocene, substantially predating the domestication event (~8–11 ka). By reconstructing historical variations in effective population size, we found evidence of a rapid population increase approximately 20–60 ka, immediately before the Last Glacial Maximum. Analyses of lineage expansions showed two sheep migratory waves at approximately 4.5–6.8 ka (lineages A and B: ~6.4–6.8 ka; C: ~4.5 ka) across eastern Eurasia, which could have been influenced by prehistoric West–East commercial trade and deliberate mating of domestic and wild sheep, respectively. A continent-scale examination of lineage diversity and approximate Bayesian computation analyses indicated that the Mongolian Plateau region was a secondary center of dispersal, acting as a “transportation hub” in eastern Eurasia: SheepfromtheMiddleEasterndomesticationcenterwere inferred to have migrated through the Caucasus and Central Asia, and arrived in North and Southwest China (lineages A, B, and C) and the Indian subcontinent ßThe Author 2015. 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 Non-Commercial License (http://creativecommons.org/licenses/by-nc/4.0/), which permits non-commercial re-use, distribution, and reproduction in any medium, provided the original work is properly cited. For commercial re-use, please contact journals.permission[email protected] Open Access Mol. Biol. Evol. 32(10):2515–2533 doi:10.1093/molbev/msv139 Advance Access publication June 16, 2015 2515 at Natural Resources Institute Finland (Luke) on October 6, 2016http://mbe.oxfordjournals.org/Downloaded from (lineages B and C) through this region. Our results provide new insights into sheep domestication, particularly with respect to origins and migrations to and from eastern Eurasia. Key words: wild ancestor, domestication, gene flow, mitogenome, Ovis aries, meta-analysis, colonization simulation. Introduction As one of the first animals ever domesticated, sheep (Ovis aries) have played an important role in human society and have spread almost globally, following human migrations (Colledge et al. 2005;Chessa et al. 2009). Early evidence implied that modern sheep breeds were first domesticated from Asian mouflon (O. orientalis) in the Fertile Crescent approximately 8–11 thousand years ago (ka) (Ryder 1984). Following domestication, as many as 1,400 sheep breeds have been developed from their wild ancestors after long-term natural and intense artificial selection (Scherf 2000). During this process, human activities have played a significant role in determining the patterns of gene flow among breeds and populations (e.g., Warmuth et al. 2012). Thus, an examination of continentwide genetic variability among modern native sheep breeds can provide a comprehensive, in-depth understanding of their genetic origins and dispersal, as well as insight into the impact of human activities on sheep throughout history. In recent decades, remarkable analytical advances in paleontological and molecular genetics have transformed our understanding of the origins and regional expansion of domestic sheep (Poplin 1979;Hiendleder, Mainz, et al. 1998;Pedrosa et al. 2005;Chessa et al. 2009;Meadows et al. 2007,2011;Kijas et al. 2009,2012;Demirci et al. 2013). Morphological change and demographic analysis implied that sheep were likely brought under domestication in a region that stretches from northern Zagros to southeastern Anatolia, approximately 10.5–11 ka or perhaps even earlier (Peters et al. 2005). In addition, a recent investigation on endogenous retroviral sequences revealed a remarkable secondary population expansion of improved domestic sheep, most likely out of Southwest Asia (i.e., the Middle East; Chessa et al. 2009). Mitochondrial DNA (mtDNA) sequence analyses have identified a general phenomenon of multiple maternal lineages (i.e., A, B, C, D, and E), some with specific geographic ranges, implying multiple maternal origins and possibly independent domestication events in sheep (Wood and Phua 1996; Hiendleder, Mainz, et al. 1998;Guo et al. 2005;Pedrosa et al. 2005;Tapio et al. 2006;Meadows et al. 2007;Singh et al. 2013). Estimates from complete and/or partial mtDNA sequences have enabled various divergence time estimates between domestic and wild sheep as well as among the five major maternal lineages of O. aries (e.g., Hiendleder, Mainz, et al. 1998; Pedrosa et al. 2005;Chen et al. 2006;Meadows et al. 2011). In general, the estimated divergence times among the five major lineages have been much earlier than the domestication period inferred from archeological evidence (Bar-Yosef and Meadow 1995;Zeder 2008). For example, the divergence time between the two most common lineages (i.e., A and B) was estimated to be as early as 1.6–1.7 Ma based on cytochrome b (Cyt-b) sequences (Hiendleder, Mainz, et al. 1998). In addition, Pedrosa et al. (2005) and Chen et al. (2006) suggested the divergence time of lineage C from lineages A and B to be approximately 0.42–0.76 Ma and approximately 0.45–0.75 Ma from the analysis of control region and Cyt-bsequences, respectively. However, a more recent study (Meadows et al. 2011) using 12 protein-coding genes from complete mitogenomes implied more recent divergence between the lineages: For example, 0.590 0.17 Ma between A and B and 0.26 0.09 Ma between C and E. So far, most ovine mtDNA investigations have only focused on one or two segments within Cyt-bgene and the control region (including the hypervariable region; e.g., Pedrosa et al. 2005); nevertheless, high levels of recurrent mutations observed in the short segment within control region in many mammal species may bias dating estimates (e.g., Achilli et al. 2009,2012; see also the reviews in Torroni et al. 2006;Taberlet et al. 2008). Moreover, previous sheep mtDNA studies have merely included breeds at a regional (e.g., Pedrosa et al. 2005;Chen et al. 2006;Wang et al. 2006; Meadows et al. 2007)orsubcontinentalscale(e.g.,Tapio et al. 2006), whereas maternal lineages of domestic sheep, particularly for breeds in Southwest, Central, East and South Asia, including the Caucasus, Iran, Pakistan, Nepal, Indonesia, Mongolia, China, and India, have been largely excluded from integrated analyses. In addition, the divergence scenarios have not been fully evaluated based on complete mitogenomes either, which could have provided refined phylogenies of maternal lineages and robust estimations of genetic variability and divergence time in domestic animals (see the review in Wang et al. 2014). Therefore, although these early mtDNA studies have provided useful insights into the history of sheep domestication in Eurasia, answers to some basic questions surrounding the domestication process are far from being settled. For example, phylogenetic relationships among wild and domestic sheep (e.g., Hiendleder, Lewalski, et al. 1998;Meadows et al. 2007), divergence times between the major maternal lineages (e.g., Pedrosa et al. 2005;Zeder 2008;Meadows et al. 2011), demographic history and population recolonization (Dobney and Larson 2006;Zeder 2008), and origins of different mtDNA lineages (Tapio et al. 2006; Meadows et al. 2007;Demirci et al. 2013;Singh et al. 2013), as well as the continent-wide patterns of gene flow from the postulated Middle Eastern domestication center to Central, East and South Asia (see, e.g., Tapio et al. 2006,2010;Cai et al. 2007,2011) remain provisional or unaddressed. The main objective of our study was to better understand the domestication and expansion of O. aries across Eurasia through a meta-analysis of complete and partial ovine mitogenomic sequences. More specifically, we aimed to refine and challenge existing paradigms on the wild origin, lineage divergence, demographic history and population recolonization of modern sheep, particularly the breeds present in eastern Eurasia. For these purposes, we sequenced the complete 2516 Lv et al. .doi:10.1093/molbev/msv139 MBE at Natural Resources Institute Finland (Luke) on October 6, 2016http://mbe.oxfordjournals.org/Downloaded from mitogenomes of 45 individuals (including O. orientalis,O. vignei, and 42 native breeds of O. aries) and the control region of a total of 875 animals (including 51 native breeds) from eastern Eurasia (fig. 1 and supplementary tables S1 and S2,Supplementary Material online). Together with the sequences retrieved from GenBank, we analyzed 85 complete mitogenomes of domestic sheep including each of the 5 lineages and 10 complete mitogenomes of O. orientalis,O. musimon,O. vignei,O. ammon,andO. canadensis using phylogenetics, molecular-dating, and demographic-reconstruction approaches. Full control region and Cyt-bsequences of seven extant wild sheep species (O. orientalis,O. musimon,O. vignei,O. ammon,O. canadensis,O. dalli,andO. nivicola)were also included in phylogenetic reconstructions. Furthermore, we carried out a meta-analysis and a simulation of colonization (e.g., approximate Bayesian computation, ABC) of mtDNA sequences, including 547 partial Cyt-band 1,470 partial control region sequences published previously (supplementary tables S2 and S3,Supplementary Material online), from native sheep breeds across Eurasia. We tried to address these questions and test two hypotheses on domestication and migrations of sheep distributed particularly in eastern Eurasia. One is the more recent origin and dispersal of lineage C when compared with those of the two widely distributed lineages A and B (Bruford 2005;Tapio et al. 2006). Another is that the arrival of some Indian sheep from the Middle Eastern domestication center could be through the Mongolian Plateau region, where archeological remains showed an early presence of domestic sheep in ancient history (e.g., Kuo et al. 1999;seealsoYang et al. 2015). Our results could help researchers better understand the demographic forces and human practice associated with animal domestication and migration in history (e.g., Hodges 1999;Larson et al. 2007,2010;Larson and Burger 2013). FIG.1. Geographic distribution of the samples in this and early ovine mtDNA studies. 2517 Ovine Mitogenomic Variations across Eastern Eurasia .doi:10.1093/molbev/msv139 MBE at Natural Resources Institute Finland (Luke) on October 6, 2016http://mbe.oxfordjournals.org/Downloaded from Results Geographic Patterns of mtDNA Variation The 45 complete domestic (GenBank accession numbers KF938317–KF938359) and wild (KF938360–KF938361) sheep mitogenomes (supplementary table S1, Supplementary Material online) sequenced in this study showed considerable sequence variability as well as variation in diversity among different regions (supplementary table S4 and fig. S1, Supplementary Material online). Also, we detected a large number of variable sites in the integrated data of partial Cyt-band control region (supplementary tables S2 and S3,Supplementary Material online). Full description of the complete mitogenome and partial mtDNA sequence variationsisinsupplementary information S1,Supplementary Material online. All control region and Cyt-bsequences analyzed in this study can be assigned to the five previously defined lineages (supplementary tables S2 and S3; see also supplementary figs. S2 and S3,Supplementary Material online). The two partial mtDNA fragments displayed similar geographic patterns (fig. 2Band C). For control region sequences, lineages A and FIG.2.Geographic distribution of the five major maternal lineages across Eurasia based on sequences obtained in this study and retrieved from GenBank. (A) Phylogenetic tree inferred from partial control region sequences (left) and lineage composition of sheep in different geographic regions at different time points (right) based on ancient specimens (Cai et al. 2007,2011;Demirci et al. 2013;Niemi et al. 2013); (B) lineage frequency distribution of partial control region sequences; previously reported lineage frequencies in 12 regions (I–XII) are detailed in supplementary table S17,Supplementary Material online; (C) lineage frequency distribution of partial Cyt-bsequences; (D) geographic distribution of fat-tailed native sheep breeds (regions with black lines) and lineage C (region colored in purple). Pie plots show the proportions of the five distinct lineages (A–E) of domestic sheep in the different geographic regions (for the details of the geographic regions, see supplementary tables S5 and S6,Supplementary Material online). In the phylogenetic tree, diagnostic mutations are showed on the branches and are named according to their nucleotide positions relative to the reference sequence AF010406; amino acid replacements are underlined and synonymous replacements are marked in black. Control region mutations (15,437–16,616 bp) are shown in blue. Insertions are indicated by a “+” after the position number and followed by the type of inserted nucleotide(s). Mutations with prefix “” indicate identical variable sites found in Meadows et al. (2007), which are used to define the five major lineages. 2518 Lv et al. .doi:10.1093/molbev/msv139 MBE at Natural Resources Institute Finland (Luke) on October 6, 2016http://mbe.oxfordjournals.org/Downloaded from B were the most common and most widely distributed, with a mean combined frequency of approximately 89% (fig. 2Band C). Lineage A was extremely frequent (~77%) in the Indian subcontinent, although its frequency was less than 10% in Europe. In contrast, lineage B was found mostly in Europe, with its highest frequency ( 490%) in Southwest Europe (fig. 2Band C). Lineage C occurred mainly in the Middle East, the Caspian Sea region, North China, and the Mongolian Plateau, with a mean frequency of approximately 18% (fig. 2B), whereas a few haplotypes of lineage C were also found in the Iberian Peninsula, India, Nepal, and Southwest China. A majority of the breeds harboring lineage C were fattailed (including fat-rump; 73.1%), higher than the proportion of fat-tailed breeds having lineage A (50.8%) or B (44.8%) (supplementary tables S5 and S6,Supplementary Material online). In addition, we found a significantly higher mean frequency of lineage C in fat-tailed breeds than in shorttailed breeds (fat-tailed: f C =0–0.50, mean f C = 0.19; shorttailed: f C = 0–0.40, mean f C = 0.05; two-sample Kolmogorov– Smirnov test: P<0.01; supplementary fig. S4,Supplementary Material online). Of the total 149 breeds studied here, 66 are fat-tailed, 78 harbor lineage C, and 57 are fat-tailed sheep carrying lineage C. Compared with the overlap expected by chance, there is a large and significant excess of breeds that are fat-tailed harboring linages C (lineage C: observed n= 57, expected by chance n= 34.65, P<0.001; supplementary fig. S5;Supplementary Material online). Lineages D and E accounted for approximately 1% of the total samples and were only found in the Middle East (see fig. 2Band C). A synthetic map across Eurasia showed that the breeds in the Mongolian Plateau region had the highest genetic variability () of control region in Asia (fig. 3A;supplementary table S7,Supplementary Material online). For lineages A and B, a relatively high level of nucleotide diversity was found in the Indian subcontinent (fig. 3Band C). In addition, the synthetic map revealed the highest level of lineage C variability in the breeds of North China, even higher than that of the breedsintheMiddleEast(fig. 3D), the presumed domestication center of modern sheep (Ryder 1984). Phylogenetic Relationships Phylogenetic relationships inferred from all the 95 complete Ovis mitogenomes (supplementary table S1,Supplementary Material online) are shown in figure 4.The85complete mitogenomes of O. aries were assigned to five major lineages (fig. 4). Ovis vignei,O. ammon, and O. canadensis clustered into three independent clades separated from O. aries, whereas O. canadensis showed the largest divergence. The clade of O. musimon and O. orientalis was closely related to O. aries. In the phylogenetic trees built from the full control region, Cyt-b, and protein gene sequences of the complete mitogenomes, the four branches of wild sheep agreed with the topology inferred from the complete mitogenomes, but domestic sheep sequences formed an unresolved group rather than the five major lineages (supplementary figs. S6–S9,Supplementary Material online). Additional phylogenetic trees obtained with the full control region and Cyt-b sequences of wild and domestic sheep (supplementary figs. FIG.3. Synthetic maps illustrating geographic variation of nucleotide variability for the total lineages and lineages A, B, and C. (A) The total lineages, (B) lineage A, (C) lineage B, and (D)lineageC. 2519 Ovine Mitogenomic Variations across Eastern Eurasia .doi:10.1093/molbev/msv139 MBE at Natural Resources Institute Finland (Luke) on October 6, 2016http://mbe.oxfordjournals.org/Downloaded from S10 and S11 and tables S8 and S9,Supplementary Material online) showed different topologies from that inferred from the complete mitogenomes (fig. 4 and supplementary fig. S6, Supplementary Material online). Specifically, instead of showing close relationships only to lineage B as inferred from the complete mitogenomes (fig. 4), the haplotypes of O. musimon and O. orientalis clustered with lineages A, B, and CofO. aries control region sequences (supplementary fig. S10, AF010406 KF938353 EF490452 EF490456 KF938341 KF938346 KF938358 KF977845 HM236176-HM236177 KF302461 KF302462 KF302460 KF302447 KF302448 KF938340 KF938351 KF302452 KF302453 KF977846 EF490451 KF938355-KF938356 KF938350 KF938347 KF938352 KF938348 KF938357 EF490455 KF302450-KF302451 KF302449 KF302454 KF302455 KF302456-KF302457 KF302458 KF302459 KF938354 KF938343 KF938344 EF490453 EF490454 KF938360 HM236184 HM236185 KF938328 KF938329 KF938339 KF938349 KF938359 KF938333 KF938335 KF938325 KF938334 KF938321 KF938324 KF938322 KF938323 KF938319 KF938317 KF938337 KF938326 KF938345 KF977847 KF302440-KF302444 KF302445 KF302446 KF938330 HM236175 KF938332 KF938336 KF938338 KF938342 KF938331 HM236174 HM236180 HM236181 KF938320 KF938318 KF938327 HM236178 HM236179 HM236182 HM236183 HM236186 HM236187 HM236189 KF938361 HM236188 JX101654 JN181255 \\ \\ \\ 0.52 (0.346-0.694) 0.69 (0.494-0.887) 0.80 (0.5831.018) 0.31 (0.200-0.418) 2.60 2.93 (2.453-3.413) 8.31 (6.182-10.436) B A D E C A’B AB’D C’E ABD’CE A1 7777 VIVIVIIIII I PLEISTOCENE PLIOCENE NEOGENE QUATERNARY I: The late MIOCENE; II: ZANCLEAN; III: PIACENZIAN; IV: GELASIAN; V: CALABRIAN; VI: IONIAN 1.00 100 1.00 73 1.00 100 1.00 100 1.00 77 1.00 97 A B C D E Ovis orientalis Ovis Ovis ammon musimon Ovis vignei Ovis canadensis O vis ammo n O v i s v i g i n g g e i O vis c a n a d e d d n Oi i i n A2 A1a A1b B2 B1b B1 B1a B1a1 B1a2 B1a3 B1a4 B1a5 B1a6 B1a7 B1a8 B1a9 B1a10 B1a12 B1a11 FIG.4. Phylogeny of domestic and wild sheep inferred from a total of 95 complete mitogenomes (supplementary table S1,Supplementary Material online) using BI and ML methods with posterior probability (the first value) and bootstrap values (the second value) on the nodes, respectively. Divergence times for the lineages (Ma) were estimated only based on the 61 complete mitogenomes of native domestic sheep breeds and wild sheep species (see supplementary table S1,Supplementary Material online). 2520 Lv et al. .doi:10.1093/molbev/msv139 MBE at Natural Resources Institute Finland (Luke) on October 6, 2016http://mbe.oxfordjournals.org/Downloaded from Supplementary Material online), and they even shared some Cyt-bhaplotypes of lineages A, B, C, and E (supplementary fig. S11,Supplementary Material online). The reduced median network analysis of partial control region sequences showed several major radiating nodes at a few mutation steps within lineages A and B. Different contributions of breeds to different regions were evident, but none of the major nodes consisted of apparent region-specific haplotypes (supplementary fig. S2,Supplementary Material online). In addition, analysis of molecular variance and pairwise-population F ST values indicated genetic differentiation between European and Asian breeds, whereas considerable maternal gene flow was found among the breeds within Asia and Europe, respectively (supplementary figs. S12–S13 and tables S10 and S11,Supplementary Material online). Selective Pressure on Different Lineages The log-likelihood values (ln L) under the one-, two-, threeand four-ratio models were ln L=18,462.81, 18,454.79, 18,419.46 and 18,413.27, respectively (table 1). The ! ratio differed between the branches under the same model and varied for the same branches under different models (table 1). The likelihood ratio tests (LRTs) revealed that the differences between two models for all the pairwise comparisons were significant (P<0.01) and that the fourratio model (free-ratio model)bestfitthedata,which indicated different !ratios among the lineages. Mean ! values for the lineages were ! A = 0.0457 , ! B = 0.0775, ! D = 0.0494, and ! C+E = 0.0496 (supplementary fig. S14, Supplementary Material online); note that these values are all much lower than 1. This observation indicates that the maternal lineages (A, B, D, and C + E) have been under strong but variable intensity of purifying selection: Purifying selection on amino acid changes in lineage B has been slightly weaker than that on the other lineages. Thus, divergence time estimation (see below) based on the protein-coding genes would be biased. Instead, using the synonymous sites might be a better choice for divergence time estimation. Divergence Times for the Nodes The estimated divergence times within the comprehensive evolutionary framework of the Cetartiodactyla are shown in supplementary figure S15,Supplementary Material online. The O. vignei/O. aries split, which is the calibration point applied to estimate the divergence times between extant O. aries lineages, was 2.6 0.9 Ma. That time is far earlier than the most recent common ancestor (TMRCA) of domestic sheep (~0.79 Ma; 95% CI: 0.64– 0.93 Ma; table 2), and even older than the O. ammon/O. aries split (2.13 0.29 Ma) estimated by Meadows et al. (2011).TheCapra/Ovis split was estimated to be 14.7 2.1 Ma (supplementary fig. S15,Supplementary Material online), and is much older than the date based on the ungulate fossil record (~5.00–7.0 Ma; Luikart et al. 2001). Using the calibration point, we obtained a substitution rate of 0.70 10 8 substitutions per nucleotide/ year for complete mitogenome, 3.12 10 8 substitutions per nucleotide/year for control region, and 0.49 10 8 per nucleotide/year for Cyt-bwithout partitions. The divergence times for each node were mostly concordant under global and local clock models when estimated from the complete mitogenomes, the synonymous mutations or the third-codon positions (table 2). The earliest split was estimated to be approximately 0.73–0.93 Ma for thedivergenceofCandEfromA,B,andD(seethenode4 in table 2), whereas the most recent split was between lineages C and E at approximately 0.29–0.36 Ma (see the node 1 in table 2), greatly predating sheep domestication (~8–11 ka; Ryder 1984).ThetimetoTMRCAofthetwomostcommon lineages (A and B) was estimated to be approximately 0.50– 0.53 Ma (see node 2 in table 2). Under the relaxed molecular clock, we also obtained similar estimates of divergence times for the nodes based on different data sets (i.e., complete mitogenomes, synonymous, and third-codon positions; see table 2). However, divergence times for nodes 2, 3, and 4 estimated from synonymous mutations by the Bayesian Inference (BI) approach were significantly (P<0.05) higher than those by the global and local maximum likelihood (ML) Table 1. Number of Parameters Fitted, dN/dSRatios, Log-Likelihood Scores, and Their Differences under Different Models. Model pln LxModels Compared 2ln L A: One xratio x 0 102 18,462.81 x 0 = 0.0563 B: Two xratios x B 103 18,454.79 x B = 0.0702 x 0 x 0 = 0.0436 A versus B 16.04** C: Three xratios x A 104 18,419.46 x A = 0.0447 A versus C 86.70** x B x B = 0.0744 A versus D 99.08** x 0 x 0 = 0.0486 B versus C 70.66** D: Four xratios x A 105 18,413.27 x A = 0.0457 B versus D 83.04** x B x B = 0.0775 C versus D 12.38** x D x D = 0.0494 x 0 x 0 = 0.0496 NOTE.—p, number of parameters in the model; ln L, log-likelihood score; !,the dN/dSratio for the branches; ! A ,! B ,and! D are the dN/dSratios for branches lineages A, B, and D, respectively (see supplementary fig. S14,Supplementary Material online); ! 0 is the background dN/dSratio for the rest branch(es); 2ln L, twice the log-likelihood difference of the models compared. **Very significant (P<0.01). 2521 Ovine Mitogenomic Variations across Eastern Eurasia .doi:10.1093/molbev/msv139 MBE at Natural Resources Institute Finland (Luke) on October 6, 2016http://mbe.oxfordjournals.org/Downloaded from approaches, respectively (supplementary fig. S16, Supplementary Material online). Prehistoric Population Expansions Bayesian skyline plot (BSP) reconstructions of historical population expansions using thecompletemitogenomesrevealed the profile of predomestic change in Ne over large time scales. Based on the estimated TMRCA for the lineages (~0.79 Ma) from complete mitogenomes, the ovine lineages showed a steep increase in Ne at approximately 20–60 ka (supplementary fig. S17,Supplementary Material online). A prehistoric steep increase in Ne was also identified in the simulations of the partial Cyt-band control region sequences (supplementary fig. S18,Supplementary Material online). However, population growth was found to have occurred at approximately 50–300 ka, much earlier than the time obtained from simulations of the complete mitogenomes (~20– 60 ka; supplementary fig. S17,Supplementary Material online). Postdomestic Lineage Expansions A synthetic map constructed with the use of interpolated  1 values, the eigenvalues for the first multidimensional scaling (MDS) plot dimension, allows us to examine the gradients of colonization out of the sheep domestication center that peak intheMiddleEast(fig. 5A;supplementary table S12, Supplementary Material online).  1 explains 69.3% of the total variation. We observe a significant correlation between the  1 eigenvalues of Asian populations and their geographic distances from the domestication center (lineages A, B, and C of Central and East Asian populations: r= 0.201; P<0.05; lineage A of Arabian and Indian populations: r= 0.547; P<0.01; see fig. 5Cand D). This suggests that the major colonization process of the Middle Eastern sheep to eastern Eurasia (including Mongolia, China, and India) was through the Caucasus and Central Asia. The interpolation map of the  2 eigenvalues suggests that the second MDS dimension could represent genetic influence from the Mongolian Plateau region in China and the Indian subcontinent (fig. 5B).  2 explains 27.3% of the total variation. Its ranking shows the Mongolian Plateau region at one extreme, whereas the Indian subcontinent at the other extreme. This is supported by a strong and significant correlation observed between geographic distances from the putative region of initial colonization (i.e., the Mongolian Plateau region) and  2 values across eastern Eurasian populations (r= 0.372; P<0.01; fig. 5E). Eigenvalues  1 and  2 for all the populations are shown in supplementary table S12,Supplementary Material online. The star-like median-joining networks (supplementary fig. S2,Supplementary Material online) and mismatch distributions (supplementary fig. S19,Supplementary Material online) revealed genetic signatures of postdomestic demographic population expansions in lineages A, B, and C. The inference was corroborated by Fs(Fu 1997), Tajima’s D(1989), and scaled effective population size statistics (2Nu;Nrepresents the effective population size and udenotes the mutation rate). Both Fu’s Fs and Tajima’s Dstatistics showed Table 2. Divergence Time Estimated by the Sequences of Complete Mitogenomes and the Protein-Coding Genes (synonymous mutation and the third-codon position) Using ML and BI Methods. Method Data Set Model Node Node 1 (T C/E ) Ma Node 2 (T A/B ) Ma Node 3 (T AB/D ) Ma Node 4 (T ABD/CE ) Ma Node 5 (T O.aries/O.vignei ) Ma T O.aries/O.ammon Ma T O.aries/O.canadensis Ma ML Mitogenome Global Time 0.36 0.51 0.74 0.88 2.60 3.00 7.72 95%(CI) (0.278–0.439) (0.402–0.616) (0.613–0.867) (0.743–1.013) — (2.673–3.323) (6.567–8.883) Local Time 0.34 0.53 0.78 0.93 2.60 3.06 8.15 95%(CI) (0.276–0.472) (0.397–0.668) (0.600–0.956) (0.721–1.131) — (2.697–3.419) (6.635–9.664) Synonymous Global Time 0.31 0.52 0.68 0.79 2.60 2.92 8.36 95%(CI) (0.217–0.405) (0.373–0.661) (0.536–0.829) (0.637–0.934) — (2.535–3.312) (6.441–10.286) Local Time 0.31 0.52 0.69 0.80 2.60 2.93 8.31 95%(CI) (0.200–0.418) (0.346–0.694) (0.494–0.887) (0.583–1.018) — (2.453–3.413) (6.182–10.436) Third codon Global Time 0.29 0.50 0.64 0.73 2.60 2.81 7.47 95%(CI) (0.190–0.390) (0.361–0.639) (0.497–0.783) (0.579–0.881) — (2.418–3.202) (6.157–8.783) Local Time 0.29 0.50 0.64 0.73 2.60 2.81 7.47 95%(CI) (0.190–0.390) (0.359–0.636) (0.498–0.783) (0.581–0.882) — (2.414–3.200) (6.158–8.786) BI Mitogenome Relaxed-molecular clock Median 0.35 0.55 0.85 0.92 2.60 2.68 6.13 95%HPD (0.130–0.641) (0.266–0.913) (0.390–1.413) (0.464–1.498) — (2.462–3.031) (2.464–11.618) Synonymous Median 0.41 0.61 0.96 1.06 2.60 2.62 5.89 95%HPD (0.147–0.772) (0.291–1.013) (0.478–1.604) (0.541–1.716) — (2.458–3.338) (5.456–11.598) Third codon Median 0.36 0.57 0.87 0.94 2.60 2.62 6.49 95%HPD (0.142–0.656) (0.27–0.912) (0.421–1.428) (0.472–1.496) — (2.461–3.083) (2.478–12.647) NOTE.—“—,” not available. 2522 Lv et al. .doi:10.1093/molbev/msv139 MBE at Natural Resources Institute Finland (Luke) on October 6, 2016http://mbe.oxfordjournals.org/Downloaded from significant (P F <0.001; P D <0.001; table 3)departuresfrom neutrality in the three lineages. Additionally, the observed mismatch distributions of the lineages were fitted to the sudden population expansion models with very low values of the sum of squared deviation (SSD 0.005; table 3) statistic and Harpending’s Raggedness index (Harpending 1994; r H = 0.021–0.031; P R <0.5; table 3). Furthermore, the estimated preand postexpansion scaled effective population sizes (2Nu) indicated an increase in the effective population size for each of the lineages (A: 0.0–77.81; B: 0.04–16.63, and C: 0.00–15.77; table 3). The postdomestic expansion time expressed in twice the number of generations multiplied by the mutation rate (=2 ut) was found to be 6.443 ka (90% CI = 4.279–7.569 ka), 6.811 ka (90% CI: 3.502–9.706 ka), and 4.549 ka (90% CI = 2.402–6.652 ka) for lineages A, B, and C, respectively, when assuming an initial expansion (i.e., lineage A involving European sheep; Tapio et al. 2006) to equal 9 ka (table 3 and fig. 6). Separate analyses for the four major geographic areas (the Middle East, India, East Asia, and Europe) resulted in wider confidence intervals than those for the combined analysis and showed somewhat different estimates of  (table 3). In particular, the expansion time for lineage C in the Middle East (3.910 ka; 90% CI: 2.818–5.120 ka) was more recent than that in East Asia (4.967 ka; 90% CI: 1.893–7.842 ka), whereas relatively earlier expansions in the Middle East were inferred for lineages A and B (table 3). Additionally, we found much later expansions of lineages A and B in India (lineage A: 4.033 ka: 90% CI: 1.517–23.100 ka; lineage B: 3.393 ka; 90% CI: 0.961–16.311 ka) than those in East Asia (lineage A: 5.877 ka; 90% CI: 5.216–6.681 ka; lineage B: 7.008 ka; 90% CI: 3.030–16.348 ka), respectively. ABC analyses based on the control region sequences identified an optimal model for each of the five sets of candidate colonization models (lineage A first-step, lineage A secondstep, lineage B first-step, lineage B second-step, and lineage C; supplementary information S4,Supplementary Material online). The optimal models exhibited much higher posterior probability and nonoverlapped 95% CIs as compared with other candidate models (table 4). These optimal models indicated that 1) lineage A first colonized from the Middle East to the Mongolian Plateau region and the Indian subcontinent separately, and later from the Mongolian Plateau region to North China, and then to Southwest China (fig. 6); 2) lineage B first colonized from the Middle East to the Mongolian Plateau region, and then from the Mongolian Plateau region to North and Southwest China and the Indian subcontinent separately (fig. 6); and 3) Lineage C first colonized from the Middle East to the Mongolian Plateau region, and later from the Mongolian Plateau region to North China, and then to the Indian subcontinent (e.g., Nepal) (fig. 6). AA B -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 0 2000 4000 6000 8000 λ1 (dimension 1) Distance (km) -0.04 -0.02 0.00 0.02 0.04 0.06 0.08 0.10 0 2000 4000 6000 λ1 (dimension 1) Distance (km) -2.0 -1.0 0.0 1.0 2.0 3.0 4.0 5.0 6.0 0 500 1000 1500 2000 2500 3000 λ2 (dimension 2) C D E Distance (km) FIG.5. Synthetic maps illustrating geographic variation of eigenvalues ()forthefirsttwoMDSdimensions( 1 and  2 ) and regression of versus geographic distance from the putative original site of colonization process. (A)syntheticmapfor 1 ,(B)syntheticmapfor 2 ,(C) regression of  1 versus geographic distances from the domestication center of sheep (represented by the geographic distance from the Kilis province of Turkey, where ancient domestic sheep are located; Demirci et al. 2013) for Asian populations (r= 0.201; P<0.05); (D) regression of  1 (based on lineage A only) versus geographic distances from the domestication center of sheep (represented by the geographic distance from the Kilis province of Turkey, where ancient domestic sheep are located; Demirci et al. 2013) for sheep populations from the Indian subcontinent (r= 0.547; P<0.01); and (E)regressionof 2 versus geographic distances from a putative “transportation hub” of the Mongolian Plateau region (represented by the geographic distance from the northernmost population [Transbaikal Finewool] sampled) for eastern Eurasian (including China, Mongolia, and India) populations (r= 0.372; P<0.01). 2523 Ovine Mitogenomic Variations across Eastern Eurasia .doi:10.1093/molbev/msv139 MBE at Natural Resources Institute Finland (Luke) on October 6, 2016http://mbe.oxfordjournals.org/Downloaded from the fit of the models to the data by comparing twice the loglikelihood difference (2ln L)toa 2 distribution with degrees of freedom equal to the difference in the number of parameters between the two models (Yang 1998). Divergence Time Estimation Phylogenetic relationships within the genus Ovis (see Results) inferred in previous analyses were used to estimate the divergence times between the major O. aries lineages using PAML v.4.7 (Yang 2007) and BEAST v.1.7.5 (Drummond and Rambaut 2007). Due to the lack of an exact fossil record between Ovis speciesforthecalibration,weusedacomprehensiveevolutionary framework (see Nomura et al. 2013;Jiang et al. 2014) to estimate the divergence time between O. aries and O. vignei. A phylogenetic tree including 24 species (supplementary table S16, Supplementary Material online) was inferred based on the 13 mtDNA protein-coding genes using the GTR + I + G model in MrBayes v.3.2.2 (Ronquist et al. 2011).Thedivergencetimeswere estimated based on five fossil calibration points (18.3–28.5 Ma between Bovinae and Caprinae, 52–58 Ma between Cetacea and hippopotamus, 434.1 Ma between baleen and toothed whales, 42.8–63.8 Ma between Caniformia and Feliformia, and 62.3–71.2 Ma between Carnivora and Perissodactyla; see Nomura et al. 2013;Jiang et al. 2014). We applied the obtained O. aries/O. vignei divergence time (2.6 Ma; see Results) and three models to estimate the divergence times between the five O. aries mtDNA lineages. Global and local clock models were implemented using the ML in PAML v.4.7 (Yang 2007) and the uncorrelated relaxed-clock model was implemented using BEAST v.1.7.5 (Drummond and Rambaut 2007). We only considered the 61 complete mitogenomes of wild sheep species and native breeds of domestic sheep (supplementary table S1;Supplementary Material online) and applied three strategies in the ML analysis: One considered the complete mitogenomes under the TN93 model, the second considered the synonymous mutations under the HKY85 model, and the third considered only the thirdcodon positions under the HKY85 model. Similarly, Bayesian Markov chain Monte Carlo (MCMC) analysis of molecular sequences was performed by applying the three strategies in the program BEAST v.1.7.5 (Drummond and Rambaut 2007). Parameters of prior distributions, including models of nucleotide substitution and the divergence time between O. aries and O. vignei, were set the same as in the ML analyses described above. Three independent runs were performed with 50 million iterations. Samples were drawn every 5,000 MCMC steps, with the first 25% samples discarded as burn-in. The results of the three independent runs were combined using the LogCombiner program (available at http:// beast.bio.ed.ac.uk/LogCombiner, last accessed October 16, 2014) from BEAST v.1.7.5 (Drummond and Rambaut 2007). Convergence was confirmed by effective sampling size (ESS) greater than 200 using the program Tracer v.1.5 (Drummond and Rambaut 2007; available at http://beast.bio.ed.ac.uk/ Tracer, last accessed December 26, 2014). BI of Population Expansions Based on the divergence time of internal nodes estimated above and the 51 complete mitogenomes of native domestic sheep breeds, we reconstructed the change in N e of O. aries through time using BSPs (Drummond et al. 2005). The analyses were also performed on the partial Cyt-band control region sequences. We ran three independent chains in each analysis using BEAST v.1.7.5 (Drummond and Rambaut 2007), with 50 million generations (after discarding the first 10% of sampled generations as burn-in) and samples drawn every 5,000 steps. We applied the HKY85 and TN93 models under relaxed-clock model for complete mitogenomes and partial Cyt-bsequences, respectively. In the analysis of control region sequences, we set similar parameter values of 200 million generations (after discarding the first 10% of sampled generations as burn-in) with samples drawn every 2,000 steps under the HKY85 and relaxed-clock models. The combination of three independent results and checks of convergencewereperformedfollowingthesameproceduresas described above. Signatures of population expansions were examined using Arlequin v.3.5 (Excoffier and Lischer 2010). First, the observed and expected mismatch distributions of pairwise differences between haplotypes were compared using Tajima’s D(Tajima 1989)andFu’sFs(Fu 1997) tests of neutrality. Furthermore, we estimated the parameters for the sudden population expansion model (Rogers 1995), and the fit of the data to the sudden population expansion model was tested. Harpending’s raggedness index (r H ;Harpending 1994)ofthe observed mismatch distribution was also calculated. Pvalues of the SSDs test to evaluate the fit and significance for the parameters (r H ,D, and Fs) were determined with 1,000 coalescent simulations using Arlequin v.3.5 (Excoffier and Lischer 2010). All the partial control region sequences were included in the calculations. Further, to corroborate our inference that the Mongolian Plateau region serves as a “transportation hub” in eastern Eurasia (see Results), we distinguished several candidate colonization scenarios for the three main O. aries lineages (A, B, and C) using the ABC (Beaumont et al. 2002) procedure in DIYABC v.2.0.4 (Cornuetetal.2014). By incorporating all the O. aries control region sequences of a 292-bp-long hypervariable fragment, we tested six, eight, and five colonization models regarding potential migration routes from the Middle Eastern domestication center to different regions in eastern Eurasia (e.g., the Mongolian Plateau region, North China, Southwest China, and the Indian subcontinent) for the lineages A, B, and C, respectively (supplementary fig. S25,Supplementary Material online). Detailed information about the ABC analyses including the candidate colonization models tested was provided in supplementary information S4,Supplementary Material online. Supplementary Material Supplementary information S1–S4,figures S1–S25,andtables S1–S19 areavailableatMolecular Biology and Evolution online (http://www.mbe.oxfordjournals.org/). 2530 Lv et al. .doi:10.1093/molbev/msv139 MBE at Natural Resources Institute Finland (Luke) on October 6, 2016http://mbe.oxfordjournals.org/Downloaded from Acknowledgments The authors thank San-Gang He, Ya-Wei Sun, Nurbi Marzanov, Mikhail Ozerov, Maciek Murawski, Tatiana Kiseleva,andthelateMirjanaCinkulovforhelpinsample collection, Anneli Virta for technical assistance, and Dr Alessandro Achilli (Universit adiPerugia,Perugia,Italy)for his comments on an earlier version of the manuscript. This work was supported by the 100-talent Program of Chinese Academy of Sciences (CAS), the National High Technology Research and Development Program of China (i.e., 863 Program, grant No. 2013AA102506), the Breakthrough Project of Strategic Priority Program of the Chinese Academy of Sciences (grant No. XDB13000000), the grants from National Natural Science Foundation of China (grants Nos. 31272413 and U1303284), and Academy of Finland (grants Nos. 250633 and 256077) as well as Chinese Government contribution to CAAS-ILRI Joint Laboratory on Livestock and Forage Genetic Resources in Beijing. The paper contributes to the CGIAR Research Program on Livestock and Fish. References Achilli A, Bonfiglio S, Olivieri A, Malus a A, Pala M, Kashani BH, Perego UA, Ajmone-Marsan P, Liotta L, Semino O, et al. 2009. The multifaceted origin of taurine cattle reflected by the mitochondrial genome. PLoS One 4:e5753. Achilli A, Olivieri A, Soares P, Lancioni H, Kashani BH, Perego UA, Nergadze SG, Carossa V, Santagostino M, Capomaccio S, et al. 2012. Mitochondrial genomes from modern horses reveal the major haplogroups that underwent domestication. Proc Natl Acad Sci U S A. 109:2449–2454. Altschul SF, Madden TL, Sch€ affer AA, Zhang J, Zhang Z, Miller W, Lipman DJ. 1997. Gapped BLAST and PSI-BLAST: a new generation of protein database search programs. Nucleic Acids Res. 25:3389–3402. Aniwashi J, Jiahan K, Hakaimofu H, Sulaiman Y, Xi S-Y, Du M, Hailati, Tuersenhali, Ayinuer. 2007. Study on hybridication of wild argali and Bashibai sheep. Xinjiang Agric Sci. 44:702–705 (in Chinese). Bandelt HJ, Forster P, R€ ohl A. 1999. Median-joining networks for inferring intraspecific phylogenies. MolBiolEvol.16:37–48. Bar-Yosef O, Meadow R. 1995. The origins of agriculture in the Near East. In:PriceT,GebauerA-B,editors.Lasthunters,firstfarmers.SantaFe: School of American Research Press. p. 39–94. Beaumont MA, Zhang W, Balding DJ. 2002. Approximate Bayesian computation in population genetics. Genetics 162:2025–2035. Bj€ ornerfeldt S, Webster MT, Vil a C. 2006. Relaxation of selective constraint on dog mitochondrial DNA following domestication. Genome Res. 16:990–994. Bruford MW. 2005. Molecular approaches to understanding animal domestication: what have we learned so far? World Poultry Science Association, 4th European Poultry Genetics Symposium. Dubrovnik, Croatia, 6–8 October, 2005 Cai DW, Han L, Zhang XL, Zhou H, Zhu H. 2007. DNA analysis of archaeological sheep remains from China. JArchaeolSci. 34:1347–1355. CaiDW,TangZW,YuHX,HanL,RenXY,ZhaoXB,ZhuH,ZhouH. 2011. Early history of Chinese domestic sheep indicated by ancient DNA analysis of Bronze Age individuals. JArchaeolSci.38:896–902. Carruthers D. 1949. Beyond the Caspian. A naturalist in Central Asia. Edinburgh: Oliver and Boyd. Chen FH, Dong GH, Zhang DJ, Liu XY, Jia X, An CB, Ma MM, Xie YW, Barton L, Ren XY, et al. 2015. Agriculture facilitated permanent human occupation of the Tibetan Plateau after 3600 B.P. Science 347:248–250. Chen SY, Duan ZY, Sha T, Xiangyu J, Wu SF, Zhang YP. 2006. Origin, genetic diversity, and population structure of Chinese domestic sheep. Gene 376:216–223. Chen WH. 1990. Index to data of agricultural archaeology-farm-tools. Agric Archaeol. 1:425–427. (in Chinese) ChessaB,PereiraF,ArnaudF,AmorimA,GoyacheF,MainlandI,Kao RR, Pemberton JM, Beraldi D, Stear MJ, et al. 2009. Revealing the history of sheep domestication using retrovirus integrations. Science 324:532–536. Colledge S, Conolly J, Shennan S. 2005. The evolution of Neolithic farming from SW Asian origins to NW European limits. Eur J Archaeol. 8:137–156. Cornuet J-M, Pudlo P, Veyssier J, Dehne-Garcia A, Gautier M, LebloisR,MarinJ-M,EstoupA.2014.DIYABCv2.0:asoftware to make approximate Bayesian computation inferences about population history using single nucleotide polymorphism, DNA sequence and microsatellite data. Bioinformatics 30:1187–1189. Crandall KA, Kelsey CR, Imamichi H, Lane HC, Salzman NP. 1999. Parallel evolution of drug resistance in HIV: failure of nonsynonymous/synonymous substitution rate ratio to detect selection. Mol Biol Evol. 16:372–382. Darriba D, Taboada GL, Doallo R, Posada D. 2012. jModelTest 2: more models, new heuristics and parallel computing. Nat Methods. 9:772. Demirci S, Koban Bas¸tanlar E, Da gtas¸ND,Pis¸kin E, Engin A, € Ozer F, Y€ unc€ uE,Do  gan S¸A, Togan _ I. 2013. Mitochondrial DNA diversity of modern, ancient and wild sheep (Ovis gmelinii anatolica)from Turkey: new insights on the evolutionary history of sheep. PLoS One 8:e81952. Dobney K, Larson G. 2006. Genetics and animal domestication: new windows on an elusive process. J Zool. 269:261–271. Drummond AJ, Rambaut A. 2007. BEAST: Bayesian evolutionary analysis by sampling trees. BMC Evol Biol. 7:214. Drummond AJ, Rambaut A, Shapiro B, Pybus OG. 2005. Bayesian coalescent inference of past population dynamics from molecular sequences. MolBiolEvol.22:1185–1192. Excoffier L, Lischer HEL. 2010. Arlequin suite ver 3.5: a new series of programs to perform population genetics analyses under Linux and Windows. MolEcolResour.10:564–567. Felsenstein J. 1985. Confidence limits on phylogenies: an approach using the bootstrap. Evolution 39:783–791. Fu Y-X. 1997. Statistical tests of neutrality of mutations against population growth, hitchhiking and background selection. Genetics 147:915–925. Gauri FN. 2013. Indo-Saudi trade relation. Arabian Journal of Business and Management Review (Nigerian Chapter) 1(2):45–57. Goldman N, Yang Z. 1994. A codon-based model of nucleotide substitution for protein-coding DNA sequences. MolBiolEvol. 11:725–736. Grant WS, Liu M, Gao T, Yanagimoto T. 2012. Limits of Bayesian skyline plot analysis of mtDNA sequences to infer historical demographies in Pacific herring (and other species). Mol Phylogenet Evol. 65:203–212. Guindon S, Gascuel O. 2003. A simple, fast, and accurate algorithm to estimate large phylogenies by maximum likelihood. Syst Biol. 52:696–704. Guo J, Du LX, Ma YH, Guan WJ, Li HB, Zhao QJ, Li X, Rao SQ. 2005. A novel maternal lineage revealed in sheep (Ovis aries). Anim Genet. 36:331–336. Hanotte O, Bradley DG, Ochieng JW, Verjee Y, Hill EW, Rege JEO. 2002. African pastoralism: genetic imprints of origins and migrations. Science 296:336–339. Harpending H. 1994. Signature of ancient population growth in a lowresolution mitochondrial DNA mismatch distribution. Hum Biol. 66:591–600. Hasegawa M, Kishino H, Yano T-A. 1985. Dating of the human-ape splitting by a molecular clock of mitochondrial DNA. J Mol Evol. 22:160–174. 2531 Ovine Mitogenomic Variations across Eastern Eurasia .doi:10.1093/molbev/msv139 MBE at Natural Resources Institute Finland (Luke) on October 6, 2016http://mbe.oxfordjournals.org/Downloaded from Heller R, Chikhi L, Siegismund HR. 2013. The confounding effect of population structure on Bayesian skyline plot inferences of demographic history. PLoS One 8:e62992. Hiendleder S, Lewalski H, Wassmuth R, Janke A. 1998. The complete mitochondrial DNA sequence of the domestic sheep (Ovis aries) and comparison with the other major ovine haplotype. J Mol Evol. 47:441–448. Hiendleder S, Mainz K, Plante Y, Lewalski H. 1998. Analysis of mitochondrial DNA indicates that domestic sheep are derived from two different ancestral maternal sources: no evidence for contributions from urial and argali sheep. J Hered. 89:113–120. Ho SYW. 2014. The changing face of the molecular evolutionary clock. Trends Ecol Evol. 29:496–503. Ho SYW, Duch^ ene S. 2014. Molecular-clock methods for estimating evolutionary rates and timescales. Mol Ecol. 23:5947–5965. Ho SYW, Lanfear R, Bromham L, Phillips MJ, Soubrier J, Rodrigo AG, Cooper A. 2011. Time-dependent rates of molecular evolution. Mol Ecol. 20:3087–3101. Ho SYW, Shapiro B. 2011. Skyline-plot methods for estimating demographic history from nucleotide sequences. Mol Ecol Resour. 11:423–434. Hodges J. 1999. Animals and values in society. Livest Res Rural Dev. 11:1–6. Hunt HV, Vander Linden M, Liu X, Motuzaite-Matuzeviciute G, Colledge S, Jones MK. 2008. Millets across Eurasia: chronology and context of early records of the genera Panicum and Setaria from archaeological sites in the Old World. Veg Hist Archaeobot. 17:5–18. Jiang Y, Xie M, Chen W, Talbot R, Maddox JF, Faraut T, Wu C, Muzny DM, Li Y, Zhang W, et al. 2014. The sheep genome illuminates biology of the rumen and lipid metabolism. Science 344:1168–1173. Jones MK, Liu X. 2009. Origins of agriculture in East Asia. Science 324:730–731. Joost S, Bonin A, Bruford MW, Despres L, Conord C, Erhardt G, Taberlet P. 2007. A spatial analysis method (SAM) to detect candidate loci for selection: towards a landscape genomics approach to adaptation. Mol Ecol. 16:3955–3969. Kijas JW, Lenstra JA, Hayes B, Boitard S, Neto LRP, San Cristobal M, Servin B, McCulloch R, Whan V, Gietzen K, et al. 2012. Genomewide analysis of the world’s sheep breeds reveals high levels of historic mixture and strong recent selection. PLoS Biol. 10:e1001258. Kijas JW, Townley D, Dalrymple BP, Heaton MP, Maddox JF, McGrath A, Wilson P, Ingersoll RG, McCulloch R, McWilliam S, et al. 2009. A genome wide survey of SNP variation reveals the genetic structure of sheep breeds. PLoS One 4:e4668. Kumar S. 2005. Molecular clocks: four decades of evolution. Nat Rev Genet. 6:654–662. Kuo F, Needham J, Chh^ eng CT. 1999. The history of zoology in China. Beijing (China): Science Press. (in Chinese) Lancioni H, Di Lorenzo P, Ceccobelli S, Perego UA, Miglio A, Landi V, Antognoni MT, Sarti FM, Lasagna E, Achilli A. 2013. Phylogenetic relationships of three Italian Merino-derived sheep breeds evaluated through a complete mitogenome analysis. PLoS One 8:e73712. Larkin MA, Blackshields G, Brown NP, Chenna R, McGettigan PA, McWilliam H, Valentin F, Wallace IM, Wilm A, Lopez R, et al. 2007. Clustal W and Clustal X version 2.0. Bioinformatics 23:2947–2948. Larson G, Albarella U, Dobney K, Rowley-Conwy P, Schibler J, Tresset A, Vigne J-D, Edwards CJ, Schlumbaum A, Dinu A, et al. 2007. Ancient DNA, pig domestication, and the spread of the Neolithic into Europe. Proc Natl Acad Sci U S A. 104:15276–15281. Larson G, Burger J. 2013. A population genetics view of animal domestication. Trends Genet. 29:197–205. Larson G, Liu R, Zhao X, Yuan J, Fuller D, Barton L, Dobney K, Fan Q, Gu Z, Liu X-H, et al. 2010. Patterns of East Asian pig domestication, migration, and turnover revealed by modern and ancient DNA. Proc Natl Acad Sci U S A. 107:7686–7691. Lippold S, Matzke NJ, Reissmann M, Hofreiter M. 2011. Whole mitochondrial genome sequencing of domestic horses reveals incorporation of extensive wild horse diversity during domestication. BMC Evol Biol. 11:328. Lu H, Zhang J, Liu K-B, Wu N, Li Y, Zhou K, Ye M, Zhang T, Zhang H, Yang X, et al. 2009. Earliest domestication of common millet (Panicum miliaceum) in East Asia extended to 10,000 years ago. Proc Natl Acad Sci U S A. 106:7367–7372. Luikart G, Gielly L, Excoffier L, Vigne J-D, Bouvet J, Taberlet P. 2001. Multiple maternal origins and weak phylogeographic structure in domestic goats. Proc Natl Acad Sci U S A. 98:5927–5932. Luo Y-Z, Cheng S-R, Lkhagva B, Badamdorj D, Hanotte O, Han J-L. 2005. Study on origin and genetic diversity of Mongolian and Chinese sheep. Acta Genet Sin. 32:1256–1265. (in Chinese) Lv F-H, Agha S, Kantanen J, Colli L, Stucki S, Kijas JW, Joost S, Li M-H, Marsan PA. 2014. Adaptations to climate-mediated selective pressuresinsheep.Mol Biol Evol. 31:3324–3343. Meadows JR, Cemal I, Karaca O, Gootwine E, Kijas JW. 2007. Five ovine mitochondrial lineages identified from sheep breeds of the Near East. Genetics 175:1371–1379. Meadows JRS, Hiendleder S, Kijas JW. 2011. Haplogroup relationships between domestic and wild sheep resolved using a mitogenome panel. Heredity 106:700–706. Muigai AT, Hanotte O. 2013. The origin of African sheep: archaeological and genetic perspectives. Afr Archaeol Rev. 30:39–50. Naderi S, Rezaei H-R, Pompanon F, Blum MGB, Negrini R, Naghash H-R, Balkız € O, Mashkour M, Gaggiotti OE, Ajmone-Marsan P, et al. 2008. The goat domestication process inferred from large-scale mitochondrial DNA analysis of wild and domestic individuals. Proc Natl Acad Sci U S A. 105:17659–17664. Niemi M, Bl€ auer A, Iso-Touru T, Nystr€ om V, Harjula J, Taavitsainen J-P, Stora ˚J, Lid en K, Kantanen J. 2013. Mitochondrial DNA and Y-chromosomal diversity in ancient populations of domestic sheep (Ovis aries) in Finland: comparison with contemporary sheep breeds. Genet Sel Evol. 45:2. Nomura K, Yonezawa T, Mano S, Kawakami S, Shedlock AM, Hasegawa M, Amano T. 2013. Domestication process of the goat revealed by an analysis of the nearly complete mitochondrial protein-encoding genes. PLoS One 8:e67775. Pe H. 1959. G. H. Luce and Pe Maung Tin (ed.): Inscriptions of Burma. Portfolio IV. Down to 702 B.E. (1340 A.D.).—Portfolio V. 703–726 B.E. (1341–1364 A.D.). (University of Rangoon Oriental Studies Publication No. 5, No. 6.) 63 pp., plates 346–462; 38 pp., plates 463–609. Oxford: University Press, 1956. Bull Sch Orient Afr Stud. 22:177. Pedrosa S, Uzun M, Arranz JJ, Gutierrez-Gil B, San Primitivo F, Bayon Y. 2005. Evidence of three maternal lineages in Near Eastern sheep supporting multiple domestication events. Proc R Soc Lond B Biol Sci. 272:2211–2217. Peng M-S, Zhang Y-P. 2011. Inferring the population expansions in peoplingofJapan.PLoS One 6:e21509. PereiraL,SilvaNM,Franco-DuarteR,FernandesV,PereiraJB,CostaMD, Martins H, Soares P, Behar DM, Richards MB, et al. 2010. Population expansion in the North African Late Pleistocene signalled by mitochondrial DNA haplogroup U6. BMC Evol Biol. 10:390. Peters J, von den Driesch A, Helmer D. 2005. The upper Euphrates-Tigris basin: cradle of agro-pastoralism. In: Vigne JD, Peters J, Helmer D, editors. The first steps of animal domestication. Oxford: Oxbow. pp. 96–124. Pickrell JK, Pritchard JK. 2012. Inference of population splits and mixtures from genome-wide allele frequency data. PLoS Genet. 8:e1002967. Poplin F. 1979. Origine du mouflon de Corse dans une nouvelle perspective pal eontologique: par marronnage. Ann G en et S el Anim. 11:133–134. Posada D. 2008. jModelTest: phylogenetic model averaging. Mol Biol Evol. 25:1253–1256. R Core Team. 2014. R: a language and environment for statistical computing. Vienna (Austria): R Foundation for Statistical Computing. Reynolds J, Weir BS, Cockerham CC. 1983. Estimation of the co-ancestry coefficient—basis for a short-term genetic-distance. Genetics 105:767–779. 2532 Lv et al. .doi:10.1093/molbev/msv139 MBE at Natural Resources Institute Finland (Luke) on October 6, 2016http://mbe.oxfordjournals.org/Downloaded from Rezaei HR, Naderi S, Chintauan-Marquier IC, Taberlet P, Virk AT, Naghash HR, Rioux D, Kaboli M, Pompanon F. 2010. Evolution and taxonomy of the wild species of the genus Ovis (Mammalia, Artiodactyla, Bovidae). Mol Phylogenet Evol. 54:315–326. Rogers AR. 1995. Genetic evidence for a Pleistocene population explosion. Evolution 49:608–615. Ronquist F, Huelsenbeck J, Teslenko M. 2011. Draft MrBayes version 3.2 Manual: tutorials and model summaries. Available from: http://mrbayes.sourceforge.net/. Rozas J, S anchez-DelBarrio JC, Messeguer X, Rozas R. 2003. DnaSP, DNA polymorphism analyses by the coalescent and other methods. Bioinformatics 19:2496–2497. Ryder ML. 1983. Sheep and man. London: Duckworth. Ryder ML. 1984. Sheep. In: Mason IL, editor. Evolution of domesticated animals. London/New York: Longman. pp. 63–85. Sambrook J, Russell DW. 2001. Molecular cloning: a laboratory manual. 3rd ed. New York: Cold Spring Harbor Laboratory Press. Scherf BD, editor. 2000. World watch list for domestic animal diversity. In: Food and Agriculture Organization of the United Nations. Rome: Food and Agriculture Organization of the United Nations. p. 58. Shaha R. 1970. Nepal, Tibet and China. JNepalCouncilWorldAff. 3:13–82. Shi N-N, Fan L, Yao Y-G, Peng M-S, Zhang Y-P. 2014. Mitochondrial genomes of domestic animals need scrutiny. Mol Ecol. 23:5393–5397. Singh S, Kumar S Jr, Kolte AP, Kumar S. 2013. Extensive variation and sub-structuring in lineage A mtDNA in Indian sheep: genetic evidence for domestication of sheep in India. PLoS One 8:e77858. Taberlet P, Valentini A, Rezaei HR, Naderi S, Pompanon F, Negrini R, Ajmone-Marsan P. 2008. Are cattle, sheep, and goats endangered species? Mol Ecol 17:275–284. Tajima F. 1989. Statistical method for testing the neutral mutation hypothesis by DNA polymorphism. Genetics 123:585–595. Tamura K, Nei M. 1993. Estimation of the number of nucleotide substitutions in the control region of mitochondrial DNA in humans and chimpanzees. MolBiolEvol.10:512–526. Tamura K, Peterson D, Peterson N, Stecher G, Nei M, Kumar S. 2011. MEGA5: molecular evolutionary genetics analysis using maximum likelihood, evolutionary distance, and maximum parsimony methods. MolBiolEvol.28:2731–2739. TapioM,MarzanovN,OzerovM,CinkulovM,GonzarenkoG,Kiselyova T, Murawski M, Viinalass H, Kantanen J. 2006. Sheep mitochondrial DNA variation in European, Caucasian, and Central Asian areas. Mol Biol Evol. 23:1776–1783. TapioM,OzerovM,TapioI,ToroMA,MarzanovN,CinkulovM, Goncharenko G, Kiselyova T, Murawski M, Kantanen J. 2010. Microsatellite-based genetic diversity and population structure of domestic sheep in northern Eurasia. BMC Genet. 11:76. Torroni A, Achilli A, Macaulay V, Richards M, Bandelt H-J. 2006. Harvesting the fruit of the human mtDNA tree. Trends Genet. 22:339–345. Wang G-D, Xie H-B, Peng M-S, Irwin D, Zhang Y-P. 2014. Domestication genomics: evidence from animals. Annu Rev Anim Biosci. 2:65–84. Wang X, Ma Y-H, Chen H. 2006. Analysis of the genetic diversity and the phylogenetic evolution of Chinese sheep based on Cyt b gene sequences. Acta Genet Sin. 33:1081–1086. WangZ,YonezawaT,LiuB,MaT,ShenX,SuJ,GuoS,HasegawaM,Liu J. 2011. Domestication relaxed selective constraints on the yak mitochodrial genomes. Mol Biol Evol. 28:1553–1556. WarmuthV,ErikssonA,BowerMA,BarkerG,BarrettE,HanksBK,LiS, Lomitashvili D, Ochir-Goryaeva M, Sizonov GV, et al. 2012. Reconstructing the origin and spread of horse domestication in the Eurasian steppe. Proc Natl Acad Sci U S A. 109:8202–8206. Wood NJ, Phua SH. 1996. Variation in the control region sequence of the sheep mitochondrial genome. Anim Genet. 27:25–33. WuGS,YaoYG,QuKX,DingZL,LiH,PalanichamyMG,DuanZY,LiN, Chen YS, Zhang YP. 2007. Population phylogenomic analysis of mitochondrial DNA in wild boars and domestic pigs revealed multiple domestication events in East Asia. Genome Biol. 8:R245. Xiang H, Gao J, Yu B, Zhou H, Cai D, Zhang Y, Chen X, Wang X, Hofreiter M, Zhao X. 2014. Early Holocene chicken domestication in northern China. Proc Natl Acad Sci U S A. 111:17564–17569. Yang X, Scuderi LA, Wang X, Scuderi LJ, Zhang D, Li H, Forman S, Xu Q, Wang R, Huang W, et al. 2015. Groundwater sapping as the cause of irreversible desertification of Hunshandake Sandy Lands, Inner Mongolia, northern China. Proc Natl Acad Sci U S A. 112:702–706. Yang X, Wan Z, Perry L, Lu H, Wang Q, Zhao C, Li J, Xie F, Yu J, Cui T, et al. 2012. Early millet use in northern China. Proc Natl Acad Sci U S A. 109:3726–3730. Yang Z. 1998. Likelihood ratio tests for detecting positive selection and application to primate lysozyme evolution. Mol Biol Evol. 15:568–573. Yang Z. 2007. PAML 4: phylogenetic analysis by maximum likelihood. MolBiolEvol.24:1586–1591. Yao Y-G, Salas A, Logan I, Bandelt H-J. 2009. mtDNA data mining in GenBank needs surveying. Am J Hum Genet. 85:929–933. Yi H. 2004. Bronze roads: a introduction to archaic cultural exchange in Eurasia. In: Department of Cultural Heritage and Museum Studies, editor. Antiquities of Eastern Asia. Beijing (China): Cultural Relics Press. (in Chinese) Zeder MA. 2008. Domestication and early agriculture in the Mediterranean Basin: origins, diffusion, and impact. Proc Natl Acad Sci U S A. 105:11597–11604. Zeder MA, Emshwiller E, Smith BD, Bradley DG. 2006. Documenting domestication: the intersection of genetics and archaeology. Trends Genet. 22:139–155. Zhang H, Paijmans JLA, Chang F, Wu X, Chen G, Lei C, Yang X, Wei Z, Bradley DG, Orlando L, et al. 2013. Morphological and genetic evidence for early Holocene cattle management in northeastern China. Nat Commun. 4:2755. 2533 Ovine Mitogenomic Variations across Eastern Eurasia .doi:10.1093/molbev/msv139 MBE at Natural Resources Institute Finland (Luke) on October 6, 2016http://mbe.oxfordjournals.org/Downloaded from