Invited review: The future of selection decisions and breeding programs: What are we breeding for, and who decides?

Genetic selection has been a very successful tool for the long-term improvement of livestock populations, and the rapid adoption of genomic selection over the last decade has doubled the rate of gain in some populations. Breeding programs seek to identify genetically superior parents of the next generation, typically as a function of an index that combines information about many economically important traits into a single number. In the United States, the data that drive this system are collected through the national dairy herd improvement program that began more than a century ago. The resulting information about animal performance, pedigree, and genotype is used to compute genomic evaluations for comparing and ranking animals for selection. However, the full expression of genetic potential requires that animals are placed in environments that can support such performance. The Agricultural Research Service of the US Department of Agriculture and the Council on Dairy Cattle Breeding collaborate to deliver state-of-the-art genomic evaluations to the dairy industry. Today, most breeding stock are selected and marketed using the net merit dollars (NM$) selection index, which evolved from 2 traits in 1926 (milk and fat yield) to a combination of 36 individual traits following the last NM$ update in 2018. Updates to NM$ require the estimation of many different values, and it can be difficult to achieve consensus from stakeholders on what should be added to, or removed from, the index at each review, and how those traits should be weighted. Over time, the majority of the emphasis in the index has shifted from yield traits to fertility, health, and fitness traits. Phenotypes for some of these new traits are difficult or expensive to measure, or require changes to on-farm habits that have not been widely adopted. This is driving interest in sensor-based systems that provide continuous measurements of the farm environment, individual animal performance, and detailed milk composition. There is also a need to capture more detailed data about the environment in which animals perform, including information about feeding, housing, milking systems, and infectious and parasitic load. However, many challenges accompany these new technologies, including a lack of standardization or validation, need for high-speed internet connections, increased computational requirements, and interpretations that are often not backed by direct observations of biological phenomena. This work will describe how US selection objectives are developed, as well as discuss opportunities and challenges associated with new technologies for measuring and recording animal performance.

Exome sequences and multi-environment field trials elucidate the genetic basis of adaptation in barley.

Broadening the genetic base of crops is crucial for developing varieties to respond to global agricultural challenges such as climate change. Here, we analysed a diverse panel of 371 domesticated lines of the model crop barley to explore the genetics of crop adaptation. We first collected exome sequence data and phenotypes of key life history traits from contrasting multi-environment common garden trials. Then we applied refined statistical methods, including some based on exomic haplotype states, for genotype-by-environment (G×E) modelling. Sub-populations defined from exomic profiles were coincident with barley's biology, geography and history, and explained a high proportion of trial phenotypic variance. Clear G×E interactions indicated adaptation profiles that varied for landraces and cultivars. Exploration of circadian clock-related genes, associated with the environmentally adaptive days to heading trait (crucial for the crop's spread from the Fertile Crescent), illustrated complexities in G×E effect directions, and the importance of latitudinally based genic context in the expression of large-effect alleles. Our analysis supports a gene-level scientific understanding of crop adaption and leads to practical opportunities for crop improvement, allowing the prioritisation of genomic regions and particular sets of lines for breeding efforts seeking to cope with climate change and other stresses.

Genome-wide patterns of homozygosity provide clues about the population history and adaptation of goats.

BACKGROUND: Patterns of homozygosity can be influenced by several factors, such as demography, recombination, and selection. Using the goat SNP50 BeadChip, we genotyped 3171 goats belonging to 117 populations with a worldwide distribution. Our objectives were to characterize the number and length of runs of homozygosity (ROH) and to detect ROH hotspots in order to gain new insights into the consequences of neutral and selection processes on the genome-wide homozygosity patterns of goats. RESULTS: The proportion of the goat genome covered by ROH is, in general, less than 15% with an inverse relationship between ROH length and frequency i.e. short ROH (< 3 Mb) are the most frequent ones. Our data also indicate that ~ 60% of the breeds display low FROH coefficients (< 0.10), while ~ 30 and ~ 10% of the goat populations show moderate (0.10 < FROH < 0.20) or high (> 0.20) FROH values. For populations from Asia, the average number of ROH is smaller and their coverage is lower in goats from the Near East than in goats from Central Asia, which is consistent with the role of the Fertile Crescent as the primary centre of goat domestication. We also observed that local breeds with small population sizes tend to have a larger fraction of the genome covered by ROH compared to breeds with tens or hundreds of thousands of individuals. Five regions on three goat chromosomes i.e. 11, 12 and 18, contain ROH hotspots that overlap with signatures of selection. CONCLUSIONS: Patterns of homozygosity (average number of ROH of 77 and genome coverage of 248 Mb; FROH < 0.15) are similar in goats from different geographic areas. The increased homozygosity in local breeds is the consequence of their small population size and geographic isolation as well as of founder effects and recent inbreeding. The existence of three ROH hotspots that co-localize with signatures of selection demonstrates that selection has also played an important role in increasing the homozygosity of specific regions in the goat genome. Finally, most of the goat breeds analysed in this work display low levels of homozygosity, which is favourable for their genetic management and viability.

