Refinement of Global Domestic Horse Biogeography Using Historic

31590_7esz032.pdf 769

Journal of Heredity

, 2019, 769-781 doi:10.1093/jhered/esz032

Original Article

Advance Access publication October 19, 2019

© The American Genetic Association 2019. All rights reserved. For per missions, please e-mail: journals.permissions@oup.com

Original Article

Re nement of Global Domestic Horse

Biogeography Using Historic Landrace Chinese

Mongolian Populations

Haige Han, Kenneth Bryan, Wunierfu Shiraigol, Dongyi Bai, Yiping Zhao, Wuyingga Bao, Siqin Yang, Wengang Zhang, DavidE. MacHugh,

Manglai Dugarjaviin, and EmmelineW. Hill

From the UCD School of Agriculture and Food Science, University College

Dublin, Beleld, Dublin D04 V1W8, Ireland

(Han, Bryan, MacHugh, and Hill); College of Animal Science, Inner Mong olia Agricultural University, Hohhot, P.R.

China (Shiraigol, Bai, Zhao, and Dugarjaviin); Vocational and Technical College of Inner Mongolia Agricultural

University, Hohhot, P.R. China (Bao); Inner Mongolia University for the Nationalities, Tongliao, P.R. China (Yang);

Institute of Animal Science, Chinese Academy of Agriculture Sciences, Be ijing, P.R. China (Zhang); and UCD Conway Institute of Biomolecular and Biomedical Research, University College Du blin, Beleld, Dublin D04 V1W8, Ireland (MacHugh). Address correspondence to E.W. Hill at the address above, or e-mail: emmeline.hill@ucd.ie . Address correspondence to M.Dugarjaviin at the address above, or e-m ail: dmanglai@163.com . Received September 27, 2018; First decision November 12, 2018; Accepted

August 1, 2019.

Corresponding Editor: Ernest Bailey

Abstract

The Mongolian horse is one of the oldest extant horse populations and although domesticated, most animals are free-ranging and experience minimal human intervention. As an ancient population originating in one of the key domestication centers, the Mongolian horse may play a key role in understanding the origins and recent evolutionary history of horses. Here we describe an analysis of high-density genome-wide single-nucleotide polymorphism (SNP) data in 40 globally dispersed horse populations ( n = 895). In particular, we have focused on new results from Chinese Mongolian horses ( n = 100) that represent 5 distinct populations. These animals were genotyped for 670K SNPs

and the data were analyzed in conjunction with 35K SNP data for 35 distinct breeds. Analyses of these

integrated SNP data sets demonstrated that the Chinese Mongolian populat ions were genetically distinct from other modern horse populations. In addition, compared to other domestic horse breeds, the Chinese Mongolian horse populations exhibited relatively high genomic diversity. These results suggest that, in genetic terms, extant Chinese Mongolian horses may be the most similar modern populations to the animals originally domesticated in this region of Asia. Chinese Mongolian horse populations may therefore retain ancestral genetic variants from the ear liest domesticates. Further genomic characterization of these populations in conjunction with archaeogenetic sequence data should be prioritized for understanding recent horse evolution and the domestication process that has led to the wealth of diversity observed in modern global horse breeds.

Subject areas:

Conservation genetics and biodiversity, Population structure and phylogeography

Keywords:

genetic variation, horse domestication, single-nucleotide polymorphism Downloaded from https://academic.oup.com/jhered/article/110/7/769/5601164 by guest on 28 August 2023 The domestication of the horse (Equus caballus) approximately

5500 years ago (

Larson and Fuller 2014

; Larson et al. 2014) had a major impact on subsequent human migrations and the cultural and economic history of Eurasia (

Anthony 2010

). Consequently, understanding the spatial and temporal origins of domestic horses is of general, as well as scienti?c interest. The earliest archaeological evidence indicates that horse domestication likely started on the Kazakh steppe with the Botai culture ~5500 years ago (

Outram et al.

2009
); although it is no longer evident that these domesticates were the ancestors of modern horses (

Gaunitz et al. 2018

). Archaeological and genetic evidence suggests that horse domestication was complex (

Librado et al. 2017

; Gaunitz et al. 2018) and likely involved mul- tiple horse domestication centers across Eurasia as well as the Iberian

Peninsula (

Cieslak et al. 2010

;

Lira et al. 2010

). During the 5 millennia since domestication, the expansion of horse populations has followed complex geographical routes, mainly related to the horse's use in warfare and for transportation. More recently, from the 18th century onwards, humans have reshaped horses through arti?cial selection and breeding for speci?c charac - teristics such as speed, athleticism, height, gait, and coat color ( Hill et al. 2010 ; Petersen et al. 2013b). Consequently, although distin- guishable from landrace populations, many of the so-called "modern breeds" have recent origins and have been developed through genetic improvement for human requirements, with the largest numbers of diverse horse breeds found principally in Europe and North America (

