Whole-genome resequencing of Hanwoo (Korean cattle) and insight into regions of homozygosity

Background Hanwoo (Korean cattle), which originated from natural crossbreeding between taurine and zebu cattle, migrated to the Korean peninsula through North China. Hanwoo were raised as draft animals until the 1970s without the introduction of foreign germplasm. Since 1979, Hanwoo has been bred as beef cattle. Genetic variation was analyzed by whole-genome deep resequencing of a Hanwoo bull. The Hanwoo genome was compared to that of two other breeds, Black Angus and Holstein, and genes within regions of homozygosity were investigated to elucidate the genetic and genomic characteristics of Hanwoo. Results The Hanwoo bull genome was sequenced to 45.6-fold coverage using the ABI SOLiD system. In total, 4.7 million single-nucleotide polymorphisms and 0.4 million small indels were identified by comparison with the Btau4.0 reference assembly. Of the total number of SNPs and indels, 58% and 87%, respectively, were novel. The overall genotype concordance between the SNPs and BovineSNP50 BeadChip data was 96.4%. Of 1.6 million genetic differences in Hanwoo, approximately 25,000 non-synonymous SNPs, splice-site variants, and coding indels (NS/SS/Is) were detected in 8,360 genes. Among 1,045 genes containing reliable specific NS/SS/Is in Hanwoo, 109 genes contained more than one novel damaging NS/SS/I. Of the genes containing NS/SS/Is, 610 genes were assigned as trait-associated genes. Moreover, 16, 78, and 51 regions of homozygosity (ROHs) were detected in Hanwoo, Black Angus, and Holstein, respectively. ‘Regulation of actin filament length’ was revealed as a significant gene ontology term and 25 trait-associated genes for meat quality and disease resistance were found in 753 genes that resided in the ROHs of Hanwoo. In Hanwoo, 43 genes were located in common ROHs between whole-genome resequencing and SNP chips in BTA2, 10, and 13 coincided with quantitative trait loci for meat fat traits. In addition, the common ROHs in BTA2 and 16 were in agreement between Hanwoo and Black Angus. Conclusions We identified 4.7 million SNPs and 0.4 million small indels by whole-genome resequencing of a Hanwoo bull. Approximately 25,000 non-synonymous SNPs, splice-site variants, and coding indels (NS/SS/Is) were detected in 8,360 genes. Additionally, we found 25 trait-associated genes for meat quality and disease resistance among 753 genes that resided in the ROHs of Hanwoo. These findings will provide useful genomic information for identifying genes or casual mutations associated with economically important traits in cattle.


Background
The bovine genome was one of the first mammalian genomes sequenced, likely because cattle are important farm animals serving as major nutritional sources for humans and because of their evolutionary position as a representative of the Ruminantia, a phylogenetically distant clade to humans and rodents [1]. The bovine genome sequencing consortium sequenced a single inbred female Hereford cow and her sire using a combination of hierarchical sequencing and whole-genome shotgun sequencing [2]; the data have been assembled into two reference genomes, Btau and UMD [3,4].
After the bovine reference genome was assembled, several bovine genomes were resequenced, providing more insight into the genetic diversity of cattle that may be associated with phenotypic differences between breeds. In 2008, Van Tassell et al. reported more than 60,000 putative single-nucleotide polymorphisms (SNPs) identified from a reduced representation DNA library of 66 cattle representing three populations [5]. In 2009, Eck et al. performed the first single cattle whole-genome resequencing and reported more than 2 million novel SNPs in a Fleckvieh bull [6]. In 2011, Kawahara-Miki et al. resequenced the genome of a single Kuchinoshima-Ushi bull, a Japanese native cattle breed whose lineage has been strictly maintained in a small island secluded from mainland Japan [7]. In that study, more than 5.5 million novel SNPs were reported, and the Kuchinoshima-Ushi bull was determined to be genetically distinct from European domestic cattle breeds. Most recently, Stothard and colleagues reported whole-genome resequencing of Black Angus and Holstein, representative beef and dairy breeds, respectively, in North America, leading to the identification of substantial numbers of SNPs and copy number variants (CNVs) that could potentially be used as genetic markers across the genome [8].
When high-density genome-wide SNP data are available, analyses can identify genetic differences between similar populations. Understanding the genetic mechanisms leading to phenotypic differentiation requires identification of the genomic regions that have been under artificial selection in cattle breeds. For example, strong artificial selection will increase the frequency of favorable alleles at loci affecting meat quality traits in meatproducing breeds such as Hanwoo or Black Angus. In this process, a small region of the genome surrounding the mutations is also selected, resulting in a small genome region that shows reduced variation. Many methods have been developed for the detection of selection signatures from genome analyses, such as the use of regions of homozygosity (ROHs) [9], the integrated haplotype score (iHS) [10], F ST [11], and the extended haplotype homozygosity (EHH) statistic [12], according to the detection of the timescale for selection signatures. ROHs are without heterozygosity in the diploid states and provide association evidence at the genome-wide scale for complex traits.
Hanwoo, a Korean cattle breed, is reported to have originated from crossbreeding between taurine and zebu cattle and migrated to the Korean peninsula through North China; their history as a draft animal dates back at least 5,000 years [13,14]. Afterward, Hanwoo was maintained without the introduction of additional germplasm. Hanwoo was raised as a draft animal until the 1970s. In the late 1970s, the Korean government initiated a Hanwoo genetic breeding program to improve meat quantity and quality.
In this study, we sequenced the genome of a Hanwoo breeding bull and identified single nucleotide polymorphisms (SNPs) based on the Bos taurus reference genome assembly (Btau4.0). SNPs of Hanwoo were compared with those of Black Angus and Holstein. Moreover, functional annotation was carried out for SNPs. We also investigated genomic regions of homozygosity in Hanwoo, Black Angus, and Holstein.