Functional SNP panel for parentage assessment and assignment in worldwide goat breeds.

BACKGROUND: International standard panels of single nucleotide polymorphisms (SNPs) have replaced microsatellites in several species for parentage assessment and assignment (PA) purposes. However, such a resource is still lacking in goats. The application of a cheap tool for PA would help the management of goat populations by improving the reliability of pedigree registration and, consequently, allow a better implementation of breeding schemes or conservation programs. RESULTS: Using data from the current GoatSNP50 chip, starting from a worldwide dataset of more than 4000 animals belonging to more than 140 breeds and populations from the AdaptMap initiative, we selected a panel of 195 SNPs. The assignment rate of this panel was up to 100% on an additional dataset that included 2000 Alpine and Saanen animals and highly related candidate sires. CONCLUSIONS: In this study, we defined a highly informative SNP panel, which will be publicly available to worldwide breeders and laboratories. Its development on such a large number of breeds and populations, together with validation on a second set of cosmopolitan breeds, makes it a promising and important genomic tool for the goat species.

Signatures of selection and environmental adaptation across the goat genome post-domestication.

BACKGROUND: Since goat was domesticated 10,000 years ago, many factors have contributed to the differentiation of goat breeds and these are classified mainly into two types: (i) adaptation to different breeding systems and/or purposes and (ii) adaptation to different environments. As a result, approximately 600 goat breeds have developed worldwide; they differ considerably from one another in terms of phenotypic characteristics and are adapted to a wide range of climatic conditions. In this work, we analyzed the AdaptMap goat dataset, which is composed of data from more than 3000 animals collected worldwide and genotyped with the CaprineSNP50 BeadChip. These animals were partitioned into groups based on geographical area, production uses, available records on solid coat color and environmental variables including the sampling geographical coordinates, to investigate the role of natural and/or artificial selection in shaping the genome of goat breeds. RESULTS: Several signatures of selection on different chromosomal regions were detected across the different breeds, sub-geographical clusters, phenotypic and climatic groups. These regions contain genes that are involved in important biological processes, such as milk-, meat- or fiber-related production, coat color, glucose pathway, oxidative stress response, size, and circadian clock differences. Our results confirm previous findings in other species on adaptation to extreme environments and human purposes and provide new genes that could explain some of the differences between goat breeds according to their geographical distribution and adaptation to different environments. CONCLUSIONS: These analyses of signatures of selection provide a comprehensive first picture of the global domestication process and adaptation of goat breeds and highlight possible genes that may have contributed to the differentiation of this species worldwide.

Genome-wide SNP profiling of worldwide goat populations reveals strong partitioning of diversity and highlights post-domestication migration routes.

BACKGROUND: Goat populations that are characterized within the AdaptMap project cover a large part of the worldwide distribution of this species and provide the opportunity to assess their diversity at a global scale. We analysed genome-wide 50 K single nucleotide polymorphism (SNP) data from 144 populations to describe the global patterns of molecular variation, compare them to those observed in other livestock species, and identify the drivers that led to the current distribution of goats. RESULTS: A high degree of genetic variability exists among the goat populations studied. Our results highlight a strong partitioning of molecular diversity between and within continents. Three major gene pools correspond to goats from Europe, Africa and West Asia. Dissection of sub-structures disclosed regional gene pools, which reflect the main post-domestication migration routes. We also identified several exchanges, mainly in African populations, and which often involve admixed and cosmopolitan breeds. Extensive gene flow has taken place within specific areas (e.g., south Europe, Morocco and Mali-Burkina Faso-Nigeria), whereas elsewhere isolation due to geographical barriers (e.g., seas or mountains) or human management has decreased local gene flows. CONCLUSIONS: After domestication in the Fertile Crescent in the early Neolithic era (ca. 12,000 YBP), domestic goats that already carried differentiated gene pools spread to Europe, Africa and Asia. The spread of these populations determined the major genomic background of the continental populations, which currently have a more marked subdivision than that observed in other ruminant livestock species. Subsequently, further diversification occurred at the regional level due to geographical and reproductive isolation, which was accompanied by additional migrations and/or importations, the traces of which are still detectable today. The effects of breed formation were clearly detected, particularly in Central and North Europe. Overall, our results highlight a remarkable diversity that occurs at the global scale and is locally partitioned and often affected by introgression from cosmopolitan breeds. These findings support the importance of long-term preservation of goat diversity, and provide a useful framework for investigating adaptive introgression, directing genetic improvement and choosing breeding targets.