Petersen et al. 2013a

). Due to the extinction of wild horses as well as the rarity of arch - aeological evidence, genetic data have emerged as a valuable tool for studying the evolutionary history of domestic horses. However, to date the absence of genetic data spanning the entire temporal and geographical range of horse domestication has prevented the establishment of ?rm conclusions regarding the de?nitive biogeo - graphic origins of domestic horses. Previous studies of extant horse genetic diversity have largely relied on modern European and North American breeds as reference populations, which have relatively re - cent histories, long after the ?rst domestic horses appear in the arch - aeological record (

Petersen et al. 2013a

; Wallner et al. 2017). Very recently, an ancient DNA-based study has reported 2 other lineages at the far eastern (Siberia) and western range of Eurasian (Iberia) apart from 2 extant lineages (domestic and Przewalski's horses), and none of them contributed signi?cantly to the ancestry of modern do - mesticates (

Fages et al. 2019

). Therefore, the origins of the modern domestic horse remain unanswered. It is generally believed that the central Asian steppe region, including Mongolia, may be one of the centers for horse domesti - cation (

Larson and Fuller 2014

). Mongolian horse is recognized as one of oldest known horse breeds in the world and has been lived in Asia for thousands of year (

Hendricks 2007

). The Mongols are known as an "ethnic group on horseback," indicating the central role of horses for their nomadic life. For thousands of years, the life of Mongols has centered on the horse for transportation, a source of food (both for horsemeat and mare's milk), for use in warfare as well as cultural signi?cance. Present-day Mongolian horses are mainly distributed in parts of northeast and north China (mainly in Inner Mongolia), the Mongolia People's Republic, and some areas of eastern Russian (

Chang 2009

). For the present study, we have focused on the Mongolian horse and derived breeds from Inner Mongolia, China, hereafter referred to as the Chinese Mongolian horse. Chinese Mongolian horses are found in the Inner Mongolia Autonomous Region western parts of the provinces of Heilongjiang,

Jilin, and Liaoning and some parts of Xinjiang Uygur Autonomous Region in China (Hendricks 2007; National Livestock and Poultry

Genetic Resources Committee 2011

). The Chinese Mongolian horse is generally classi?ed as a breed in toto but as a result of long-ter m selection for adaptation to local environments by herdsmen, several phenotypically distinct subtypes have evolved: Wushen (desert type), Wuzhumuqin (steppe type), and Baicha Iron Hoof (mountain type). Other breeds (e.g., Sanhe and Abaga Black) found in Inner Mongolia have been developed from local Mongolian stock following admix - ture with imported horses from Europe and elsewhere. To date, Chinese Mongolian horses are mostly exploited as sources of meat and milk, serve as transportation for the daily work of the nomads, and are used for horse racing. Mongolian horses representing landrace populations from cen - tral Asia have large census population sizes and have not experi - enced intense selection like many modern breeds (

Hendricks 2007

). This demographic history would likely lead to retention of ancestral predomestic alleles at multiple loci across the genome. Therefore, we hypothesized that high-resolution population genomics and ?ne- scale characterization of the genetic structure of Mongolian horse populations would shed light on domestication events in this region. The recent advent of high-throughput and cost-effective genotyping technologies makes it possible to assess genetic diversity within and among global horse populations at a genome-wide scale. For the present study, we present an analysis of horse genetic di - versity using samples from 40 modern horse populations ( n = 895). We have focused speci?cally on new data generated for Chinese

Mongolian horses (

n = 100) from 5 distinct populations that were genotyped for 670K single-nucleotide polymorphisms (SNPs). These data were integrated with publicly available 50K SNP data from 35 global breeds (

Petersen et al. 2013a

). The primary aims of this study were to test the following hypoth - eses: 1) genetic variation is higher in Chinese Mongolian landrace populations than other modern global breeds; and 2) this genetic variation re?ects an ancient ancestry associated with domestication events in the Mongolian steppe region.

Methods

Chinese Mongolian Horse Populations, DNA

Sampling, and DNA Extraction

Samples were collected from

n = 155 horses representing 5 popu - lations in Inner Mongolia, China (Abaga Black [AB], Baicha Iron Hoof [BCIH], Sanhe [SH], Wushen [WS], Wuzhumuqin [WZ]) (

Supplementary Table S1

). WS, WZ, and BCIH represent 3 subtypes of Chinese Mongolian horse. AB is recognized as an indigenous breed distinct from other Mongolian subtypes and its breeding his - tory can be traced back to the 13th century. On the other hand, the SH population has only been recognized as an established breed since 1986 (

National Livestock and Poultry Genetic Resources

Committee 2011

). Horses were sampled in a range of locations across the Inner Mongolia Autonomous Region, China and sampling focused on ani - mals bred locally (

Figure 1

). To obtain a representative sample of genetic variation in these populations, hair samples were collected in 22 different villages. Since there was no pedigree information available for the horses, sampling of related individuals was avoided using data from a questionnaire completed by horse owners that documented, where possible, known relationships. Genomic DNA was extracted using the Qiagen DNeasy Blood & Tissue Kit 250 (Qiagen GmbH, Hilden, Germany) following manufacturer's instructions. Genomic DNA samples were quanti?ed, 770

Journal of Heredity, 2019, Vol. 110, No. 7Downloaded from https://academic.oup.com/jhered/article/110/7/769/5601164 by guest on 28 August 2023

and quality checked using a NanoDrop 1000 Spectrophotometer (Thermo Fisher Scienti?c, Waltham, MA) and stored at -20 °C.

SNP Genotyping and Quality Control

The 155 Chinese Mongolian horses were genotyped with the Axiom Equine Genotyping Array (Axiom MNEC670; Affymetrix- Thermo Fisher Scienti?c). Each DNA sample was diluted to a con - centration of 400 ng in 40 L and the samples were genotyped in

2 separate batches (Batch 1, n = 93; Batch 2, n = 63), including

1 pair of duplicate samples for quality control (QC) purposes

(100% genotyping concordance was observed for this pair of samples). Since no pedigree information was available, genomic relationships among horses were estimated using autosomal iden - tity by descent (pi-hat) values in PLINK v1.9 (

Chang et al. 2015

). Individual animals with pi-hat values > 0.3 were removed (after pruning for minor allele frequency [MAF] < 0.05) and also where results of a genetic sex check (PLINK --sex-check) did not match the sex description in the questionnaire. Individual animal samples with SNP call rates < 90% were also excluded from the ?nal data set. For SNPs pruning, autosomal 611K SNPs were retained using the command chr 1-31 in PLINK. Then SNPs were also pruned for genotype rate (<0.99), MAF (<0.05), and Hardy-Weinberg equilibrium (<10 -5 ). Following QC, genotypes for n = 100 horses remained with a total of 358 824 autosomal SNPs that could be used for population genomics analyses of the Chinese Mongolian horses (without comparator populations). Assembly of Comparative SNP Data Sets, QC, and Filtering of SNPs For comparative population genomics analyses, data for animals ( n = 795) from 35 distinct populations with 50 042 autosomal SNPs (PLINK-formatted, Illumina top-stranded) were obtained ( www. animalgenome.org/repository/pub/UMN2013.0125 ). To integrate the new Affymetrix SNP data with the public Illumina data, the Affymetrix SNPs for the Chinese Mongolian horses were ?rst con - verted to Illumina top strand format using a custom python script. To eliminate ascertainment bias, as far as possible, we preprocessed the public Illumina 50K data using a method described by

Petersen

et al. (2013a) , such that horses from the SNP discovery breeds (Akhal Teke, Andalusian, Arabian, Icelandic, Quarter Horse, Standardbred, and Thoroughbred) were removed from the Illumina data set, which was then pruned to exclude SNPs with MAF < 0.05. All remaining horses were then included, and SNPs that did not pass QC were re - moved from all analyses. Finally, the converted Affymetrix SNPs and

Figure 1. Chinese Mongolian horses representing ve distinct populations and sampling locations in Inner Mongolia, China. Map modied from original by Uwe