Results and discussion
Genome sequencing, SNP/indel detection, and genotype concordance Whole-genome sequencing of a Hanwoo bull was performed using the ABI SOLiD platform. Approximately 6.04 billion reads were produced from three independently prepared libraries. Using BFAST,~3.68 billion reads of 120 gigabases (Gb) were aligned to Btau4.0 and filtered for redundant sequence reads. In total, 98.3% of the reference genome sequence was covered with an average mapping depth of~45.6-fold (Table 1, Additional files 1, 2 and 3). This up-to-date Hanwoo sequence coverage was the highest in the bovine genomes sequenced until now, which could facilitate more reliable SNP identification [6][7][8]. Sequencing data from Black Angus and Holstein were reanalyzed with modified parameters to compare the sequencing data of Hanwoo to Black Angus and Holstein from a previous report [8]. The mapping depth of coverage in the Black Angus and Holstein were 9.8-fold and 10.8-fold, respectively, slightly lower than that in a previous report [8]. This inconsistency may be due to a difference in the application programs and algorithms used for analysis. However, in spite of relatively low read depths of Black Angus and Holstein bulls, 97.4% and 97.7% of the reference genome was covered by the sequenced reads at the minimum read depth of 1, respectively (Additional file 3), higher than the 93% coverage with 15.8-fold mapping depth reported in Kuchinoshima-Ushi [7].
In total, 4,781,758 SNPs were identified in the Hanwoo genome using the Genome Analysis Tool Kit (GATK) 1.0.5974 [15,16] To evaluate the SNP calling from our high-throughput genome sequencing data, concordance analysis was performed between Hanwoo genome resequencing and the SNP chip data. The same genomic DNA from Hanwoo used for deep resequencing was genotyped for 54,001 SNPs using BovineSNP50 BeadChip (Illumina). All probe sequences were mapped against the Btau4.0 reference genome assembly, and 50,411 positions were identified as unique genomic loci. In total, 1,061 (2.8%) of 38,049 homozygous calls by the SNP chip have been identified as heterozygous by NGS. In total, 526 (4.3%) of 12,362 heterozygous calls by the SNP chip were identified as homozygous by NGS (Additional file 4). The overall genotype concordance was 96.2%. The nonreference sensitivity and non-reference discrepancy rates were 97.1% and 7.0%, respectively. Non-reference sensitivity is the fraction of sites called variants (A/B or B/B) in comparison to those that are also called variants in evaluation data. The non-reference discrepancy rate, which is a good measure for testing the accuracy of genotype calls, can show the accuracy of genotype calling at sites called by both sites by excluding concordant genotypes (http://gatkforums.broadinstitute.org/discussion/48/using-varianteval).

Functional annotation of genomic variation
The SNPs in genic regions were annotated using 20,955 genes from the NCBI Reference Sequence Database (RefSeq). In total, 1,663,599 SNPs (34.8%) identified in the Hanwoo genome were located in genic regions: 1,591,380 SNPs were located in introns, 21,507 SNPs were located in untranslated regions (UTRs), and 460 SNPs were located in splice sites. In total, 47,823 coding SNPs including 22,752 non-synonymous nucleotide substitutions such as missense and nonsense/read-through SNPs were also found ( Figure 1 and Additional file 3). In total, 142,297 indels (36.4%) were in genic regions, of which 2,163 indels were identified as variations that may change amino acid sequences such as frameshift, nonsense, and splice-site SNPs, which may have the potential to cause functional differences. Non-synonymous SNPs, splice-site variants, and coding indels within a coding DNA sequence (NS/SS/I), which may affect gene function, were detected in Hanwoo (24,915 in 8,360 genes), Black Angus (15,107 in 6,563 genes), and Holstein (16,963 in 6,692), respectively (Additional files 3 and 5). The Hanwoo genome contained more NS/SS/Is than those of Black Angus and Holstein. This suggests that Hanwoo is a more genetically distant breed than Black Angus and Holstein based on the reference genome of Hereford, which is consistent with a previous report [17]. Of all reference genes (20,955), 10,906 genes contained NS/SS/I genes and 737 genes revealed more than 10 NS/ SS/Is in all breeds (Additional file 5). ATP-binding cassette subfamily C member 4 (ABCC4) and zinc-finger protein 280B (ZNF280B) genes showed more than 100 NS/SS/Is. Four isoforms (copies) of the ABCC4 gene are located on BTA12 in tandem with each other (ENSBTAG0000003 2603, ENSBTAG00000047764, ENSBTAG00000023309, and ENSBTAG00000047383). Fifty-four variations (NS/ SS/Is) in four isoforms are recorded in Ensembl. However, the ZNF280B gene is a single-copy gene (ENSBTAG 00000001005) located on BTA17 and 83 NS/SS/Is exist in Ensembl, although ZNF280B has a smaller genome span (8.463 kb) and transcript (1.980 kb) compared to the genome spans (87.521 to 165.199 kb) and transcripts (2.529 to 3.930 kb) of ABCC4 gene copies. These findings show that these two genes surely belong to the gene group of more NS/SS/Is rather than other common genes. A study has reported that the number of copies of the ABCC4 gene increases and the gene is overexpressed in the process of selection for resistant mouse cells against antibiotics such as ciprofloxacin [18]. Therefore, this suggests that genes containing several NS/SS/Is may have evolved into multi-copy genes for environmental adaptation, or that NS/SS/Is may be distorted due to an incorrect reference genome sequence. However, this is necessary for experimental validation based on phenomena such as CNV or segmental duplication. Alternatively, the possibility of the presence of pseudogenes should not be excluded for genes containing several NS/SS/Is. Among 10,906 genes containing NS/SS/Is, the number of genes containing specific NS/SS/Is was 1,983 in Hanwoo, 1,199 in Black Angus, and 900 in Holstein. In Hanwoo, 1,045 genes contained reliable specific NS/SS/Is with more than tenfold depth. Furthermore, of 1,045 genes containing specific NS/SS/Is, 293 genes were revealed in Hanwoo only and 109 genes contained more than one novel damaging NS/SS/I in the functions among them (Additional file 6). Seven NS/SS/Is and six novel damaging NS/SS/Is were found in Hanwoo specifically within the raftlin lipid raft linking protein 1 (RFTN1) gene, which is important in the formation or maintenance of membrane lipid rafts [19] and is overexpressed in smooth muscles (Gene Expression Atlas in EBI). Next, we investigated whether NS/SS/I-containing genes were associated with economic traits and then categorized them into meat, disease resistance, growth, milk, and fecundity. We used previously reported information on trait-associated genes [7,20,21]. In total, 619 genes were assigned as trait-associated genes: 464 genes for meat quality, 144 genes for disease resistance, 25 genes for milk production, 8 genes for fecundity, and 6 genes for growth rate (Additional file 5). Of the 464 genes for meat quality, 228 contained more than one NS/SS/I. The titin (TTN) gene has 62 NS/SS/Is, the largest number among the genes related to meat quality. The bovine major histocompatibility complex (MHC) class I heavy chain isoform 1 precursor (BOLA, ENSB TAG00000002069) gene, which contain 32 NS/SS/Is, has the largest number of NS/SS/Is among the genes related to disease resistance (Additional file 5). The higher number of NS/SS/Is in TTN than BOLA may be due to the difference of gene size; 274.866 kb and 3.788 kb of TTN and BOLA genes, respectively. The TTN gene encodes the titin protein, the largest protein, which consists of 317 exons in 274.866 kb of genomic DNA in BTA2 (Ensembl database UMD3.1). The TTN gene plays a role in myofibrillogenesis and is associated with marbling [22]. Moreover, the 231054C>T variant within the promoter region of TTN is associated with a marbling trait and is differentially expressed between high-and lowmarbling muscle samples [23]. However, many NS/SS/Is likely affect the function of titin, which acts as a molecular spring for the passive elasticity of muscles [24]. These NS/ SS/Is within the TTN gene may be informative variants for understanding the effects of steric changes in the TTN protein.

S N P I N D E L
Of the 144 genes for disease resistance, 74 also contained NS/SS/Is, and many novel damaging NS/SS/Is were detected in several genes including BCL2-like 1 (BCL2L1), nitric oxide synthase 1 (NOS1), nucleotidebinding oligomerization domain-containing protein 2 (NOD2), granzyme A (GZMA), and semaphorin-5A (SEMA5A), as well as the BOLA gene (Additional file 5). Among 109 genes containing more than one novel damaging specific NS/SS/Is in Hanwoo, the BCL2L1, GZMA, and CD5 genes are known as candidate genes for the disease resistance trait (Additional file 6). We suggest that the exonic variation identified in this study will provide valuable information for functional studies as well as marker development associated with economic traits in cattle.

Regions of homozygosity within the three breeds
A ROH is a continuous or uninterrupted stretch of DNA without heterozygosity in the diploid state. A discrepancy has existed in the minimum standard of definition of ROH among the groups that have been studied for ROHs to date [25]. Most previous ROH studies have been performed with SNP chip results, an average of 50 SNPs of 5 Mb in size with an average distance of 100 kb between them, and an allowance of up to 2% heterozygous SNPs within a ROH [25]. However, at present, no standardized criteria have been established for defining ROHs [25]. In this study, using mass genotype data derived from whole-genome resequencing, we shortened the detection window of ROHs and loosened the permissible ratio of heterozygous SNPs (Additional file 7). Our criteria were as follows: the ROH detection window was 400 kb and 20% of heterozygous SNPs were allowed for Hanwoo, Black Angus, and Holstein ( Figure 2). We defined 16, 78, and 51 ROHs in Hanwoo, Black Angus, and Holstein, respectively (Table 2 and Additional file 8).
Angus and Holstein were bred for meat and milk production, respectively. In contrast, Hanwoo was raised as a draft animal until the 1970s. Since 1979, Hanwoo has been bred as beef cattle according to the Hanwoo genetic improvement national program organized by the government. Here, we suggest that the total lengths of ROHs in Holstein and Black Angus are longer than those of Hanwoo because Holstein and Black Angus have been artificially selected for a longer period of time. Overall, the dispersing pattern of ROHs in chromosomes was variable and also differed in the overlapping pattern of ROHs between breeds; we found two overlapping regions between Hanwoo and Black Angus, three overlapping regions between Hanwoo and Holstein, and 14 overlapping regions between Black Angus and Holstein (Additional file 8). These patterns would result from different origins and breeding strategies among the three breeds because Black Angus and Holstein originated in Aberdeen, Scotland and the Netherlands, respectively, and have been bred as beef and dairy cattle, respectively, while Hanwoo was bred independently as beef cattle in the Korean peninsula since 1979.
A total of 753 genes resided in the ROHs of Hanwoo, whereas 1,320 and 2,482 genes existed in the ROHs of Black Angus and Holstein, respectively. Among them, 77 and 30 common genes were located within the overlapping ROHs between Hanwoo and Angus and between Hanwoo and Holstein, respectively (Additional file 9). Among the 753 genes in the ROHs of Hanwoo, 505 genes contained no NS/SS/Is. A total of 2,158 genes (10.2%) contained no NS/SS/Is in the ROHs in any of the three breeds (Additional file 5). Moreover, we performed functional enrichment analysis using gene ontology (GO) for genes in the ROHs of the three breeds (Additional file 10). In Hanwoo, one significant GO term was 'regulation of actin filament length related to muscle metabolism' (GO:0030832, p-value = 0.044), including actin-related protein 3 homolog (ACTR3), actin-related protein 2/3 complex, subunit 2 (ARPC2), villin 1 (VIL1), and destrin (DSTN) genes. Meat tenderness is generated by the disruption of actin filaments and by breaking down the interaction between the actin and myosin filaments [26]. Notably, a significant GO term of 'striated muscle cell differentiation' (GO:0051146, p-value = 0.034) was found in Holstein, including the retinoid X receptor, alpha (RXRA) gene, which inhibits adipogenesis [27] and plays a negative role in marbling in Hanwoo [28]. Because Hanwoo and Black Angus were bred as beef cattle, 77 genes in overlapping ROHs between the breeds were used to analyze GO and the KEGG pathway. Although eight significant GO terms were detected, most were related to the immune system, such as T cell activation and lymphocyte activation, rather than meat traits. The presence of many immune system-related genes in the identified ROHs could  reflect selection (natural or artificial) for disease resistance. According to functional enrichment analysis using KEGG pathway terms, vitamin B6 metabolism (bta00750, p-value = 0.025) was significantly enriched, including the aldehyde oxidase 1 (AOX1) gene in Hanwoo (Additional file 10). Vitamin B6 induces the differentiation of adipocytes from pre-adipocytes and facilitates fat accumulation [29]. In particular, the AOX1 gene is a target of peroxisome proliferator-activated receptors alpha and gamma (PPARα and PPARγ) as a key gene in adipogenesis [30]. Melanogenesis (bta04916, p-value = 0.009) was also detected in Holstein. Sixteen genes, including tyrosinase (TYR) and the melanocortin 1 receptor (MC1R), exist in this term. TYR is the rate-limiting enzyme in the melanogenesis pathway. Tyrosinase activity is regulated by the MC1R. Recently, Kuhn and Weikard reported the dilution of black pigment (eumelanin) in an F 2 Holstein × Charolais population [31]. These genes may also be partially responsible for the coat color pattern of Holsteins. These observations suggest that genes within the ROH accumulate the biological functions for characteristics of each breed during a process of artificial selection. We found trait-associated genes containing NS/SS/Is in the ROHs of each breed (Table 3). Twenty-five traitassociated genes for meat quality and disease resistance were found in the ROHs of Hanwoo. Among these, nine were associated with meat quality traits and contained NS/SS/Is: TTN, acetylcholine receptor subunit alpha (CHRNA1), isocitrate dehydrogenase 1 (IDH), amylotrophic lateral sclerosis 2 (ALS2), Sp1 transcription factor (SP1), retinoic acid receptor gamma (RARG), collagen type IX alpha 3 (COL9A3), fatty acid binding protein 4 (FABP4), and insulin-like growth factor 1 receptor (IGF1R; Table 3 and Additional file 5). In humans, the homozygosity association approach has been applied to identify causal mutations for autosomal recessive disorders in consanguineous families [32][33][34][35][36][37][38][39][40], as well as the genome-wide investigation of candidate genes for complex phenotypes such as schizophrenia [41], late-onset Alzheimer's disease [42], and height [43]. Unlike humans, farm animals are artificially selected to improve genetic performance for economic traits such as meat quality and milk production. Therefore, genes related to economically important traits are gradually fixed with a dominant allele as a result of artificial selection. In livestock, fixed genes should be considered major genes that control a certain trait for any breed because genes can become fixed, dominant alleles through artificial selection or evolution. However, we suggest that genes containing NS/SS/Is in ROHs are still strong candidate genes for meat quality traits in the Hanwoo population because they may be in the selection process for breeding. According to previous reports, a splice-site mutation within the FABP4 gene in Australian cattle is significantly associated with intramuscular fat content [44]; SNPs within the FABP4 gene are associated with palmitoleic acid and linoleic acid content in intramuscular fat in Japanese cattle [45] and with backfat thickness [46], marbling score, and carcass weight in Hanwoo [47]. The IDH1 gene, which is responsible for ketoglutarate, CO 2 , and NADPH production from isocitrate in the cytosol and associated with body weight and fat deposition [48], had one damaging NS/SS/I in Hanwoo.
In addition, to detect common ROHs between the genome sequence and SNP chip, we calculated the ROHs from 40 Hanwoo bulls as well as 20 Angus and 19 Holstein individuals using the Bovine 50K SNP chip ( Figure 2). We identified four, eight, and eight common ROHs between both data sets in Hanwoo, Black Angus, and Holstein, respectively (Additional file 9). In Hanwoo, 43 genes located in common ROHs were shared between genome sequencing and SNP chips in BTA2, 10, and 13 ( Figure 2). Of 43 genes, 22 genes contained no NS/SS/I (Additional file 5). Moreover, four common ROHs in Hanwoo coincided with quantitative trait loci (QTLs) for meat fat traits ( Figure 2). Specifically, two regions in BTA2 (95.3-96.4 Mb and 100.9-101.4 Mb in Btau4.0) were common ROHs between genome sequencing and SNP chips in Hanwoo and Black Angus. Of the 18 genes that resided in these regions, WD repeat domain 12 (WDR12), amyotrophic lateral sclerosis 2 (juvenile) chromosome region, candidate 8 (ALS2CR8), cytochrome P450, family 20, subfamily A, polypeptide 1 (CYP20A1), and cAMP responsive element binding protein 1 (CREB1) genes belonged to a significant GO term of metabolic processes in Hanwoo. Among them, the CREB1 gene has been shown to be related to fat metabolism. In 2012, Lee et al. reported that the expression of the cAMP responsive element binding protein (CREB1) gene is higher in muscle with high IMF content in Hanwoo [49]. CREB1 is a transcription factor containing a basic leucine zipper. The CREB protein is phosphorylated in response to increased cAMP, allowing it to efficiently interact with the transcriptional co-activator protein, CREB binding protein, to stimulate the transcription of cAMP target genes [50]. Moreover, Casimir and Ntambi reported that intracellular cAMP activates the expression of the stearoyl-CoA desaturase gene, a key enzyme involved in monounsaturated fatty acid synthesis through activation of the CREB protein [51]. In 2009, Wang et al. observed that messenger RNA expression of a lipogenesis-related gene, stearoyl-coA desaturase (SCD), peaked at 20 to 25 months in crosses between Wagyu and Hereford, which was highly correlated with intramuscular fat content in these animals [52]. These findings suggest that elevated CREB expression may stimulate genes involved in the lipid biosynthesis pathway such as SCD [51] and HMG-Co synthase [53], resulting in an increase in IMF content within muscles. Also, the ALS2 gene, which is related to meat traits, as well as cytotoxic T-lymphocyte-associated protein 4 (CTLA4) and CD28  Figure 2). In livestock animals bred by an improvement scheme for economic traits, the use of ROHs will be a good genomic strategy for tracking and planning improvements in breeding.