Use of SNP genotypes to identify carriers of harmful recessive mutations in cattle populations.

BACKGROUND: SNP (single nucleotide polymorphisms) genotype data are increasingly available in cattle populations and, among other things, can be used to predict carriers of specific mutations. It is therefore convenient to have a practical statistical method for the accurate classification of individuals into carriers and non-carriers. In this paper, we compared - through cross-validation- five classification models (Lasso-penalized logistic regression -Lasso, Support Vector Machines with either linear or radial kernel -SVML and SVMR, k-nearest neighbors -KNN, and multi-allelic gene prediction -MAG), for the identification of carriers of the TUBD1 recessive mutation on BTA19 (Bos taurus autosome 19), known to be associated with high calf mortality. A population of 3116 Fleckvieh and 392 Brown Swiss animals genotyped with the 54K SNP-chip was available for the analysis. RESULTS: In general, the use of SNP genotypes proved to be very effective for the identification of mutation carriers. The best predictive models were Lasso, SVML and MAG, with an average error rate, respectively, of 0.2 %, 0.4 % and 0.6 % in Fleckvieh, and 1.2 %, 0.9 % and 1.7 % in Brown Swiss. For the three models, the false positive rate was, respectively, 0.1 %, 0.1 % and 0.2 % in Fleckvieh, and 3.0 %, 2.4 % and 1.6 % in Brown Swiss; the false negative rate was 4.4 %, 7.6 %1.0 % in Fleckvieh, and 0.0 %, 0.1% and 0.8 % in Brown Swiss. MAG appeared to be more robust to sample size reduction: with 25 % of the data, the average error rate was 0.7 % and 2.2 % in Fleckvieh and Brown Swiss, compared to 2.1 % and 5.5 % with Lasso, and 2.6 % and 12.0 % with SVML. CONCLUSIONS: The use of SNP genotypes is a very effective and efficient technique for the identification of mutation carriers in cattle populations. Very few misclassifications were observed, overall and both in the carriers and non-carriers classes. This indicates that this is a very reliable approach for potential applications in cattle breeding.

A practical approach to detect ancestral haplotypes in livestock populations.

BACKGROUND: The effects of different evolutionary forces are expected to lead to the conservation, over many generations, of particular genomic regions (haplotypes) due to the development of linkage disequilibrium (LD). The detection and identification of early (ancestral) haplotypes can be used to clarify the evolutionary dynamics of different populations as well as identify selection signatures and genomic regions of interest to be used both in conservation and breeding programs. The aims of this study were to develop a simple procedure to identify ancestral haplotypes segregating across several generations both within and between populations with genetic links based on whole-genome scanning. This procedure was tested with simulated and then applied to real data from different genotyped populations of Spanish, Fleckvieh, Simmental and Brown-Swiss cattle. RESULTS: The identification of ancestral haplotypes has shown coincident patterns of selection across different breeds, allowing the detection of common regions of interest on different bovine chromosomes and mirroring the evolutionary dynamics of the studied populations. These regions, mainly located on chromosomes BTA5, BTA6, BTA7 and BTA21 are related with certain animal traits such as coat colour and milk protein and fat content. CONCLUSION: In agreement with previous studies, the detection of ancestral haplotypes provides useful information for the development and comparison of breeding and conservation programs both through the identification of selection signatures and other regions of interest, and as indicator of the general genetic status of the populations.

SNPchiMp v.3: integrating and standardizing single nucleotide polymorphism data for livestock species.

