Genetic structure and relationships of 16 Asian and European cattle populations using DigiTag2 assay

Abstract In this study, we genotyped 117 autosomal single nucleotide polymorphisms using a DigiTag2 assay to assess the genetic diversity, structure and relationships of 16 Eurasian cattle populations, including nine cattle breeds and seven native cattle. Phylogenetic and principal component analyses showed that Bos taurus and Bos indicus populations were clearly distinguished, whereas Japanese Shorthorn and Japanese Polled clustered with European populations. Furthermore, STRUCTURE analysis demonstrated the distinct separation between Bos taurus and Bos indicus (K=2), and between European and Asian populations (K=3). In addition, Japanese Holstein exhibited an admixture pattern with Asian and European cattle (K=3‐5). Mongolian (K=13‐16) and Japanese Black (K=14‐16) populations exhibited admixture patterns with different ancestries. Bos indicus populations exhibited a uniform genetic structure at K=2‐11, thereby suggesting that there are close genetic relationships among Bos indicus populations. However, the Bhutan and Bangladesh populations formed a cluster distinct from the other Bos indicus populations at K=12‐16. In conclusion, our study could sufficiently explain the genetic construction of Asian cattle populations, including: (i) the close genetic relationships among Bos indicus populations; (ii) the genetic influences of European breeds on Japanese breeds; (iii) the genetic admixture in Japanese Holstein, Mongolian and Japanese Black cattle; and (iv) the genetic subpopulations in Southeast Asia.


INTRODUCTION
In Asia, 258 cattle breeds have been reported, and some of these breeds (11%) are classified as being at the risk of extinction (Scherf 2000). In addition to these breeds, diverse native cattle, which are not classified as breeds, are raised in the fields of Asia. Therefore, understanding the genetic diversity, structure, relationships and evolutionary history of Asian cattle are important for improving future breeding and for the conservation of genetic resources.
The bovine whole genome sequence has been determined and over 2.2 million putative single nucleotide polymorphisms (SNPs) were identified (Bovine HapMap Consortium et al. 2009). Recently, high-density SNP arrays, such as the Illumina Bovine SNP50 BeadChip (50 K) and the Illumina BovineHD BeadChip (770 K) (Illumina Inc., San Diego, CA, USA), were developed, providing powerful tools for genetic studies in cattle. However, using a smaller number of SNP markers is preferable to decrease the experimental costs in large-scale investigations. DigiTag2 assay is a SNP typing system with a high conversion rate (>90%), high accuracy and low cost (Nishida et al. 2007), which has been successfully used as a typing method based on a moderate number of SNPs (Hijikata et al. 2012;Nishida et al. 2012;Shimogiri et al. 2012;Srilohasin et al. 2014).
In this study, we assessed the genetic diversity, structure and relationships of 470 unrelated cattle from 16 cattle populations using 117 autosomal SNPs for genotyping in a DigiTag2 assay. We investigated representative samples from native cattle in Southeast Asia and Northeast Asia, four Wagyu breeds, and three European cattle to elucidate the genetic structure and relationships of Asian cattle.

SNPs selection
The DigiTag2 assay is comprised of 96 SNPs in each SNP panel. In this study, two panels were developed and total of 192 SNPs were selected from the 2641 SNPs described by McKay et al. (2008) and Illumina BovineSNP50 BeadChip (Illumina Inc., San Diego, CA, USA). We employed an interval of >8.2 Mbp between two markers on the same autosome to avoid the linked loci (Khatkar et al. 2008). Among the 192 SNPs, 117 SNPs (52 SNPs on the first panel and 65 SNPs on the second panel) had a SNP call rate ≥95% and a minor allele frequency (MAF) ≥0.05, and were therefore utilized for analysis. Further information related to these SNP markers is provided in Table S1. In addition, 10 samples with individual call rates <0.05 (one Tosa Japanese Brown, three Higo Japanese Brown, one Japanese Shorthorn, one Japanese Black, two Myanmar cattle and two Bhutanese cattle) were excluded from the present study.

Statistical analyses
The genotype and allele frequencies, the proportions of polymorphic markers and the expected (He) and observed (Ho) heterozygosities were calculated as indices of genetic diversity. Analyses of pairwise population differentiation (Fst) and within population differentiation (Fis) were performed with ARLEQUIN ver 3.5.1.2 (Excoffier et al. 2005).
To investigate the phylogenetic relationships among populations, the standard genetic distance of Nei (1972) was calculated from allele frequencies using PHYLIP 3.6 (Felsenstein 1989). Two types of phylogenetic trees were constructed from the distance matrix using the unweighted pair group method with arithmetic mean (UPGMA) (Sneath & Sokal 1973) and neighbor-joining tree (NJ) algorithms (Saitou & Nei 1987) in MEGA 5.03 (Tamura et al. 2007). Principal component analysis (PCA) was performed with MVSP 3.1 (Kovach 2005) and the principal components were calculated from all the allele frequencies in each population.
The population structure and degree of admixture were inferred using STRUCTURE 2.3.4 (Pritchard et al. 2000), including prior information on populations. STRUCTURE uses Bayesian clustering of multi-locus genotypes to assign individuals to populations, thereby estimating individual admixture proportions and inferring the number of parental populations (K) for a given sample. To obtain a representative value of K for modeling the data, we performed 16 independent runs of the Gibbs sampler for each K (1 ≤ K ≤ 17) with a 20 000 initial burn-in used to minimize the effect of the starting configurations, followed by 100 000 Markov Chain Monte Carlo iterations, as recommended by Falush et al. (2007). We used the default settings, and the admixture model with correlated allele frequencies and the parameter of individual admixture alpha were set to be the same for all clusters.