Conclusions
In this study, we sequenced the whole genome of a Hanwoo bull and newly identified 2,454,142 SNPs and 342,287 small indels by comparison with the Hereford reference genome sequence. We also found 1,663,599 SNPs and 142,297 indels that were located in genic regions of 20,955 genes in the NCBI Reference Sequence Database (RefSeq), of which 22,752 SNPs and 2,163 indels were non-synonymous, frameshift, nonsense, or splice-site SNPs potentially capable of affecting protein functions. This suggests that genes containing several NS/ SS/Is may have evolved into multi-copy genes for environmental adaptation, or that NS/SS/Is may be distorted due to an incorrect reference genome sequence. A ROH is a continuous or uninterrupted stretch of DNA without heterozygosity in the diploid state. In this study, we defined 16 ROHs in Hanwoo using a detection window of 400 kb and 20% of heterozygous SNPs using genotype data derived from whole-genome resequencing. The cumulative lengths of ROHs per genome, as well as the number of ROHs in Hanwoo, were smaller than those in Black Angus and Holstein. This suggests that the total lengths of ROHs in Holstein and Black Angus are longer than those of Hanwoo due to a longer period of time for artificial selection in those breeds. In addition, the dispersing pattern of ROHs in chromosomes was different between breeds. We suggest that these patterns would result from the different origins and breeding strategies among these three breeds. Moreover, 753 genes were observed in the ROHs of Hanwoo, of which 25 genes were associated with meat quality and disease resistance traits. In addition, we observed common ROHs between the genome sequence and high-density SNP chip data. This combinatorial ROH survey approach may be another effective method for identifying domestication genes. The findings of this study will provide valuable information for functional studies, as well as for marker development associated with economically important traits in cattle.