BACKGROUND: In recent years, the use of genomic information in livestock species for genetic improvement, association studies and many other fields has become routine. In order to accommodate different market requirements in terms of genotyping cost, manufacturers of single nucleotide polymorphism (SNP) arrays, private companies and international consortia have developed a large number of arrays with different content and different SNP density. The number of currently available SNP arrays differs among species: ranging from one for goats to more than ten for cattle, and the number of arrays available is increasing rapidly. However, there is limited or no effort to standardize and integrate array- specific (e.g. SNP IDs, allele coding) and species-specific (i.e. past and current assemblies) SNP information. RESULTS: Here we present SNPchiMp v.3, a solution to these issues for the six major livestock species (cow, pig, horse, sheep, goat and chicken). Original data was collected directly from SNP array producers and specific international genome consortia, and stored in a MySQL database. The database was then linked to an open-access web tool and to public databases. SNPchiMp v.3 ensures fast access to the database (retrieving within/across SNP array data) and the possibility of annotating SNP array data in a user-friendly fashion. CONCLUSIONS: This platform allows easy integration and standardization, and it is aimed at both industry and research. It also enables users to easily link the information available from the array producer with data in public databases, without the need of additional bioinformatics tools or pipelines. In recognition of the open-access use of Ensembl resources, SNPchiMp v.3 was officially credited as an Ensembl E!mpowered tool. Availability at http://bioinformatics.tecnoparco.org/SNPchimp.

The assessment of inter-individual variation of whole-genome DNA sequence in 32 cows.

Despite the growing number of sequenced bovine genomes, the knowledge of the population-wide variation of sequences remains limited. In many studies, statistical methodology was not applied in order to relate findings in the sequenced samples to a population-wide level. Our goal was to assess the population-wide variation in DNA sequence based on whole-genome sequences of 32 Holstein-Friesian cows. The number of SNPs significantly varied across individuals. The number of identified SNPs increased with coverage, following a logarithmic curve. A total of 15,272,427 SNPs were identified, 99.16 % of them being bi-allelic. Missense SNPs were classified into three categories based on their genomic location: housekeeping genes, genes undergoing strong selection, and genes neutral to selection. The number of missense SNPs was significantly higher within genes neutral to selection than in the other two categories. The number of variants located within 3'UTR and 5'UTR regions was also significantly different across gene families. Moreover, the number of insertions and deletions differed significantly among cows varying between 261,712 and 330,103 insertions and from 271,398 to 343,649 deletions. Results not only demonstrate inter-individual variation in the number of SNPs and indels but also show that the number of missense SNPs differs across genes representing different functional backgrounds.

AffyPipe: an open-source pipeline for Affymetrix Axiom genotyping workflow.

UNLABELLED: The Affymetrix Axiom genotyping standard and 'best practice' workflow for Linux and Mac users consists of three stand-alone executable programs (Affymetrix Power Tools) and an R package (SNPolisher). Currently, SNP analysis has to be performed in a step-by-step procedure. Manual intervention and/or programming skills by the user is required at each intermediate point, as Affymetrix Power Tools programs do not produce input files for the program next-in-line. An additional problem is that the output format of genotypes is not compatible with most analysis software currently available. AffyPipe solves all the above problems, by automating both standard and 'best practice' workflows for any species genotyped with the Axiom technology. AffyPipe does not require programming skills and performs all the steps necessary to obtain a final genotype file. Furthermore, users can directly edit SNP probes and export genotypes in PLINK format. AVAILABILITY AND IMPLEMENTATION: https://github.com/nicolazzie/AffyPipe.git.

Relative extended haplotype homozygosity signals across breeds reveal dairy and beef specific signatures of selection.

BACKGROUND: A number of methods are available to scan a genome for selection signatures by evaluating patterns of diversity within and between breeds. Among these, "extended haplotype homozygosity" (EHH) is a reliable approach to detect genome regions under recent selective pressure. The objective of this study was to use this approach to identify regions that are under recent positive selection and shared by the most representative Italian dairy and beef cattle breeds. RESULTS: A total of 3220 animals from Italian Holstein (2179), Italian Brown (775), Simmental (493), Marchigiana (485) and Piedmontese (379) breeds were genotyped with the Illumina BovineSNP50 BeadChip v.1. After standard quality control procedures, genotypes were phased and core haplotypes were identified. The decay of linkage disequilibrium (LD) for each core haplotype was assessed by measuring the EHH. Since accurate estimates of local recombination rates were not available, relative EHH (rEHH) was calculated for each core haplotype. Genomic regions that carry frequent core haplotypes and with significant rEHH values were considered as candidates for recent positive selection. Candidate regions were aligned across to identify signals shared by dairy or beef cattle breeds. Overall, 82 and 87 common regions were detected among dairy and beef cattle breeds, respectively. Bioinformatic analysis identified 244 and 232 genes in these common genomic regions. Gene annotation and pathway analysis showed that these genes are involved in molecular functions that are biologically related to milk or meat production. CONCLUSIONS: Our results suggest that a multi-breed approach can lead to the identification of genomic signatures in breeds of cattle that are selected for the same production goal and thus to the localisation of genomic regions of interest in dairy and beef production.