RESULTS
We estimated the allele frequencies of 117 SNP markers in 16 populations. The mean MAF in each population ranged from 0.077 (Mongol) to 0.419 (Bangladesh) ( Table 1). The mean MAF in Bos taurus populations (0.267) was higher than that in Bos indicus (0.195). We also calculated Ho and He within populations (Table 1). In the Bos taurus populations, Ho and He ranged from 0.315 (Higo Japanese Brown) to 0.394 (Mongol) and from 0.308 (Higo Japanese Brown) to 0.389 (Mongol), respectively. In contrast, those in Bos indicus populations were relatively low ranging from 0.229 (Laos) to 0.287 (Bhutan) for Ho and from 0.250 (Bangladesh) to 0.294 (Bhutan) for He. No deviations from the Hardy-Weinberg equilibrium between Ho and He were observed in all populations (P > 0.05).
We estimated Nei's genetic distance and Fst among all the cattle populations ( Table 2). The genetic distances between Bos taurus and Bos indicus were high (0.140-0.326), whereas those among Mongol cattle and Bos indicus populations were relatively low (0.140-0.186). The highest value was observed between Japanese Polled and Bangladesh (0.326) and   Table 1. the lowest value was observed between Cambodia and Vietnam or Laos (0.008). In Bos taurus, the distance between Tosa Japanese Brown and Japanese Shorthorn was highest (0.203), and that between Mongol and Hanwoo was lowest (0.042). The distances among Bos indicus populations were relatively low (0.008-0.041). Two types of phylogenetic constructions, NJ tree (Fig. 1A) and UPGMA tree (Fig. 1B), among 16 cattle populations based on Nei's genetic distance illustrated the main divergence between the Bos taurus and Bos indicus clades. The Bos indicus populations were clustered close together. The Japanese Shorthorn and Japanese Polled were clustered with the European populations.
According to the PCA (Fig. 2), Bos taurus and Bos indicus were distinguished by PC1 (17.6%). European and Asian breeds were separated by PC2 (4.5%), whereas Japanese Polled and Japanese Shorthorn were clustered with European breeds. Japanese Holstein were located in an intermediate position between the European and Asian populations. Bos indicus populations were more closely clustered than Bos taurus populations as well as results of phylogenetic constructions. Figure 3 shows the results of the STRUCTURE analysis with representative K values. We inferred the optimum number of genetic ancestral populations (K) by computing the log likelihood of the data (LnP(D)) and the second order rate of change of the likelihood function with respect to K (ΔK) (Evanno et al. 2005).
Bos taurus and Bos indicus were clearly distinguished at K=2 with the maximum ΔK value. At K=3, the separation between European and Asian Bos taurus populations was observed, although Japanese Shorthorn and Polled clustered with the European populations. At K=5, Tosa and Higo Japanese Brown were separated into two distinct groups. In addition, Japanese Holstein revealed admixture patterns with European and Asian populations at K=3-5. At K=11, each Bos taurus population formed an independent cluster. At K=2-11, the Bos indicus populations clustered into one group, but at K=12-16, the Bhutan and Bangladesh cattle populations formed a separate cluster from the other Bos indicus populations. Mongolian, Japanese Black and Korean cattle showed admixture patterns that originated from different ancestries at K=13-16, K=14-16 and K=16, respectively.

DISCUSSION
In this study, we evaluated the genetic diversity, relationships and structure of 16 cattle populations using 117 autosomal SNPs. The genetic indices (Table 2), phylogenetic constructions (Fig. 1), PCA (Fig. 2) and STRUCTURE analysis (Fig. 3) determined the distinct subdivisions of Bos taurus and Bos indicus subspecies. These results agree with the historical background, which states that the cattle originated from independent domestication events in different locations in the Fertile Crescent (Bos taurus) and on the Indian subcontinent (Bos indicus).  Table 1. Bos taurus population indicated by blue and Bos indicus by green letters. Scale bar indicates the standard genetic distance of Nei (1972).  Table 1. Bos taurus population indicated by blue and Bos indicus by green letters.