DNA samples
We sequenced the genome of a proven Hanwoo bull (27223) obtained from the Hanwoo Experiment Station, National Institute of Animal Science, Rural Development Administration, Korea. Bull 27223 was selected for mapping for its representativeness of the population at the Hanwoo Experiment station. Bull 27223 is a descendent of KPN369, which was one of the most frequently used Hanwoo bulls for artificial insemination in Korea during the early 2010s. Also, bull 27223 was selected for its superiority in growth performance with superior genetic potential in carcass quality. Therefore, many calves born since 2010 have been sired by this bull. The study protocol and standard operating procedures were reviewed Briefly, gDNA was fragmented using Covaris S2 (Covaris) and HydroShear (Genomic Solutions) at the proper settings for targeted sizes. A QIAquick Gel Extraction Kit (Qiagen) was used for subsequent purification of sheared DNA, enzymatic reactions, and size-selected DNA in agarose gels according to the manufacturer's instructions. To repair damaged DNA ends and obtain 5′-phosphorylated blunt-ends (5′P), the fragments were end-repaired using the End-It DNA End-Repair Kit (Epicentre Biotechnologies) according to the manufacturer's instructions. Ligations for the adaptor attachment and circularization were accomplished using the Quick Ligation Kit (New England BioLabs). DNA quantitations were performed using a NanoDrop ND 1000 Spectrophotometer (Thermo Fisher Scientific), except for those followed by library amplification for emulsion PCR (ePCR).
In chronological order, the sheared gDNA fragments were end-repaired and the LMP CAP Adaptors (missing the 5′ phosphate from one oligonucleotide resulting in a nick on each strand when the DNA is circularized at a later step) were ligated to the end-repaired DNA fragments. The adaptor-ligated products were separated on a 1% agarose gel and excised from the gel at the appropriate positions for span size ranges (600-700 bp, 1-2 kb, and 0.6-2.2 kb). Size-selected DNA fragments were circularized with a biotinylated internal adaptor. Uncircularized DNA fragments were eliminated using Plasmid-Safe ATP-Dependent DNase (Epicentre Biotechnologies). Nick translation was performed for 14 min at 0°C in an icewater bath using Escherichia coli DNA polymerase I with the circularized DNA fragments. The nick-translated products were cleaved at the nicks using T7 exonuclease and S1 nuclease, and end-repaired as described above. P1 and P2 adaptors (used for library amplification, ePCR, and ligation sequencing) were ligated to ends of the endrepaired DNA. Then the ligated DNA underwent nick translation with DNA polymerase I. The completed library was amplified using Library PCR primers 1 and 2 with Cloned Pfu polymerase (Stratagene) or Platinum® PCR Amplification Mix (SOLiD Long Mate-Paired Library Construction Kit, ABI). The amplified library was ran on a 4% agarose gel and the correct-sized band (275-300 bp) was excised and eluted, and quantitated by Qubit IT (Invitrogen). ePCR was carried out according to the Applied Biosystems SOLiD System: Template Bead Preparation Guide. The concentration of each library for ePCR was designed to range from 1.0 to 1.5 pM.