Searching new signals for production traits through gene-based association analysis in three Italian cattle breeds.

Genome-wide association studies (GWAS) have been widely applied to disentangle the genetic basis of complex traits. In cattle breeds, classical GWAS approaches with medium-density marker panels are far from conclusive, especially for complex traits. This is due to the intrinsic limitations of GWAS and the assumptions that are made to step from the association signals to the functional variations. Here, we applied a gene-based strategy to prioritize genotype-phenotype associations found for milk production and quality traits with classical approaches in three Italian dairy cattle breeds with different sample sizes (Italian Brown n = 745; Italian Holstein n = 2058; Italian Simmental n = 477). Although classical regression on single markers revealed only a single genome-wide significant genotype-phenotype association, for Italian Holstein, the gene-based approach identified specific genes in each breed that are associated with milk physiology and mammary gland development. As no standard method has yet been established to step from variation to functional units (i.e., genes), the strategy proposed here may contribute to revealing new genes that play significant roles in complex traits, such as those investigated here, amplifying low association signals using a gene-centric approach.

SNPchiMp: a database to disentangle the SNPchip jungle in bovine livestock.

BACKGROUND: Currently, six commercial whole-genome SNP chips are available for cattle genotyping, produced by two different genotyping platforms. Technical issues need to be addressed to combine data that originates from the different platforms, or different versions of the same array generated by the manufacturer. For example: i) genome coordinates for SNPs may refer to different genome assemblies; ii) reference genome sequences are updated over time changing the positions, or even removing sequences which contain SNPs; iii) not all commercial SNP ID's are searchable within public databases; iv) SNPs can be coded using different formats and referencing different strands (e.g. A/B or A/C/T/G alleles, referencing forward/reverse, top/bottom or plus/minus strand); v) Due to new information being discovered, higher density chips do not necessarily include all the SNPs present in the lower density chips; and, vi) SNP IDs may not be consistent across chips and platforms. Most researchers and breed associations manage SNP data in real-time and thus require tools to standardise data in a user-friendly manner. DESCRIPTION: Here we present SNPchiMp, a MySQL database linked to an open access web-based interface. Features of this interface include, but are not limited to, the following functions: 1) referencing the SNP mapping information to the latest genome assembly, 2) extraction of information contained in dbSNP for SNPs present in all commercially available bovine chips, and 3) identification of SNPs in common between two or more bovine chips (e.g. for SNP imputation from lower to higher density). In addition, SNPchiMp can retrieve this information on subsets of SNPs, accessing such data either via physical position on a supported assembly, or by a list of SNP IDs, rs or ss identifiers. CONCLUSIONS: This tool combines many different sources of information, that otherwise are time consuming to obtain and difficult to integrate. The SNPchiMp not only provides the information in a user-friendly format, but also enables researchers to perform a large number of operations with a few clicks of the mouse. This significantly reduces the time needed to execute the large number of operations required to manage SNP data.

Signatures of selection in five Italian cattle breeds detected by a 54K SNP panel.

In this study we used a medium density panel of SNP markers to perform population genetic analysis in five Italian cattle breeds. The BovineSNP50 BeadChip was used to genotype a total of 2,935 bulls of Piedmontese, Marchigiana, Italian Holstein, Italian Brown and Italian Pezzata Rossa breeds. To determine a genome-wide pattern of positive selection we mapped the F st values against genome location. The highest F st peaks were obtained on BTA6 and BTA13 where some candidate genes are located. We identified selection signatures peculiar of each breed which suggest selection for genes involved in milk or meat traits. The genetic structure was investigated by using a multidimensional scaling of the genetic distance matrix and a Bayesian approach implemented in the STRUCTURE software. The genotyping data showed a clear partitioning of the cattle genetic diversity into distinct breeds if a number of clusters equal to the number of populations were given. Assuming a lower number of clusters beef breeds group together. Both methods showed all five breeds separated in well defined clusters and the Bayesian approach assigned individuals to the breed of origin. The work is of interest not only because it enriches the knowledge on the process of evolution but also because the results generated could have implications for selective breeding programs.