Dedering (Creative Commons licence

https://creativecommons.org/licenses/by-sa/3.0/deed.en CC BY-SA 3.0). Satellite data from Google Maps (credits: Google LLC, SK telecom Co., Ltd., ZENRIN Co., Ltd.). See online version for full colours.

Journal of Heredity,

2019, Vol. 110, No. 7

771Downloaded from https://academic.oup.com/jhered/article/110/7/769/5601164 by guest on 28 August 2023

the preprocessed Illumina data were merged by using the bmerge command in PLINK. SNP QC and ?ltering were performed across populations on the merged data set: individual SNPs with >1% missing data and MAF < 0.05 were removed for evaluation of runs of homozygosity (ROH), linkage disequilibrium (LD) analysis, principal component analysis (PCA), construction of phylogenetic trees, and investigation of popu - lation structure. Consequently, a total of 35 414 SNPs were used for these downstream genetic analyses. An additional LD pruning step was applied for estimates of individual inbreeding coef?cients ( F), expected heterozygosity ( H e ), and pairwise F ST estimation, by pruning for r 2 < 0.1 in PLINK using 100 SNP windows and incrementally advancing 25 SNPs per set (--indep-pairwise 100 25 0.1). After LD pruning, a total of 14 012 SNPs were used for these genetic analyses. We evaluated the performance of the SNP integration pipeline by estimating the correlation coef?cients for the mean F and H e values in the present study, and that of

Petersen et al. (2013a)

, which were r 2 = 0.97 ( P < 0.001) and r 2 = 0.94 ( P < 0.001), respectively (

Supplementary Figure S1

). We obtained publicly available whole genome sequence data (BAM-formatted) from 17 Przewalski's horses ( https://www.ebi. ac.uk/ena/data/view/PRJEB10098 ). SNPs were called using the mpileup command in SAMtools (1.3.1) (

Li et al. 2009

). SNPs were then called using bcftools with the output from mpileup and subse - quently 35K common SNPs were extracted using the positions com - mand from VCFtools (

Danecek et al. 2011

).

Within-Breed Diversity

Expected heterozygosity (

H e ) and individual inbreeding coef?cients ( F ) were computed using PLINK. H e within each population was esti - mated using the Hardy-Weinberg test statistics command (--hardy). F was estimated based on the observed versus expected number of homozygote counts for each sample (--het). Mean F and H e estimates were determined for each population. Genome-wide LD was esti - mated by calculating squared correlations ( r 2 ) between all pairs of SNPs for a particular population (a minimum of 15 horses per popu - lation). Pairwise LD was estimated for every pair of variants with a maximum of 100 000 variants between variant pairs (--ld-window

100000), ≤4000 kb apart (--ld-window-kb 4000), and with

r 2 value > 0.01 (--ld-window- r 2 0.01) within populations. Genomic inbreeding was evaluated using ROH with the --homozyg command. ROH were de?ned as tracts of homozygous genotypes that were >1000 kb in length identi?ed for 1 SNP per

120 kb on average and 2 consecutive SNPs <1000 kb apart. No

more than 2 missing genotypes and 1 possible heterozygous geno - type were allowed. The following parameters were set: --homozyg; --homozyg-kb 1000; --homozyg-snp 50; --homozyg-gap 1000; --homozyg-window-het 1; --homozyg-density 120; --homozyg- window-missing 2. The individual sum of total ROH per animal was calculated initially. The F ROH statistic proposed by (

McQuillan

et al. 2008 ) was then calculated, whereby the total length of ROH covering an individual animal's genome ( L ROH ) is divided by the length of the autosomal genome ( L AUTO ); F ROH = L ROH / L AUTO . Here, we used the length of the equine autosomal genome as 2 242 960 kb ( www.ncbi.nlm.nih.gov/genome/145?genome_assembly_id=22878 ).

Interbreed Analyses

PCA and the estimation of pairwise

F ST values were conducted using smartPCA from the EIGENSOFT package (version 4.2) (

Price et al. 2006

). To include all the horses in the analysis, an option, outliersigmathresh,

was set as 10. All other parameters were set to default values.Maximum likelihood (ML) topology trees were generated using

the TreeMix (version 1.12) software package (

Pickrell and Pritchard

2012
). The Przewalski's horse population ( n = 17) was used as an outgroup. TreeMix was run without sample size correction (-noss), but with 50 SNP blocks (-k 50) as described in the TreeMix soft - ware documentation. TreeMix was used ?rst with all 40 popula - tions to visualize the general branching patterns in the topology generated. We next built a ML tree of 33 established populations excluding breeds (FM, FT, HAN, MARM, PT, QH, and SZWB) that are at a relatively early stage in the process of breed development or that have been in?uenced by gene ?ow from the Thoroughbred population. We then identify populations that are poor ?ts to the tree model and model migration events (migration edges) relating to these populations. The phylogenies and phylogenetic network were visualized by using an R script provided in the TreeMix software package. Population structure was investigated using the fastSTRUCTURE software package (

Raj et al. 2014

) with the simple prior approach.

To identify the "true"

K value for the number of ancestral popula - tions, a series of fastSTRUCTURE runs with prede?ned K values were examined using the chooseK.py script. Outputs from the fastSTRUCTURE analyses were visualized using the DISTRUCT software program (

Rosenberg et al. 2002

).

Results

Sample Summary FollowingfiQC

Following QC, 358 824 SNPs were used for samples from 5 Chinese

Mongolian horse populations (

n = 100). A total of 35 414 SNPs were used for 41 populations ( n = 912). Przewalski's horse represents a population that did not contribute genetically to modern horses and was therefore used in the TreeMix analysis as an outgroup. Table 1 provides detailed information on the horse populations.

Intrapopulation Analyses

The individual inbreeding coef?cient (

F ) and expected heterozy - gosity ( H e ) were estimated for each population using 14 012 auto - somal SNPs that remained after pruning for MAF, genotyping rate, and LD. The F and H e values for each population are provided in

Table 1

. The highest F and the lowest H e values were observed in the CLYD, EXMR, and MNGP breeds, while the landrace populations (AB, BCIH, CSP, MON, NFST, SH, TUVA, WS, and WZ) and the populations that have had appreciable admixture (HAN, MARM,

PT, and QH) exhibited the lowest