Library sequencing of template beads
Sequencing was performed according to the Applied Biosystems SOLiD System: Instrument Operation Guide. Templated beads were deposited onto two slides and sequencing was carried out to 50 bases using SOLiD v3.0 chemistry, with the exception that the library prepared from 0.6-2.2 kb-sheared DNA fragments was used for four slides and sequencing was carried out to 50 bases using SOLiD v3 plus chemistry.

Short-read alignment, variant calling, and annotation
Paired-end 50 bp reads from Hanwoo, Black Angus, and Holstein were mapped to the Btau4.0 reference genome assembly using BFAST 0.7.0a [54], with options bfast match "-A 1 -z -K 100 -M 500," bfast localalign "-A 1 -o 10," and bfast postprocess "-A 1 -a 3 -Y 2 -z -O 1." Aligned reads considered to be PCR duplicates were removed using the MarkDuplicates algorithm in Picard tools 1.57. This algorithm identifies the 5′ coordinates and mapping orientations of each read pair by considering gaps and jumps. The reads that mapped to the same position and orientation are marked as duplicates except the best scored read pair. The score of a read pair is defined as the sum of base qualities >15. Next, the IndelRealigner module in the Genome Analysis Toolkit (GATK) 1.0.5974 [15] was used to perform local realignment around indels to produce an accurate alignment and CountCovariates and TableRecalculation modules to recalibrate the base quality score. An in-house script was applied to modify the read quality, which was generated by BFAST before the GATK recalibration step. The quality scale generated by BFAST presented up to~63 and was skewed to the maximum value. Such an overestimated quality scale prevented the filtration of false-positive variations while GATK runs genotyping. The in-house script scaled down the overestimated quality values to~40. SNP and small indel calling were performed using GATK UnifiedGenotyper [16] with a minimum base quality of Q17 (phred score 17) with "-stand_call_conf 0 --stand_ emit_conf 0 --max_deletion_fraction 1.00" and a minimum mapping quality of Q30 (phred score 30) with "-stand_call_conf 0 --stand_emit_conf 0 --genotype_likeli hoods_model INDEL --minIndelCnt 3". Hanwoo, Black Angus, and Holstein were genotyped separately using GATK UnifiedGenotyper. Then, the variants identified in three breeds were merged by genomic position for downstream analysis. A novel variant was defined as one that was not present in the cattle dbSNP 133.
Annotations of variants were based on the 34,577 Cow RefSeq in NCBI (downloaded April 2, 2012). The cattle RefSeqs were aligned against Btau4.0 using BLAT with the 'fine' option to obtain the genomic positions of genes, exons, and coding regions. In total, 33,080 RefSeqs were aligned against the reference genome. Among the aligned RefSeqs, the sequences with >90% coverage and a <1% error rate were selected. Then one representative RefSeq was selected from the RefSeqs derived from the same gene. As the result, we selected 29,197 RefSeqs for variant annotation. We identified 2-base canonical splice sites (GU/AG) at the end of an intron as a splice site. The genomic locations of some trait-associated genes that were not obtained from NCBI RefSeqs were defined from previously reported gene information [7]. The selected genes were used to predefine the annotation data of all possible variants and pre-calculate the SIFT [55] predictions and scores. We selected the coding indels, splice-site variants, and non-synonymous SNPs (NS/SS/Is) that showed SIFT scores of <0.05 as the potentially damaging variants.
Specific NS/SS/I variants were detected by the following criteria: We first selected the NS/SS/Is for which at least 10 reads were aligned and an allele was 50% more abundant than the other alleles for all three breeds at the position.

