Implications
In this study, we identified genetic associations of marbling score (MS) with two intragenic nucleotide variants at selection signals of Korean cattle. This implied that selection for MS has resulted in selection signals on the bovine genome. Favorable alleles of the variants would be incorporated to marker-assisted selection for MS. Further efforts are needed to confirm the associations of the markers with MS using other breeds and to make fine maps in and around the putative candidate genes responsible for MS.
Introduction
Meat quality has been a major concern in the genetic improvement of Korean cattle (Lee and Pollak, Reference Lee and Pollak2002). In particular, marbling is the primary factor for determining meat quality and price. A selection index has been applied since the 1980s to improve Korean cattle reflecting pedigree-based genetic merits of MS (Kim and Lee, Reference Kim and Lee2000). Recent efforts have been made to identify MS genetic markers for marker-assisted selection. Genetic markers were discovered within candidate genes encoding mitogen-activated protein kinase 6, peroxisomal proliferator-activated receptor –γ coactivator-1α and sterol regulatory element binding factor 1 (P<0.05; Ryu et al., Reference Ryu, Kim, Kim, Kim and Lee2012; Lee et al., Reference Lee, Oh, Lee, La and Yeo2013). Genetic association was examined more intensively through a genomewide association study (GWAS), but rs109593638 was the only nucleotide variant associated with MS (P<1.27×10−6; Kim et al., Reference Kim, Ryu, Woo, Kim, Kim and Lee2011).
Recently, a study employing composite log likelihood with a sliding window of five consecutive nucleotide variants revealed 147 selection signals of Korean cattle as regions of causal variants and variants in their vicinity with allele frequencies augmented by selection (Ryu and Lee, Reference Ryu and Lee2014). The selection signals may be attributed largely to artificial selection undergone to improve meat quality in recent decades. Unexpectedly, few genetic associations at the most probable selection signal of each chromosome were found (17%, 5 out of 30 selection signals; Ryu and Lee, Reference Ryu and Lee2014).
Nowadays, cutting-edge technologies reduce sequencing costs, leading to remarkable progress in the discovery of genes and nucleotide variants in the bovine genome. Genetic variants of newly annotated genes in the bovine genome would be helpful to overcome limited findings of intragenic variants associated with meat quality at selection signals in the previous study. In this study, we examined genetic associations of MS with updated genes and variants at the most probable selection signal of each chromosome in Korean cattle.
Material and methods
Phenotypes
A total of 411 Korean cattle were used to examine genetic associations of MS with intragenic nucleotide variants in the most significant selection signal on each chromosome. All animals used as subjects were steers raised under a controlled feeding program in the farm of the Korean Cattle Improvement Center and slaughtered at the same slaughter house from 2008 to 2010. They were slaughtered at 871 days, and MS was scored from 1 to 9 (mean: 6.7) using the left carcass cut across the vertebra between the last thoracic vertebra and the first lumbar vertebra. A larger MS score signifies more abundant intramuscular fat. We used genotypes of 309 Korean cattle from a previous GWAS (Kim et al., Reference Kim, Ryu, Woo, Kim, Kim and Lee2011). Tissue samples of longissimus muscle were consecutively collected from the other 102 Korean cattle and genotyped in the current study.
Genotypes
Genomewide nucleotide variants were genotyped with Illumina BovineSNP50 BeadChip using Beadstudio software (Illumina, San Diego, CA, USA). Animals were filtered with a call rate >98%. Variants were filtered with a call rate >98%, minor allele frequency >0.05 and Hardy–Weinberg disequilibrium (P<0.01). All the variants filtered by the criteria were utilized for background polygenic effects. Intragenic variants selected at selection signals were individually tested for genetic association with MS. The candidate genes and the corresponding intragenic variants were determined based on a recently updated database, the bovine genome Btau_4.6.1 build. The current study utilized a new set of variants that had not been dealt with in the study of Ryu and Lee (Reference Ryu and Lee2014).
Genetic association analysis
Genetic associations were analyzed employing the mixed model with random polygenic effect based on identical-by-state genetic relationships among animals in order to control for spurious associations produced by population structure. The analytical model incorporated the leaving-one-chromosome-out approach as follows:
where y is the vector of observed phenotypes; µ the overall mean; 1 the vector of 1’s; β the fixed effect for the minor allele of the candidate nucleotide variant; x the vector with elements of 0, 1 and 2 for the homozygote of the minor allele, heterozygote and homozygote of the major allele; and g¯ the vector of random polygenic effects using the genome without the chromosome in which the candidate variant is located. That is,
$${\bf{g}} ^{{\bf{{\minus}}} } {\bf{\,\sim\,}} N(0,{\bf{G}} ^{{\bf{{\minus}}} } \sigma _{{{\rm g}^{{\rm {\minus}}} }}^{2} )$$
where G¯ is the n×n genomic relationship matrix with elements consisting of pairwise genetic relationship coefficients estimated excluding one chromosome and
$\sigma _{{{\rm g}^{{\minus}} }}^{2} $
the polygenic variance component that should be re-estimated whenever a specific chromosome is excluded from the calculation of the genetic relationship matrix. The pairwise relationship coefficient was calculated as:
where N
c
is the number of chromosomes (=30), N
m
the number of variants in the mth chromosome, p
ij
the frequency of the minor allele at the jth variant in the ith chromosome and n
ijk
(n
ijl
) represents the number (0, 1 or 2) of minor alleles for the kth (lth) individuals. ε the vector of random errors
$({\bf {\epsilon}}\,\sim\,N(0,{\mib{I}} \sigma _{{\epsilon}}^{2} ))$
, where
$\sigma _{{\rm {\epsilon}}}^{2} $
the error variance component and
I
the identity matrix. The polygenic and error variance components were estimated using the restricted maximum likelihood (REML). We first estimated polygenic and error variance components by EM-REML, then the estimates were used as initial values in order to obtain their AI-REML estimates. The fixed candidate nucleotide variant effect was then obtained with the estimated variance components under the mixed model equations. Multiple testing by Bonferroni correction was applied to the association analysis. The statistical analysis was conducted using Genome-wide Complex Trait Analysis freeware (version 2.0; Yang et al., Reference Yang, Lee, Goddard and Visscher2011).
Replication analysis
Genetic associations were validated with an independent data set of Korean cattle. We used 327 Korean cattle slaughtered between 2007 and 2009. The MS was scored according to the same method as described above. Their tissue samples were obtained from the Hongcheon Federation of Livestock Cooperation in Korea, and genomic DNA was extracted from the tissue samples using a QIAamp Tissue Kit (Qiagen, Hilden, Germany). Nucleotide variants were genotyped using the TaqMan polymerase chain reaction assay (Applied Biosystems, Foster City, CA, USA). Reactions were carried out following the manufactures protocol, and the products resulted from the reactions were analyzed using ABI PRISM 7900HT (Applied Biosystems). Genetic associations with MS were analyzed using linear regression under the assumption of an additive contribution of minor allele. Multiple testing by Bonferroni correction was also applied to the association analysis.
Results
A total of 39 092 nucleotide variants of 407 Korean cattle with an average call rate of 99.98% were obtained after quality control of the data. Genetic associations of MS were examined with a total of 129 intragenic markers: 21 updated markers and 108 markers in 88 newly updated genes (Table 1). Six variants located within cache domain containing 1 (CACHD1), like-glycosyltransferase (LARGE), cell division cycle 42 binding protein kinase alpha (CDC42BPA), snurportin 1 (SNUPN), protein tyrosine phosphatase, non-receptor type 9 (PTPN9) and chondroitin sulfate proteoglycan 4 (CSPG4) genes were identified at selection signals of Korean cattle (P<3.88×10−4, Table 2). In particular, MS was independently associated with each of the genes encoding SNUPN, PTPN9 and CSPG4 within the most probable selection signal of chromosome 21. All variants were located in introns except the CSPG4 variant, which was located in an exon. Further genetic association analysis using independent Korean cattle confirmed the association of MS with the intragenic variants of CDC42BPA and LARGE genes (P<8.33×10−3, Table 3).
Genes located at the most probable selection signal of each chromosomeFootnote a
BTA=Bos taurus autosome.
a Genes are listed based on Btau 4.6.1. Bold type indicates newly updated genes on Btau 4.6.1.
b No genes are located at selection signals of chromosomes 1 and 15. No selection signals were identified on chromosome 12 by employing a chromosome-wide significance threshold.
Genetic associations of marbling score with nucleotide markers at selection signalsFootnote a