F and the highest H e values. Using the 35 414 autosomal SNPs, pruned only for MAF and genotyping rate, the level of genomic inbreeding was calculated using the ROH method. Figure 2 illustrates the total combined summed ROH per individual genome estimated for a selection of populations, including 7 Asian populations (AB, BCIH, MON, SH, TUVA, WS, and WZ) and 5 other representative comparator popula - tions (CLYD, FIN, NORF, QH, and TB). There were distinguishable patterns with respect to breed origins. Genome-wide summed ROH was observed to be lowest in the 7 Asian populations (AB, BCIH, MON, SH, TUVA, WS, and WZ). Conversely, the CLYD and TB breeds exhibited the highest sums of ROH for individual animals.

Furthermore, the mean

F ROH statistics estimates also re?ected the summed ROH results described above. As shown in

Supplementary

Table S2

, mean F ROH was lowest in the BCIH, MON, SH, WS, and

WZ populations and the highest mean

F ROH values were observed in CLYD and TB breeds. Additional summary ROH data are provided in

Supplementary Table S2

. 772

Journal of Heredity, 2019, Vol. 110, No. 7Downloaded from https://academic.oup.com/jhered/article/110/7/769/5601164 by guest on 28 August 2023

Population-speci?c decay of LD essentially mirrored the results of the F ROH analysis (

Supplementary Figure S2

).

Interpopulation Analyses

PCA and

F ST Results PCA was performed using the 2 different data sets, 5 Chinese Mongolian populations and 40 domestic horse populations. A PCA based on a set of SNPs (358 824 SNPs) for the 5 Chinese Mongolian populations was performed to investigate ?ne-scale population struc - ture and genomic diversity. PC1 and PC2 accounted for 15.58% and

13.44% of the total variation, respectively (

Figure 3

) and PC1 distin - guished AB from the other populations. For the PC2 axis, BCIH was partitioned into 2 clusters, one of which was similar to SH, WS, and WZ. While SH is an improved breed it was not distinguishable from

WS and WZ for either PC1 or PC2.

PCA was also performed to examine genome-wide SNP vari -

ation in the 40 populations. Based on geographic origin, the 40 populations grouped into 5 groups in total, one of them is the Asian group having 7 populations (5 Chinese Mongolian popu-lations and 2 Asian populations available using publicly avail-able data sets - MON and TUVA) (Table 1). The results of the PCA for PC1 and PC2 are plotted in Figure 4. Here, the PC1 axis, which accounted for 33.61% of total variance, separated the TB from the other populations. Horses from Europe and Middle East/Iberia were distinguished from the TB cluster on the PC1 axis. The majority of Asian animals were centrally located, for both the PC1 and PC2 axes, and were closest to Scandinavian populations (FIN, ICE, NORF, and NSWE). On the PC2 axis, which explained 12.68% of the total variance, 2 European populations (SHR and CLYD) were clearly separated from the other populations. The Asian populations were also located centrally in this dimension. Based on the PC1 and PC2 dimensions, animals from WS, WZ, and MON clustered together, and horses from the FIN breed were also located within this cluster. Individuals from SH, BCIH, and

Table 1. Horse populations, population code, sample size, mean individual inbreeding coefcient (F), expected heterozygosity (He), and

geographical groupings

PopulationsCodeSample size (n)FH

e

Geographical groupings

Przewalski's horsePRZE17

Abaga BlackAB150.0510.344Asia (ASIA)

Akhal TekeAKTK190.1030.326Middle East/Iberia (MIDI) AndalusianAND180.1280.316Middle East/Iberia (MIDI)

ArabianARR240.1280.314Middle East/Iberia (MIDI)

Baicha Iron HoofBCIH190.0380.349Asia (ASIA)

BelgianBEL300.1190.319Europe (EURO)

ClydesdaleCLYD240.2690.265Europe (EURO)

CaspianCSP180.050.344Middle East/Iberia (MIDI)

ExmoorEXMR240.2480.272Europe (EURO)

Fell PonyFELL210.1260.317Europe (EURO)

FinnhorseFIN270.0610.340Europe (EURO)

Florida CrackerFLCR70.1730.304Americas (AMER)

Franches-MontagnesFM190.0970.327Thoroughbred-derived (THOR) French TrotterFT170.1030.325Thoroughbred-derived (THOR) HanoverianHAN150.0410.347Thoroughbred-derived (THOR)

IcelandicICE250.1070.324Europe (EURO)

LusitanoLUST240.0990.327Middle East/Iberia (MIDI)

MaremmanoMARM240.0330.350Thoroughbred-derived (THOR)

MiniatureMINI210.0920.329Americas (AMER)

Mangalarga PaulistaMNGP150.2440.274Thoroughbred-derived (THOR)

MongolianMON190.0260.353Asia (ASIA)

MorganMOR400.0950.328Americas (AMER)

New Forest PonyNFST150.0380.349Europe (EURO)

Norwegian FjordNORF210.1360.313Europe (EURO)

North Swedish HorseNSWE190.1430.311Europe (EURO)

PercheronPERC230.0980.327Europe (EURO)

Peruvian PasoPERU210.0650.338Americas (AMER)

Puerto Rican Paso FinoPRPF200.1150.320Americas (AMER)

PaintPT250.0350.349Thoroughbred-derived (THOR)

Quarter HorseQH400.0450.346Thoroughbred-derived (THOR)

SaddlebredSB250.0990.326Americas (AMER)

SanheSH230.0160.357Asia (ASIA)

ShetlandSHET270.1960.291Europe (EURO)

ShireSHR230.1610.304Europe (EURO)

StandardbredSTBD400.1370.313Americas (AMER)

Swiss WarmbloodSZWB140.0490.344Thoroughbred-derived (THOR) ThoroughbredTB360.1230.317Thoroughbred-derived (THOR)

TuvaTUVA150.0410.348Asia (ASIA)

WushenWS220.0580.341Asia (ASIA)

WuzhumuqinWZ210.0480.345Asia (ASIA)

Total41912

Journal of Heredity,

2019, Vol. 110, No. 7

773Downloaded from https://academic.oup.com/jhered/article/110/7/769/5601164 by guest on 28 August 2023

AB populations were not clustered together and overlapped with some CSP horses.

Pairwise F

ST values were calculated for the 40 populations (

Supplementary Table S3

). As well as the lowest level of differenti - ation between the PT and QH populations ( F ST = 0.002), relatively low F ST values were also observed among the Asian populations. For example, the lowest level of genetic differentiation ( F ST = 0.005) was observed between WS and WZ as well as MON and WS.

Furthermore, low values of

F ST were also detected between the Asian populations and other landrace populations (CSP and NFST), as well as populations with high intrapopulation diversity (FIN).