ROHs
To measure the genome-wide pattern of selection of a breed, we defined a ROH as follows. The minimum ROH size was set to 400 kb; each chromosome was divided into 400 kb bins, and the ratio of homozygous SNPs per bin was employed as the degree of homozygosity of the bin. To look for a series of high-degree bins rather than separated, one-point peak bins, a degree was smoothed by an average of the two neighbor bins on each side. A continuous extension of bins with a high degree of homozygosity was defined as a ROH. In this study, a 0.8 degree was imposed to determine the ROHs. One breed may contain a ROH that shows a high degree of homozygosity while the others do not. This helps to explain the breed-specific selection pressure. We defined a subset of ROHs that was not duplicated in the other breed's ROHs as specific ROHs (sROHs). ROHs were identified from SNP chip data using HomozygosityMapper [56].

SNP genotyping
To evaluate the accuracy of SNP calling from resequencing of the Hanwoo genome, the same genomic DNA sample was applied to SNP chip analysis. We used BovineSNP50 BeadChip (Illumina) [57] to genotype the Hanwoo genome. In total, 40 proven bulls in the 45th Hanwoo Performance and Progeny Test Program in Korea, as well as 20 Angus and 19 Holstein individuals, were used for SNP genotyping with the same platform to investigate the ROHs. A consensus SNP genotype was obtained by selecting a maximally expressed genotype from the same location in a breed. Over 90% of the consensus genotypes appeared in more than half of the individuals for all three breeds (data not shown). A ROH was computed from the consensus SNP genotypes with the same criteria that were applied to calculate the ROH of whole-genome SNPs.

