Whole Genome Resequencing of Heugu (Korean Black Cattle) for the Genome-Wide SNP Discovery

Abstract

Heugu (Korea Black Cattle) is one of the indigenous cattle breeds in Korea; however there has been severe lack of genomic studies on the breed. In this study, we report the first whole genome resequencing of Heugu at higher sequence coverage using Illumina HiSeq 2000 platform. More than 153.6 Giga base pairs sequence was obtained, of which 97% of the reads were mapped to the bovine reference sequence assembly (UMD 3.1). The number of non-redundantly mapped sequence reads corresponds to approximately 28.9-fold coverage across the genome. From these data, we identified a total of over six million single nucleotide polymorphisms (SNPs), of which 29.4% were found to be novel using the single nucleotide polymorphism database build 137. Extensive annotation was performed on all the detected SNPs, showing that most of SNPs were located in intergenic regions (70.7%), which is well corresponded with previous studies. Of the total SNPs, we identified substantial numbers of non-synonymous SNPs (13,979) in 5,999 genes, which could potentially affect meat quality traits in cattle. These results provide genome-wide SNPs that can serve as useful genetic tools and as candidates in searches for phenotype-altering DNA difference implicated with meat quality traits in cattle. The importance of this study can be further pronounced with the first whole genome sequencing of the valuable local genetic resource to be used in further genomic comparison studies with diverse cattle breeds.

Introduction

Heugu (Korean black cattle) is one of the four Korean native cattle breeds that were registered in the Food and Agricultural Organization (Choi, 2009; FAO, 2013; Jo et al., 2012). It is widely known that the population size of Heugu has been drastically diminished partly due to the policy of unifying coat color at the beginning of the 20th century. During recent decades, the depleted population size has not been restored, while the beef industry has mainly concentrated on Korean Brown Cattle as a representative beef cattle breed in Korea. The current population of the Heugu is very low and it is currently raised in very limited areas such as Chungcheong and Kangwon areas in Korea peninsula. Recently, the value of diverse genetic resources has been of great interest to conserve valuable local genetic resources for biodiversity and pioneering an additional beef market in Korea (Choi et al., 2012). As recognizing the importance of conserving the local genetic resources, National Institute of Animal Science in Republic of Korea initiated the animal genetic resource information management system to manage the phenotypic as well as genotypic data of livestock (Choi et al., 2012).

Improving meat quality in cattle can be achieved mainly by management systems and improving genetic potential. To improve the genetic potential in cattle, selective breeding based on quantitative genetic theory has been traditionally applied, leading to a substantial improvement of economically important traits in cattle. Particularly, major meat quality traits are known to have an intermediate to high heritability, suggesting that improving meat quality can be achieved by genetic selection (Dikeman et al., 2005). However, selective breeding is costly and time consuming, and is limited to traits easily measured. Improving the meat quality traits is very difficult by selective breeding because the trait can be usually measured after slaughter. To overcome the limitation of traditional breeding schemes, genomic technologies have received intensive attention in livestock research communities (Dekkers, 2004). Currently, there are several commercial DNA markers available to predict meat quality such as Igenity TenderGENETM and GeneSTAR marbling for tenderness and marbling, respectively. Most meat quality traits are complex traits that are affected by lots of genes with most of each gene contributing a small amount to the meat trait. Therefore, numerous numbers of genetic markers such as single nucleotide polymorphisms (SNPs) need to be evenly spread across the whole genome to locate each casual gene affecting the traits of interest.

The SNP is the most abundant genetic variation that accounts for genetic variations in traits of interest in cattle. It has been widely used to identify promising genomic regions affecting economically important traits in cattle by genome wide association analysis and signatures of selection studies (Barendse et al., 2009; Bolormaa et al., 2011; Qanbari et al., 2011; Sherman et al., 2010). Large-scale SNP has become publically available since the international bovine genome sequencing and HapMap projects (Elsik et al., 2009; Gibbs et al., 2009). Furthermore, the recent technological advances in massively parallel sequencing has expedited to increase SNP density from diverse cattle breeds (Canavez et al., 2012; Choi et al., 2013; Kawahara-Miki et al., 2011; Lee et al., 2013; Stothard et al., 2011). However, there is currently a drastic shortage of the genome sequencing studies on native cattle breeds in Asia, for example we could not locate any internationally peer-reviewed article for genome-wide investigations on Heugu (search terms in PubMed: Heugu and Korean Black cattle). To the best of our knowledge, this is the first article to present a result from whole genome sequencing on Heugu. The main objective of this study is to identify large scale genome-wide SNPs by whole genome sequencing at higher coverage using Illumina HiSeq 2000 sequencing platform, followed by further annotation to exhibit functional categories of all the detected SNPs.

 