An alternative interpretation of residual feed intake by phenotypic recursive relationships in dairy cattle.

There has been increasing interest in residual feed intake (RFI) as a measure of net feed efficiency in dairy cattle. Residual feed intake phenotypes are obtained as residuals from linear regression encompassing relevant factors (i.e., energy sinks) to account for body tissue mobilization. By rearranging the single-trait linear regression, we showed a causal RFI interpretation underlying the linear regression for RFI. It postulates recursive effects in energy allocation from energy sinks on dry matter intake, but the feedback or simultaneous effects are nonexistent. A Bayesian recursive structural equation model was proposed for directly predicting RFI and energy sinks and estimating relevant genetic parameters simultaneously. A simplified Markov chain Monte Carlo algorithm was described. The recursive model is asymptotically equivalent to one-step linear regression for RFI, yet extends the analytical capacity to multiple-trait analysis.

Domestication of cattle: Two or three events?

Cattle have been invaluable for the transition of human society from nomadic hunter-gatherers to sedentary farming communities throughout much of Europe, Asia and Africa since the earliest domestication of cattle more than 10,000 years ago. Although current understanding of relationships among ancestral populations remains limited, domestication of cattle is thought to have occurred on two or three occasions, giving rise to the taurine (Bos taurus) and indicine (Bos indicus) species that share the aurochs (Bos primigenius) as common ancestor ~250,000 years ago. Indicine and taurine cattle were domesticated in the Indus Valley and Fertile Crescent, respectively; however, an additional domestication event for taurine in the Western Desert of Egypt has also been proposed. We analysed medium density Illumina Bovine SNP array (~54,000 loci) data across 3,196 individuals, representing 180 taurine and indicine populations to investigate population structure within and between populations, and domestication and demographic dynamics using approximate Bayesian computation (ABC). Comparative analyses between scenarios modelling two and three domestication events consistently favour a model with only two episodes and suggest that the additional genetic variation component usually detected in African taurine cattle may be explained by hybridization with local aurochs in Africa after the domestication of taurine cattle in the Fertile Crescent. African indicine cattle exhibit high levels of shared genetic variation with Asian indicine cattle due to their recent divergence and with African taurine cattle through relatively recent gene flow. Scenarios with unidirectional or bidirectional migratory events between European taurine and Asian indicine cattle are also plausible, although further studies are needed to disentangle the complex human-mediated dispersion patterns of domestic cattle. This study therefore helps to clarify the effect of past demographic history on the genetic variation of modern cattle, providing a basis for further analyses exploring alternative migratory routes for early domestic populations.

New Insights on Water Buffalo Genomic Diversity and Post-Domestication Migration Routes From Medium Density SNP Chip Data.

The domestic water buffalo is native to the Asian continent but through historical migrations and recent importations, nowadays has a worldwide distribution. The two types of water buffalo, i.e., river and swamp, display distinct morphological and behavioral traits, different karyotypes and also have different purposes and geographical distributions. River buffaloes from Pakistan, Iran, Turkey, Egypt, Romania, Bulgaria, Italy, Mozambique, Brazil and Colombia, and swamp buffaloes from China, Thailand, Philippines, Indonesia and Brazil were genotyped with a species-specific medium-density 90K SNP panel. We estimated the levels of molecular diversity and described population structure, which revealed historical relationships between populations and migration events. Three distinct gene pools were identified in pure river as well as in pure swamp buffalo populations. Genomic admixture was seen in the Philippines and in Brazil, resulting from importations of animals for breed improvement. Our results were largely consistent with previous archeological, historical and molecular-based evidence for two independent domestication events for river- and swamp-type buffaloes, which occurred in the Indo-Pakistani region and close to the China/Indochina border, respectively. Based on a geographical analysis of the distribution of diversity, our evidence also indicated that the water buffalo spread out of the domestication centers followed two major divergent migration directions: river buffaloes migrated west from the Indian sub-continent while swamp buffaloes migrated from northern Indochina via an east-south-eastern route. These data suggest that the current distribution of water buffalo diversity has been shaped by the combined effects of multiple migration events occurred at different stages of the post-domestication history of the species.