BTA=Bos taurus autosome; SNP=single nucleotide polymorphism; MAF=minor allele frequency; P HWE=P-value estimated for Hardy–Weinberg equilibrium test; CDC42BPA=cell division cycle 42 binding protein kinase alpha; LARGE=like-glycosyltransferase; SNUPN=snurportin 1; PTPN9=protein tyrosine phosphatase, non-receptor type 9; CSPG4=chondroitin sulfate proteoglycan 4; CACHD1=cache domain containing 1.
a Only significant SNPs by Bonferroni correction are presented (P<3.88×10−4).
b Major/minor allele.
c Regression coefficient estimated for minor allele using an additive model.
Validation of nucleotide markers associated with marbling score using an independent data setFootnote a

BTA=Bos taurus autosome; SNP=single nucleotide polymorphism; MAF=minor allele frequency; P HWE=P-value estimated for Hardy–Weinberg equilibrium test; CDC42BPA=cell division cycle 42 binding protein kinase alpha; LARGE=like-glycosyltransferase; SNUPN=snurportin 1; PTPN9=protein tyrosine phosphatase, non-receptor type 9; CSPG4=chondroitin sulfate proteoglycan 4; CACHD1=cache domain containing 1.
a Significant SNPs by Bonferroni correction (P<8.33×10−3) are in bold.
b Major/minor allele.
c Regression coefficient estimated for minor allele using an additive model.
Discussion
The current study revealed two intragenic markers associated with MS at selection signals of Korean cattle. They were identified with genes defined by rapidly updated bovine genome annotations after a discovery of five markers from a previous study (Ryu and Lee, Reference Ryu and Lee2014). As a result, a total of seven intragenic markers have been discovered to be associated with MS at selection signals of Korean cattle. The post-selection signal association analyses were more efficient than the genome-wide association analysis (Kim et al., Reference Kim, Ryu, Woo, Kim, Kim and Lee2011), which revealed only one nucleotide marker (rs109593638, human immunodeficiency virus type I enhance binding protein 3 (HIVEP3) intron 2, P<1.01×10−6) associated with MS. This is practical evidence that intensive directional selection for MS has resulted in selection signals on the bovine genome by augmenting allele frequencies of functional variants and their proximity variants. Nevertheless, no genetic associations were observed at selection signals of the other 21 chromosomes. We suspect that some signals might be associated with physiological and morphometric traits of Korean cattle since Korean cattle moved from China and have been native to the Korean peninsula without introducing additional germplasm for over 2000 years (Lee and Pollak, Reference Lee and Pollak2002).
The genes in which nucleotide variants were associated with MS in this study may be considered putative candidates responsible for the meat quality of Korean cattle. In particular, CDC42BPA has been reported as having a genetic association with porcine meat quality. Expression of the gene was correlated with meat quality of crossbred produced by Duroc and Pietrain, known to differ in meat quality (Ponsuksili et al., Reference Ponsuksili, Murani, Phatsara, Schwerin, Schellander and Wimmers2009). The gene was identified as a selection signal when comparing domestic pigs to wild boars after candidate genes had been screened from a gene coexpression network analysis using lean and obese pig breeds (Zhao et al., Reference Zhao, Liu and Liu2015). Although any other associations between LARGE gene and meat quality of cattle had not been found in previous studies, the gene also might be also a candidate for the meat quality considering its role of glycosylation of α-dystroglycan in skeletal muscle (Kanagawa et al., Reference Kanagawa, Saito, Kunz, Yoshida-Moriguchi, Barresi, Kobayashi, Muschler, Dumanski, Michele, Oldstone and Campbell2004). The glycosylation by LARGE is a key process to form a link between subsarcolemmal cytoskeleton and laminin of extracellular matrix in maintaining skeletal muscle. Lack of LARGE induced loss of laminin-binding activity, resulting in muscular dystrophy of mice (Barresi et al., Reference Barresi, Michele, Kanagawa, Harper, Dovico, Satz, Moore, Zhang, Schachter, Dumanski, Cohn, Nishino and Campbell2004).
The present study employed a mixed model with random polygenic effects based on identical-by-state genetic relationships. The analytical method explains population stratification better than the fixed model which many GWAS have used with genomic inflation factors (Kang et al., Reference Kang, Sul, Service, Zaitlen, Kong, Freimer, Sabatti and Eskin2010; Ryoo and Lee, Reference Ryoo and Lee2014; Zhang et al., Reference Zhang, Ersoz, Lai, Todhunter, Tiwari, Gore, Bradbury, Yu, Arnett, Ordovas and Buckler2010). Most of all, mixed model methodology can reduce false positive and negative associations. Simulation studies showed that the mixed model performed better in terms of false discovery rate and statistical power than the fixed model (Shin and Lee, Reference Shin and Lee2015a and Reference Shin and Lee2015b).
Acknowledgments
The valuable comments from Dr Aduli Malau-Aduli and two anonymous reviewers on an earlier version of this article are greatly appreciated. This study was supported by the Basic Science Research Program of the National Research Foundation of Korea (NRF) funded by the Ministry of Education, Science, and Technology (Grant No. 2012002096). The authors thank our laboratory personnel who participated in collecting tissues of Korean cattle and extracting genomic DNA from the tissues.