Materials and Methods

DNA sampling

For whole genome sequencing (WGS) reaction, we obtained genomic DNA from a horned Heugu bull having birth at April 2th, 2011 (animal ID:002065936033), and the selected bull was accessed at Chungbuk Institute of Livestock and Veterinary Research, Cheongwon, Republic of Korea. Whole blood was initially extracted and used to isolate the genomic DNA for further sequencing reaction, according to the protocols by the Institutional Animal Care and Use Committee of the National Institute of Animal Science (Suwon, Korea). Because the quality of genomic DNA is critical to get quality sequence data for the massively parallel sequencing, we applied the careful quality checks of the DNA. The quantity and quality of the extracted DNA were examined by Nano-Drop spectrophotometer (Nano-Drop Technologies, USA) and gel electrophoresis using 1% agarose gel. Further concentration check of the double-stranded DNA (dsDNA) was followed using Qubit dsDNA HS Assay (Invitrogen, USA).

Library preparation and sequencing

Using the purified genomic DNA, libraries were constructed to have template insert sizes in 300-400 base-pair (bp) using Illumina TruSeq DNA Sample Preparation Kit according to the manufacturer's instructions (Illumina, USA). A total of 153 Giga base pairs in 101 bp paired-end short reads were generated by massively parallel sequencing reaction using the Illumina HiSeq 2000 sequencer. Out of the total initial reads, over 97% of the reads were mapped to the bovine reference genome sequence assembly, UMD 3.1 (Zimin et al., 2009). In order to call high confidence SNPs, a further filtering step was applied to remove potential duplicated reads, generating approximately fifty percent of the sequenced reads (76.9 Giga base pairs) to be used in the further SNP calling.

Read mapping and SNP genotyping