Modeling Topologies of Relationships Between

Populations

To explore the relationships among different horse populations, we ?rst generated a ML tree without migration edges using the 40 popu - lations, which also included the Przewalski's horse population ( n =17) as an outgroup (

Figure 5

). The 7 Asian populations (AB, BCIH, MON, SH, TUVA, WS, and WZ) are closely related based on the pos - ition of these populations in the ML tree. However, these populations were not part of a discrete clade and emerged as a series of branches between the 2 large groups evident in the PCA plot shown in

Figure

4. The clade consisting of MON, WS, and WZ populations formed a

short branch at the base of the tree and the AB, BCIH, SH, and TUVA populations formed individual short branches. On the other hand, 2 most drifted populations, TB and CLYD had the longest branches. To focus the investigation on gene ?ow between Chinese Mongolian horses and other populations, we ?rst built a ML tree

using 33 established populations including the Przewalski's horse population as an outgroup (Figure 6). Populations (i.e., AB, AKTK, ARR, CSP, MNGP, MOR, and SB) that are poor ?t to tree model were identi?ed by visualize the residues from the ?t of the model

to the data (Supplementary Figure S3). We then sequentially added migration events to the tree. In Figure 6, we show the inferred graph with 6 migration edges. The ?rst migration edge (m = 1) explained modeled gene ?ow from the ARR population to the CSP popula-tion. As the number of modeled migration events increased, gene ?ow from the root of the TB and ARR populations to the Chinese Mongolian breeds (i.e., AB, BCIH, and SH) was also evident.

Population Structure Analysis

An analysis of population structure among the 40 populations showed clustering of the 7 Asian populations (AB, BCIH, MON, SH,

TUVA, WS, and WZ) throughout the values of

K = 40. K = 22-29 emerged as the most likely number of "true" clusters, corresponding to the likely number of ancestral populations. All samples from the 7 Asian populations clustered together throughout the all values of K . Evidence of clustering of the Asian populations (hereafter the Asian cluster) and several populations was observed throughout the all values of K (

Figure 7

and

Supplementary Table S4

). Clustering of the Asian horses with FIN horse was observed until K = 20. From this point, the FIN breed represented a distinct cluster until K = 40, with the exception of K = 27 and 29 at which point the FIN remained in the Asian cluster. Compared to the FIN breed, clustering of NFST with the Asian populations was observed up to K = 40 albeit at a moderate portion of assignment. Similar to the NFST, CSP appears to be in the Asian cluster throughout K = 40 with moderate frequencies. Finally, a small proportion of assignment of

Figure 2. Individual sums of total ROH per animal in 7 Asian populations (AB, BCIH, MON, SH, TUVA, WS, and WZ) and 5 other comparator breeds (CLYD, FIN,

NORF, QH, and TB). See online version for full colors.

774 Journal of Heredity, 2019, Vol. 110, No. 7Downloaded from https://academic.oup.com/jhered/article/110/7/769/5601164 by guest on 28 August 2023

TB-crossed populations (HAN, MARM, PT, QH, and SZWB) to the

Asian cluster was observed across a range of

K values.

Discussion

We have performed an analysis based on genome-wide SNP data generated for 5 distinct Chinese Mongolian horse populations (total sample size of 100 with a mean of n = 20 for each population). These data were integrated with SNP data from 35 widely distrib - uted modern horse populations to provide a global picture of horse genomic diversity. By comparing intrapopulation genetic diversity and establishing relationships with other modern breeds, we have, for the ?rst time, comprehensively described sequence diversity in

Chinese Mongolian horses.

Our newly sampled horses were genotyped using the higher density SNP genotyping array (670K) and then downsampled to

35K SNPs to integrate into public data. However, the potential im

- pact of SNP ascertainment bias on diversity calculations must still needed to be considered since the initial horse SNP discovery for SNP variants included on the genotyping arrays focused on a rela - tively limited panel of breeds (

McCue et al. 2012

; Schaefer et al. 2017
). In order to mitigate ascertainment bias as much as possible, F , H e , and F ST were estimated using genotypes following stringent LD pruning since a previous study suggested that an LD-based pruning

strategy could help to produce results which are most comparable to those obtained from whole genome sequencing data (Malomane et al. 2018).

PCA plots, modeling topologies, population structure analysis, and estimation of F ST values were used to reconstruct genetic af?n - ities among Chinese Mongolian horses and other modern popula - tions, which largely re?ected historic and established relationships. Our results indicate that the Chinese Mongolian populations, as well as the MON and TUVA populations, are closely related, but genet - ically distinct from other global populations. Low divergence was observed between the Asian populations and other global landrace populations. In addition, the in?uence of Thoroughbred (TB) and TB-related populations on Chinese Mongolian horses was observed and this was also observed in our Y chromosome diversity analysis (

Han et al. 2019

).

Asian Horse Populations

In addition to the 5 Chinese Mongolian horse populations sampled for the present study, 2 other Asian populations were also available from the publicly available Illumina data set (MON and TUVA). Therefore, we considered these 7 populations as an Asian horse population group. Because of their geographic location, it is con - sidered likely that Asian horse populations most closely resemble early domestic horses that subsequently spread across Eurasia (

Bjørnstad et al. 2003

; Hendricks 2007; Warmuth et al. 2013). These Asian populations exhibited close genetic relationships based on fi0.4fi0.3fi0.2fi0.10.0 0.30.20.10.00.1

PC1 (15.58%)

PC2 (13.44%)

Breed AB BCIH SH WS WZ Figure 3. Principal component analysis (PCA) plot for PC1 and PC2 calculated for horses from five Chinese Mongolian populations (n = 100). The PCA was

performed using data from 358 824 SNPs. X and Y axes show the relative variance contributions for the PC1 and PC2. See online version for full colours.

Journal of Heredity,

2019, Vol. 110, No. 7

775Downloaded from https://academic.oup.com/jhered/article/110/7/769/5601164 by guest on 28 August 2023