Trait-associated genes and QTL regions
We obtained information on trait-associated genes from previous reports to analyze the Kuchinoshima-Ushi breed genome [7]. The genes were categorized into five economic traits: meat, disease resistance, growth, milk, and fecundity. Some genes that did not appear in NCBI RefSeq were added to the gene set for further analysis. QTL regions were identified from information on Cattle QTLs in the Animal QTLdb (Release 17; http://www. animalgenome.org/cgi-bin/QTLdb/BT/index) [20,21]. QTL locations by bp (UMD 3.1) were downloaded and three types of QTLs were selected: meat fat, meat tenderness, and milk traits. The associated names of these three QTL types described in QTLdb are as follows: intramuscular fat, marbling score, and marbling score (EBV) for meat fat; shear force and tenderness score for meat tenderness; and milk yield, milk yield (daughter deviation), milk yield (EBV), milk yield (PTA), dairy capacity composite index, and dairy form for milk traits.
Because the QTL were based on the UMD 3.1 genome, we converted the locations to coordinates from the Btau4.0 genome. Sequences of the selected QTL were extracted from UMD 3.1 genome sequences and aligned to the Btau4.0 reference genome sequences using LASTZ with the following options: seed = 14 of 22; chain = gapped, step = 5. The alignments were filtered with a minimum of 1,000 bases, 99% average identity, and 5% coverage. The syntenic locations were merged into a large location allowing gaps of 10% at the syntenic locations at most.

Functional enrichment analysis
We determined genes whose genomic positions overlapped partially or completely with the ROH for each breed. We performed functional enrichment analysis against the candidate genes that were within a ROH region within the Gene Ontology and KEGG pathway terms using the Database for Annotation Visualization and Integrated Discovery (DAVID) tool (http://david. abcc.ncifcrf.gov/). Only the enriched GO terms with raw p-values <0.05 were used for further interpretation in this study. The functional relationships of the genes of interest were used in the Pathway studio program (Stratagene) [58].