The sequenced reads were mapped to the bovine reference genome (UMD 3.1) using a fast short-read alignment program, Burrows Wheeler Aligner (BWA) version 0.6.2 with default option (Li and Durbin, 2009). Among the initially aligned reads, we further filtered out short reads which were potentially generated by PCR duplicate, using the Mark Duplicates in Picard software package version 1.59 (http://picard.sourceforge.net/). Mate information of the reads was re-synchronized using the Picard Fix Mate Information. Base-pair quality scores were recalibrated with Base-Recalibrator and Print-Reads in the Genome Analysis Tool kit (GATK) version 2.4 (McKenna et al., 2010), and then alignments around small insertion and deletions (InDels) were re-aligned using GATK Indel-Realigner.

The GATK Unified Genotyper (DePristo et al., 2011) was applied to generate an initial set of SNP. The SNPs were called under the conditions as the 30.0 of standard call confidence and SNP genotype likelihood model. In order to identify high quality variants, those initially called SNPs were evaluated with a Gaussian mixture model that was built based on 13,106,312 known bovine variants in single nucleotide polymorphism database (dbSNP) version 137 and then the resulting outliers were discarded. Furthermore, five annotations (Hardy Weinberg, Homopolymer Run, Transmission Disequilibrium Test, Allele Balance and Variant Type) were applied to each SNP call using GATK Variant Annotator. SNP filtration was performed using GATK Variant Filtration by following criteria: (1) SNP with overall quality (QUAL) less than 100, (2) quality divided by non-reference depth (QD) score less than 2.0, (3) Fisher exact test to detect strand bias (FS) more than 60.0, (4) root mean square of mapping quality (MQ) less than 40.0, (5) consistency of the site with two segregating haplotypes (Haplotype Score) more than 13.0, (6) u-based z-approximation from the Mann-Whitney Rank Sum test for mapping quality (MQ Rank Sum) less than -12.5 or (7) the rank sum test for the distance from the end of the read for reads with the alternative allele (Read Pos Rank Sum) less than -8.0. 6,043,690 SNPs which passed all above criteria were collected using GATK Select Variants.

To estimate the functional effects of SNPs in genes, we annotated and predicted the effects using SnpEff version 3.3f (Cingolani et al., 2012b). We used the UMD 3.1.71 (based on genome version UMD 3.1 and ensemble gene annotation version 71) and the bovine dbSNP version 137 to annotate each of the SNPs detected from Heugu. All the SNPs were assigned with diverse range of functional categories based on genomic coordinates, functional class, codon change, amino acid change, amino acid length, gene name, transcript biotype (if available), gene coding, transcript ID, exon rank and corresponding genotype (for further details, refer to the ‘http://snpeff.sourceforge.net/SnpEff_manual.html’). Subsequently, SnpSift (Cingolani et al., 2012a) was applied to convert the resulting variant calling format file to a plain table-formatted text file, and we further retrieved the non-synonymous SNPs (nsSNPs) from the file for the further downstream analysis. In addition, Ensembl protein ID obtained from BioMart (Ensembl version 72) was added to the text file (Kasprzyk, 2011).

Fig. 1.The Heugu bull that was whole-genome resequenced in this study, sampled at Chungbuk Institute of Livestock and Veterinary Research, Cheongwon, Republic of Korea. A picture in the left side presents the front face of the Heugu, while the right side one shows a side of the Heugu.

 

Results and Discussion

Sequencing and mapping

Genomic DNA was isolated from whole blood extracted from a selected Heugu bull, and it was used to perform massively parallel sequencing using the Illumina HiSeq 2000 sequencing platform. To our knowledge, this study is the first whole-genome sequencing (WGS) result of Heugu. The WGS reaction initially generated 1,520,994,456 sequence reads in length (totally 153.6 Giga base-pairs: Gb), of which 1,477,464,807 reads (97.1%) were mapped to the bovine reference sequence assembly (UMD 3.1) using BWA (Table 1). In order to detect reliable SNPs, we applied further steps to the reads, including removing potential PCR duplicates and realigning reads around insertion and deletion regions. After removing the duplicated reads, 761,330,287 sequence reads were obtained, corresponding to a sequence depth coverage of 28.9-fold over the bovine UMD 3.1 reference assembly. On average, the sequence reads covered approximately 98.9% of the UMD 3.1 reference genome, and they were used for further SNP calling process (Table 1). As shown in a previous study (Choi et al., 2013), the sequence depth applied in this study is high enough to call the high-confidence SNPs across the genome.

Table 1.Abbreviations in this table are: bp, base pair; x, fold.

Table 2.Abbreviations in this table are: Ti/Tv, Transition to Transversion; UTR, untranslated region.

SNP detection

A total of 6,043,690 SNPs were identified throughout the genome, of which 29.4% (1,774,965 SNPs) were found to be novel by comparing against the dbSNP version 137 (Table 2). The proportion of the novel SNPs in this study is lower than previous cattle WGS studies such as 87% and 45% for a Japanese native breed, Kuchinoshima-Ushi and a Korean native breed, Chikso respectively (Choi et al., 2013; Kawahara-Miki et al., 2011). This result is not deviated from our expectation because we used the most up-to-date dbSNP version that further includes extensive sets of SNPs particularly from the international 1000 bull genomes project. Despite the lower proportion, the status of novelty is not trivial, suggesting that further sequencing need to be performed from multiple numbers of animals and breeds for an extensive catalogue of SNPs in cattle. The heterozygous to homozygous ratio for all the detected SNPs was calculated to show ~1.63 (3,744,590:2,299,100 SNPs). Furthermore, transition to transversion ratio of the SNPs was computed as 2.21 (Table 2), which approximated the transition to transversion ratio value empirically obtained from the WGS study in Human using over a thousand individuals (Abecasis et al., 2012). Because the transition to transversion ratio can be one of the indicatives to access the putative SNP quality, the result suggest that the SNPs detected in this study have reasonable quality which can be used in further genomic investigations in cattle.

Annotation of the identified SNPs

To access potential functional roles of the putative SNPs, an extensive annotation was performed on each of the SNPs, using GATK Variant Annotator (McKenna et al., 2010). Through the annotation, a diverse range of functional class terms were assigned to each of the SNPs:downstream, exon (non-coding), intergenic, intron, exon (mitochondria), non-synonymous coding, non-synonymous start, start gained, start lost, stop gained, stop lost, synonymous coding, synonymous stop, upstream, untranslated region at three prime, and untranslated regions at five prime (Table 2). Most of the SNPs were located between genes or within intron across the genome (5,764,920:~95% of the total SNPs); the 4,275,170 and 1,489,750 SNPs were located in intergenic and intronic region respectively (Table 2). The result is well coincided with previous WGS studies in cattle (Choi et al., 2013; Kawahara-Miki et al., 2011; Lee et al., 2013; Stothard et al., 2011). Notably, a significant number of SNPs were identified in coding regions across the genome in this study. The SNPs in the coding regions were more deeply annotated to exhibit 18,849 synonymous SNPs obtained from synonymous-coding and -stop, while we located 13,979 nonsynonymous SNPs (nsSNPs) from missense and nonsense mutation that could affect traits of interest in cattle (Table 2). Another functionally important SNP would be the one located in splice sites which could affect exon configurations, and we located a total of 261 SNPs in the splice site in this study.

Functional annotation of nsSNPs and its potential functional significance

As noted, we identified numerous numbers of nsSNPs that potentially lead to genetic variations in economically important traits in cattle. The nsSNPs have an impact on protein stability and degradation by amino acid sequence change, and it is widely known that many phenotypes are implicated with the nsSNPs in other organisms such as Human (Stenson et al., 2003). In total, 13,979 nsSNPs were detected in 5,999 genes that were based on Ensembl gene ID, and several nsSNP were instantly located in the genes which have been known to be implicated with carcass and meat quality traits in cattle.

For example, we identified a SNP at 44,069,063 bp in chromosome 29 for calpain 1 gene (CAPN1) and five SNPs were located in Calpastatin (CAST) gene (98,485,261, 98,485,273, 98,535,683, 98,535,716, and 98,554,459 bp in chromosome 7). Beef tenderness is the most important palatability aspect for consumers especially in Canada and the US, and these two genes have been widely known as the most critical genes implicated with the proteolytic tenderisation process (Koohmaraie and Geesink, 2006). There are currently several reports available to show that SNPs in the CAPN1 gene or the CAST are associated with tenderness scores in cattle (Casas et al., 2006; Drinkwater et al., 2006; Schenkel et al., 2006). So, the SNP detected in this study could benefit the further association studies to dissect the dynamics of genetic variant for the tenderness trait in cattle. Another example is the SNPs residing at a gene that is potentially implicated with marbling. Marbling refers to the amount of intramuscular fat present within the meat sample, and it is the most highly appreciated factor to evaluate the beef quality particularly in Japan and Korea. Currently, several candidate genes were found to have association with the intramuscular fat content in beef meat. Among the candidates, thyroglobulin (TG) encodes a dimeric glycoprotein hormone which is one of the largest proteins, and it was reported to be significantly associated with marbling score in Wagyubased cattle breeds (Casas et al., 2007). Among the nsSNP in this study, a total of five nsSNPs were found to reside at the TG (9,296,282; 9,377,082; 9,443,813; 9,469,822; and 9,483,749bp at chromosome 14). In this study, we identified numerous numbers of SNPs that could have functional implication with traits of interest in cattle, and these valuable resources are expected to be used in further studies for dissecting how each SNP affects meat quality traits in cattle.

 

Conclusions

In this article, a Korean native cattle breed, Heugu has been the first whole-genome resequenced using massively parallel sequencing technology. Throughout the genome, approximately 6 million SNPs were identified, and we bring a significant number of novel SNP markers (~1.8 million SNPs) that can be used for improving genetic potential of the meat quality in cattle. Further annotation was applied to all the detected SNPs and it shows numerous genes containing nsSNPs that could affect phenotypic characteristics of Heugu. The annotated SNPs identified in this genome sequencing study can serve as potential candidates in searches for phenotype-altering DNA difference implicated with meat quality traits in cattle. Notably, another main benefit of this study would be the raw sequence data of Heugu that can be used in further wholegenome comparisons between diverse native cattle breeds or different organisms.