the results obtained using PCA, hierarchical clustering/structure as- sessment, and ML topology tree; they also displayed relatively low intrapopulation divergence as estimated using F ST values. These re - sults indicate that these horse populations have a shared ancient ancestry, which has been preserved for thousands of years. Despite some gene ?ow, they have not been signi?cantly perturbed by gene ?ow from breeds elsewhere in Eurasia, which is also evident for our Y chromosome data (

Han et al. 2019

). However, these popula - tions are genetically differentiated to a similar level as some distinct modern breeds. For example, an F ST value of 0.037 was observed between the MINI and SHET breeds, which is comparable to the F ST value of 0.038 between the AB and BCIH populations. In add - ition, by comparison, the closest relationship observed among the global breeds was for the PT and QH populations ( F ST = 0.002); this compares to pairwise F ST values of 0.005 for the most closely re - lated of the Asian populations: WS versus WZ and WS versus MON (

Supplementary Table S3

).

Asian Populations and Other LandraceBreeds

Genetic similarities between the Asian populations and FIN as well as 2 other landrace populations (CSP and NFST) were ob - served. While FIN horses are unlikely to have originated from na - tive Finnish wild horses, there is some historical evidence that FIN horses have explicit Mongolian origins (

Kantanen 2014

). Our data support this hypothesis because the pairwise F ST

values between the FIN and the Asian populations indicate less divergence than that observed between FIN and other northern European populations of similar cultural and geographical provenance (i.e., the Scandinavian ICE, NORF, and NSWE breeds). In addition, the in?uence of the Asian populations on the FIN horse ranged from 23.5% to 97.9% until K = 20 (Supplementary Table S4). This close relationship was

also supported by the PCA analysis (

Figure 4

), as some of the FIN horses clustered with WS, WZ, and MON. These results suggest that contemporary FIN horses may indeed have had direct ancestors of

Mongolian origin.

Although the ARR and AKTK populations are geographically proximal to the CSP population, the CSP breed has more genetic af?nities with Asian horse populations. As well as low pairwise F ST values observed between CSP and the Asian populations, a relatively high percentage assignment of CSP to the Asian cluster was observed in the genetic structure analysis (

Supplementary Table S4

). These re - sults are likely due to an ancient shared ancestry. The CSP popula - tion originates in Iran and is believed to be one of the oldest horse breeds in the world with a history that has been traced back to 5000 BP (

Shasavarani and Rahimi-Mianji 2010

). Similarly, low divergence and a close genetic relationship between the NFST and the Asian populations was observed. The NFST is rec - ognized as an old population with high levels of diversity and is gen - erally free-ranging in Great Britain (

Petersen et al. 2013a

). While there is no literature documenting the direct genetic relationship between 0.1

5 0.10 0.050.000.05

0.1

00.050.000.05

PC1 (33.61%)

PC2 (12.68%)

Breed FLCR MINI MNGP MOR PERU PRPF SB STBD AB BCIH MON SH TU VA WS WZ BE L CLYD EXMR FELL FI NICE NFST NORF NSWE PERC SHET SHR AKT K AND ARR CS P LUST FM FT HA N MARM PT QH SZWB TB Figure 4. Principal component analysis (PCA) plot for PC1 and PC2 comprising 40 p opulations split into ve groups according to geographic origin of the

populations. The PCA was performed using data from 35 414 SNPs and individuals from the same group are represented by the same shape. X and Y axes show

the relative variance contributions for the PC1 and PC2. See online version for full colours.

776 Journal of Heredity, 2019, Vol. 110, No. 7Downloaded from https://academic.oup.com/jhered/article/110/7/769/5601164 by guest on 28 August 2023

the Asian populations and NFST, it has been reported that domestic horses in UK may have been imported from the eastern steppes, the

Iberian refugium, or both (

Warmuth et al. 2011

). Therefore, this ob - servation may indicate that the NFST were imported from the eastern steppes where the ancestors of the Chinese Mongolian horse lived.

Asian Populations and TB/TB-Derived Populations

A small proportion of TB-crossed populations assigned to the Asian populations (i.e., AB, BCIH, and SH) was observed throughout the K values examined in the analysis of population structure (

Figure 7

and

Supplementary Table S4

). Along with the observation of migration edges generated from the topologies (

Figure 6

), all of male horses from AB, BCIH, and SH (except for 1 SH) carry haplotypes represent most stallion lines in modern horse breeds (

Han et al. 2019

). These patterns likely re?ect recent admixture due to the recent importation of TB and/or ARR to improve local herds in Inner Mongolia. There was no observation of this gene ?ow in the other Chinese Mongolian populations (WS, WZ), which is in agreement with our

Y chromosome data (

Han et al. 2019

) and local knowledge that indicates that these populations have not been subjected to breed improvement from exotic stock.

Within-Breed Diversity in Chinese

Mongolian Horses

The highest genetic diversity (

H e ) was observed in the 5 Chinese Mongolian populations (AB, BCIH, WS, WZ, and SH), 2 other Asian populations (MON and TUVA), and 2 old landrace popula - tions (CSP and NFST). High genetic diversity was also observed in breeds that are in the early stages of breed development and allowed to crossbreed, such as QH and PT. In contrast, low genetic diversity was observed in populations that have undergone severe population bottlenecks, such as the CLYD and EXMR. While the 5 Chinese Mongolian horse populations are generally described as "Mongolian," they appear to have different demographic histories. Therefore, we evaluated genetic diversity separately within each population. Among them, the highest diversity was observed in the SH. Demographic factors including a large population size, a di - versity of founding stock, and substantial phenotypic diversity within Thoroughbred (TB, THOR)Andalusian (AND, MIDI)Puerto Rican Paso Fino (PRPF, AMER )

Arabian (ARR, MIDI

)Fell Pony (FELL, EURO)

Przewalski's horse (PRZE)

New Forest Pony (NFST, EURO

) Franches-Montagnes (FM, THOR)Percheron (PERC, EURO)

Swiss Warmblood (SZWB, THOR)

Mangalarga Paulista (MNGP, THOR

)French Troer (FT, THOR)North Swedish Horse (NSWE, EURO)

Miniature (MINI, AMER

)

Exmoor (EXMR, EURO)Icelandic (ICE, EURO)

Maremmano (MARM, THOR)Lusitano (LUST, MIDI)Tuva (TUVA, ASIA)

Morgan (MOR, AMER

)Shetland (SHET, EURO)

Quarter Horse (QH, THOR)

Paint (PT, THOR)

Wuzhumuqin (WZ, ASIA)Norwegian Fjord (NORF, EURO)

Hanoverian (HAN, THOR

)

Clydesdale (C

LYD, EURO)

Sanhe (SH, ASIA)

Mongolian (MON, ASIA)

Abaga Black (AB, ASIA)Shire (SHR, EURO)

Standardbred (STBD, AMER)Caspian (CSP, MIDI)Baicha Iron Hoof (BCIH, ASIA)Belgian (BEL, EURO)

Akhal Teke (AKTK, MIDI)Finnhorse (FIN, EURO)

Peruvian Paso (PERU, AMER)Wushen (WS, ASIA)

Florida Cracker (FLCR, AMER)Saddlebred (SB, AMER)

Dri parameter

0.000.010.020.03

10 s.e.

Figure 5. Maximum likelihood (ML) tree with no migration edges for 40 horse populations. The ML tree was generated using 35 414 autosomal SNPs and the

Przewalski"s horse (PRZE) is included as an outgroup. Horizontal branch lengths are proportional to the amount of genetic drift that has occurred on the branch.

The scale bar shows 10 times the mean standard error of the estimated entries in the sample covariance matrix. See online version for full colours.

Journal of Heredity,

2019, Vol. 110, No. 7

777Downloaded from https://academic.oup.com/jhered/article/110/7/769/5601164 by guest on 28 August 2023

the population contribute to the high level of genetic diversity ob- served in this population. Since the 20th century SH horses have been developed from local Mongolian stock and outcrossed with a range of other horse breeds (e.g., Orlov, Anglo-Norman horse, Anglo-Arabian horse, Arabian horse, Shetland, Thoroughbred, and Standardbred). Furthermore, in 1955 it was reported that subpopulation structure was evident in this population based on physical characteristics, and light, medium, and heavy horses were described (

National Livestock

and Poultry Genetic Resources Committee 2011 ). However, such population substructure was not detectable in our sample ( n = 21).

The relatively high

H e , low

F, short blocks of LD, and a low frequency

of ROH re?ect recent outcrossing in SH. WZ and WS are the 2 principal subtypes of Chinese Mongolian

horses and have relatively high census population sizes: 24 587 for WZ and 5000 for WS in 2005, respectively (National Livestock and Poultry Genetic Resources Committee 2011). They were the 2 popu-lations most similar to the previously genotyped MON population. Since these horse populations have been free-ranging, are less man-aged than modern horse breeds, and are not subject to intense selective breeding, high within-breed genetic diversity would be expected.

The BCIH is a subtype of Chinese Mongolian horse known for its tough hooves adapted for traversing rugged mountain terrain. They are therefore also known as "Iron Hoof" horses. This subtype was thought to be extinct since horse power was supplanted in the early

20th century and the number of BCIH decreased to fewer than 100 in

2005 (

National Livestock and Poultry Genetic Resources Committee 2011
). Recently, a number of BCIH were "rediscovered" by local herdsmen and identi?ed based on their physical characteristics and

Figure 6. Maximum likelihood (ML) tree with migration edges for 33 populations and the Przewalski's horse population as an outgroup. This topology was

generated with six migration edges with the second migration edge (labe lled 2) hypothesized to be gene flow from ancestors of TB and ARR into the ancestry

to the SH and BCIH. The sixth migration edge (labelled 6) is hypothesized to be another gene flow event from the ancestral population of the

TB and ARR into

the lineage ancestral to the AB. Migration arrows are coloured according to their weight. Horizontal branch lengths are proportional to the amount of genetic

drift that has occurred on the branch. The scale bar shows 10 times the mean standard error of the estimated entries in the sample covariance matrix. See online

version for full colours. 778

Journal of Heredity, 2019, Vol. 110, No. 7Downloaded from https://academic.oup.com/jhered/article/110/7/769/5601164 by guest on 28 August 2023

Downloaded from https://academic.oup.com/jhered/article/110/7/769/5601164 by guest on 28 August 2023

Supplementary Material

Supplementary data are available at

Journal of Heredity

online.

Data Availability

SNPs for 100 Chinese Mongolian horses are available at Open

Science Framework (

https://osf.io/2xvqf/quickles ).

Funding

This work was supported by Science Foundation Ireland (grant number 11/PI/1166) and the China Scholarship Council.

Acknowledgments

We would like to thank the horse owners who kindly supported sample collection.

Ethics Statement

Samples were collected with informed, signed, owners' consent. Approval for collection and the movement of genetic material for the research was granted by University College Dublin Animal Research Ethics Committee, the Inner Mongolia Entry-Exit Inspection and Quarantine of the People's Republic of China, the Inner Mongolia Agricultural University, and the Department of Agriculture, Ireland.

References

Anthony DW. 2010.

The horse, the wheel and language: how Bronze-Age riders from the Eurasian steppes shaped the modern world . Princeton (NJ): Princeton University Press. Bjørnstad G, Nilsen NØ, Røed KH. 2003. Genetic relationship between Mon - golian and Norwegian horses?

Anim Genet

. 34:55-58.

Chang H. 2009.

Animal genetic resources

. Beijing: Science Press. Chang CC, Chow CC, Tellier LC, Vattikuti S, Purcell SM, Lee JJ. 2015. Second-generation PLINK: rising to the challenge of larger and richer datasets.

Gigascience

. 4:7. Cieslak M, Pruvost M, Benecke N, Hofreiter M, Morales A, Reissmann M, Ludwig A. 2010. Origin and history of mitochondrial DNA lineages in domestic horses.

PLoS One

. 5:e15311. Danecek P, Auton A, Abecasis G, Albers CA, Banks E, DePristo MA,

Handsaker RE, Lunter G, Marth GT, Sherry ST,

etal. ; 1000 Genomes Project Analysis Group. 2011. The variant call format and VCFtools. Bio - informatics . 27:2156-2158. Fages A, Hanghøj K, Khan N, Gaunitz C, Seguin-Orlando A, Leonardi M, McCrory Constantz C, Gamba C, Al-Rasheid KAS, Albizuri S, etal. 2019. Tracking ve millennia of horse management with extensive ancient genome time series. Cell . 177:1419-35 e31. Gaunitz C, Fages A, Hanghøj K, Albrechtsen A, Khan N, Schubert M, Seguin-

Orlando A, Owens IJ, Felkel S, Bignon-Lau O,

etal. 2018. Ancient gen
- omes revisit the ancestry of domestic and Przewalski's horses.

Science

.

360:111-114.

Han H, Wallner B, Rigler D, MacHugh DE, Manglai D, Hill EW. 2019. Chinese Mongolian horses may retain early domestic male genetic lineages yet to be discovered. LID . doi:10.1111/age.12780

Hendricks BL. 2007.

International encyclopedia of horse breeds

. Norman (OK): University of Oklahoma Press. Hill EW, Gu J, Eivers SS, Fonseca RG, McGivney BA, Govindarajan P, Orr N, Katz LM, MacHugh DE, MacHugh D. 2010. A sequence polymorphism in MSTN predicts sprinting ability and racing stamina in thoroughbred horses.

PLoS One

. 5:e8645.Kantanen J. 2014. Mongolian horse is genetically closer relative to Finnhorse than other Nordic horse breeds. Available from: URL https:// finnhorseblog.com/2014/09/27/mongolian-horse-is-genetically-closer- relative-to-nnhorse-than-other-nordic-horse-breeds/ . Larson G, Fuller DQ. 2014. The evolution of animal domestication.

Annu Rev

Ecol Evol Syst

. 45:115-136. Larson G, Piperno DR, Allaby RG, Purugganan MD, Andersson L, Arroyo- Kalin M, Barton L, Vigueira CC, Denham T, Dobney K, etal. 2014. Cur-
rent perspectives and the future of domestication studies.

Proc Natl Acad

Sci USA

. 111:6139-6146. Li H, Handsaker B, Wysoker A, Fennell T, Ruan J , Homer N, Marth G, Abecasis G, Durbin R; 1000 Genome Project Data Processing Subgroup.

2009. The sequence alignment/map format and SAMtools.

Bioinformatics

.

25:2078-2079.

Librado P, Gamba C, Gaunitz C, Der Sarkissian C, Pruvost M, Albrechtsen A,

Fages A, Khan N, Schubert M, Jagannathan V,

et al. 2017. Ancient
genomic changes associated with domestication of the horse.

Science

.

356:442-445.

Lira J, Linderholm A, Olaria C, Brandström Durling M, Gilbert MT, Ellegren H, Willerslev E, Lidén K, Arsuaga JL, Götherström A. 2010. An - cient DNA reveals traces of Iberian Neolithic and Bronze Age lineages in modern Iberian horses.

Mol Ecol

. 19:64-78. Malomane DK, Reimer C, Weigend S, Weigend A, Shari AR, Simianer H.

2018. Efciency of different strategies to mitigate ascertainment bias whe

n using SNP panels in diversity studies.

BMC Genomics

. 19:22. McCue ME, Bannasch DL, Petersen JL, Gurr J, Bailey E, Binns MM, Distl O, Guerin G, Hasegawa T, Hill EW, Leeb T, Lindgren G, Penedo MC, Roed KH, Ryder OA, Swinburne JE, Tozaki T, Valberg SJ, Vaudin M, Lindblad-Toh K, Wade CM, Mickelson JR 20 12. A high density SNP array for the domestic horse and extant Perissodactyla: utility for asso - ciation mapping, genetic diversity, and phylogeny studies. Plos Genetics.

8:e1002451.

McQuillan R, Leutenegger AL, Abdel-Rahman R, Franklin CS, Pericic M, Barac-Lauc L, Smolej-Narancic N, Janicijevic B, Polasek O, Tenesa A, etal.

2008. Runs of homozygosity in European populations.

Am J Hum Genet

.

83:359-372.

National Livestock and Poultry Genetic Resources Committee. 2011.

Annals

of livestock and poultry genetic resources in China (donkeys, horses and camels) . Beijing (China): China Agricultural Press. Outram AK, Stear NA, Bendrey R, Olsen S, Kasparov A, Zaibert V, Thorpe N, Evershed RP. 2009. The earliest horse harnessing and milking.

Science

.

323:1332-1335.

Petersen JL, Mickelson JR, Cothran EG, Andersson LS, Axelsson J, Bailey E,

Bannasch D, Binns MM, Borges AS, Brama P,

etal. 2013a. Genetic diver-
sity in the modern horse illustrated from genome-wide SNP data. PLoS One . 8:e54997. Petersen JL, Mickelson JR, Rendahl AK, Valberg SJ, Andersson LS, Axelsson J,

Bailey E, Bannasch D, Binns MM, Borges AS,

etal. 2013b. Genome-wide
analysis reveals selection for important traits in domestic horse breeds .

PLoS Genet

. 9:e1003211. Pickrell JK, Pritchard JK. 2012. Inference of population splits and mixtures from genome-wide allele frequency data.

PLoS Genet

. 8:e1002967. Price AL, Patterson NJ, Plenge RM, Weinblatt ME, Shadick NA, Reich D.

2006. Principal components analysis corrects for stratication in genome-

wide association studies.

Nat Genet

. 38:904-909. Raj A, Stephens M, Pritchard JK. 2014. fastSTRUCTURE: variational infer- ence of population structure in large SNP data sets. Genetics. 197:573- 589.
Rosenberg NA, Pritchard JK, Weber JL, Cann HM, Kidd KK, Zhivotovsky LA, Feldman MW. 2002. Genetic structure of human populations.

Science

.

298:2381-2385.

Schaefer RJ, Schubert M, Bailey E, Bannasch DL, Barrey E, Bar-Gal GK, Brem G, Brooks SA, Distl O, Fries R, etal. 2017. Developing a 670k genotyping array to tag ~2M SNPs across 24 horse breeds.

Bmc Genomics

18: 565.

780

Journal of Heredity, 2019, Vol. 110, No. 7Downloaded from https://academic.oup.com/jhered/article/110/7/769/5601164 by guest on 28 August 2023

Shasavarani H, Rahimi-Mianji G. 2010. Analysis of genetic diversity and es- timation of inbreeding coef?cient within Caspian horse population using microsatellite markers.

Afr J Biotechnol

. 9:293-299. Wallner B, Palmieri N, Vogl C, Rigler D, Bozlak E, Druml T, Jagannathan V, Leeb T, Fries R, Tetens J, et al. 2017. Y chromosome uncovers the recent oriental origin of modern stallions.

Curr Biol

. 27:2029-2035.e5.Warmuth V, Eriksson A, Bower MA, Cañon J, Cothran G, Distl O, Glowatzki-

Mullis ML, Hunt H, Luís C, do Mar Oom M,

et al. 2011. European do
- mestic horses originated in two holocene refugia.

PLoS One

. 6:e18194. Warmuth V, Manica A, Eriksson A, Barker G, Bower M. 2013. Autosomal gen - etic diversity in non-breed horses from eastern Eurasia provides insight s into historical population movements.

Anim Genet

. 44:53-61.

Journal of Heredity,

2019, Vol. 110, No. 7

781Downloaded from https://academic.oup.com/jhered/article/110/7/769/5601164 by guest on 28 August 2023