Haplotype-Based Genome-Wide Association Study and Identification of Candidate Genes Associated with Carcass Traits in Hanwoo Cattle

Hanwoo, is the most popular native beef cattle in South Korea. Due to its extensive popularity, research is ongoing to enhance its carcass quality and marbling traits. In this study we conducted a haplotype-based genome-wide association study (GWAS) by constructing haplotype blocks by three methods: number of single nucleotide polymorphisms (SNPs) in a haplotype block (nsnp), length of genomic region in kb (Len) and linkage disequilibrium (LD). Significant haplotype blocks and genes associated with them were identified for carcass traits such as BFT (back fat thickness), EMA (eye Muscle area), CWT (carcass weight) and MS (marbling score). Gene-set enrichment analysis and functional annotation of genes in the significantly-associated loci revealed candidate genes, including PLCB1 and PLCB4 present on BTA13, coding for phospholipases, which might be important candidates for increasing fat deposition due to their role in lipid metabolism and adipogenesis. CEL (carboxyl ester lipase), a bile-salt activated lipase, responsible for lipid catabolic process was also identified within the significantly-associated haplotype block on BTA11. The results were validated in a different Hanwoo population. The genes and pathways identified in this study may serve as good candidates for improving carcass traits in Hanwoo cattle.


Introduction
Hanwoo, the premium beef cattle of South Korea, is well known for its juiciness, texture, and flavor. In order to maintain consumer preferences and profitability carcass and meat quality must be continuously enhanced [1]. Consumers prefer highly-marbled meat, as it increases tenderness and palatability, therefore, identifying markers associated with muscular fat deposition is imperative [2]. Traits such as carcass traits, reproductive traits, milk production traits, and many more can be easily analyzed by genomic region mapping or by identification of quantitative trait locus (QTL). Genomic loci associated with complex traits can be explored through technologies such as SNP (single nucleotide polymorphism) genotyping and high-throughput sequencing, together with GWAS (genome-wide association study). Many studies have explored the effect of SNPs on carcass traits in Hanwoo [3][4][5][6][7][8][9][10][11] including one [3], which explored genetic contribution of exonic SNPs on carcass traits in Hanwoo using linear mixed model and Bayesian networks. They also suggested that a large proportion of genetic variance was contributed by synonymous variants. Unlike single SNP GWAS, the use of haplotype blocks where each block is tagged with an effective SNP (SNP tag) in GWAS, enables high-resolution association and progeny test score [39]. Animals were slaughtered at approximately 24 months of age. Four carcass traits-back fat thickness (BFT), carcass weight (CWT), eye muscle area (EMA), and marbling score (MS) were analyzed in this study. CWT was measured in kilograms ranging from 158 to 518 kg. BFT was measured in millimeters ranging from 2 to 25 mm, EMA was measured in centimeter square ranging from 22 to 130 cm 2 and marbling score was assessed on a 1-9-point scale according to the Korean beef marbling standard (KAPE) ranging from 1 to 9 in given data. Descriptive statistics of the four carcass traits for 1161 Hanwoo cattle are given in Table 1 with boxplot for phenotypic traits in Figure S1. A total of 611,799 SNPs was genotyped with the Illumina Bovine SNP700K BeadChip for 898 animals. Quality control was performed through PLINK ver. 1.9 [40]. For further analyses, SNPs having low minor allele frequency (<0.01), low genotyping rate (<0.90), significant deviation from Hardy-Weinberg equilibrium (<0.001) were discarded; only biallelic format was consider for further study. Individuals with low genotyping call rate (<0.95) and missing slaughter age were excluded from the study. After filtering for quality, 887 animals with 547,836 SNPs remained for further analyses.

Construction of Haplotype Blocks
Haplotype blocks were created for 887 animals with 547,836 SNPs. For construction of haplotype blocks, two steps were followed: (1) phasing of chromosomes through SHAPEIT.v2 [41] and (2) construction of haplotype blocks through GHap, package in R.

Phasing of Chromosomes
The genotype data was segregated based on chromosomes with creation of individual chromosome ped and map files as explained below (Scheme 1). of slaughtered animals, so ethical committee approval was not required. Available information included sire and dam information along with the individual animals' date of birth, age at slaughter and progeny test score [39]. Animals were slaughtered at approximately 24 months of age. Four carcass traits-back fat thickness (BFT), carcass weight (CWT), eye muscle area (EMA), and marbling score (MS) were analyzed in this study. CWT was measured in kilograms ranging from 158 to 518 kg. BFT was measured in millimeters ranging from 2 to 25 mm, EMA was measured in centimeter square ranging from 22 to 130 cm 2 and marbling score was assessed on a 1-9-point scale according to the Korean beef marbling standard (KAPE) ranging from 1 to 9 in given data. Descriptive statistics of the four carcass traits for 1161 Hanwoo cattle are given in Table 1 with boxplot for phenotypic traits in Figure S1. A total of 611,799 SNPs was genotyped with the Illumina Bovine SNP700K BeadChip for 898 animals. Quality control was performed through PLINK ver. 1.9 [40]. For further analyses, SNPs having low minor allele frequency (<0.01), low genotyping rate (<0.90), significant deviation from Hardy-Weinberg equilibrium (<0.001) were discarded; only biallelic format was consider for further study. Individuals with low genotyping call rate (<0.95) and missing slaughter age were excluded from the study. After filtering for quality, 887 animals with 547,836 SNPs remained for further analyses.

Construction of Haplotype blocks
Haplotype blocks were created for 887 animals with 547,836 SNPs. For construction of haplotype blocks, two steps were followed: (1) phasing of chromosomes through SHAPEIT.v2 [41] and (2) construction of haplotype blocks through GHap, package in R.

Phasing of Chromosomes
The genotype data was segregated based on chromosomes with creation of individual chromosome ped and map files as explained below (Scheme 1). Scheme 1. Splitting of Genotypic data (SNP) based on chromosome.
The genotypes for each chromosome were phased using SHAPEIT.v2 with 200 states and window size of 0.5 Mb [41]. SHAPEIT was used for estimating haplotypes from either genotypes or sequencing data. This created the haps and sample file which was used to generate haplotype blocks through GHap. The genotypes for each chromosome were phased using SHAPEIT.v2 with 200 states and window size of 0.5 Mb [41]. SHAPEIT was used for estimating haplotypes from either genotypes or sequencing data. This created the haps and sample file which was used to generate haplotype blocks through GHap.

Haplotype Blocks (Haploblocks) Construction
Haplotype blocks were constructed based on LD, Len, and nsnp methods and corresponding phase, sample, and marker files for each chromosome were created. The genomic positions in haplotype blocks were based on the distance from the first SNP. After phasing, haplotype blocks were constructed. The blocks were based on: (1) the number of SNP haplotype blocks containing 2, 5, 10, 20, 30, and 50 SNPs per block, respectively, (2) genomic window sizes of 5 kb, 10 kb, 20 kb, 50 kb, 100 kb, and 200 kb, respectively, and (3) LD with r 2 value of 0.2, 0.3, 0.4, 0.5, 0.6, and 0.7 respectively. A total of 979,942 haplotype blocks were created using the above-mentioned three methods: LD, LEN, and nsnp. We then used GHap (genome-wide haplotyping) package in R to identify haplotype alleles from these haplotypes' blocks. After construction of haplotypes for each chromosome based on the above categories, all chromosomes were merged as shown below (Scheme 2): haplotype blocks were based on the distance from the first SNP. After phasing, haplotype blocks were constructed. The blocks were based on: (1) the number of SNP haplotype blocks containing 2, 5, 10, 20, 30, and 50 SNPs per block, respectively, (2) genomic window sizes of 5 kb, 10 kb, 20 kb, 50 kb, 100 kb, and 200 kb, respectively, and (3) LD with r 2 value of 0.2, 0.3, 0.4, 0.5, 0.6, and 0.7 respectively. A total of 979,942 haplotype blocks were created using the above-mentioned three methods: LD, LEN, and nsnp. We then used GHap (genome-wide haplotyping) package in R to identify haplotype alleles from these haplotypes' blocks. After construction of haplotypes for each chromosome based on the above categories, all chromosomes were merged as shown below (Scheme 2): Scheme 2. Haplotype blocks constructed based on nsnp, len and LD for each chromosome are merged to form one data set of that particular category. These data set are further used for association study.

Genome-Wide Association Study (GWAS)
Haplotype alleles determined by GHap were used for association studies using linear mixed model implemented in GEMMA (genome-wide efficient mixed model association) software [38]. As it is computationally efficient in performing large scale GWAS studies. Phenotype and genotype data for 887 animals were used to perform the studies for all four traits: BFT, CWT, EMA, and MS. The model, was used as follows where y is an n-vector of quantitative traits for n individuals and W = (w1; ….. ; wc) is an n x c matrix of covariates (fixed effects). Slaughter age was used as fixed effect in this study, as all animals were steers and reared in same slaughterhouse. Hence, sex and slaughter place are not considered as fixed effects. α is a c vector of the corresponding coefficient, x is a n-vector of haplotype, β is the effect size of the marker, and µ is an n-vector of random effects. GEMMA obtains maximum likelihood estimate (MLE) of and outputs the corresponding p-value. These p-values were used to filter out the significant haplotype blocks for each trait. Significant haplotype blocks were identified based on a Bonferronicorrected p-value (0.05/979942 = 5.10 × 10 −8 ).

Validation
To validate our result, the genotype and phenotype data of a further sample of 468 Hanwoo cattle, reported in a previous study, was used [11]. Similar steps were followed to construct haplotype blocks and to study association analyses as used for the original data set.

Gene Identification in Significant Haplotype Blocks
To identify candidate genes within the associated haplotype blocks, it was necessary to investigate the genomic region. Significant haplotype blocks were scanned from start to end position to identify the genes present. Ensembl database [42] (Ensembl 97: July 2019) with ARS-UCD1.2 was used as the reference genome and detailed gene information was obtained from the NCBI (National Center for Biotechnology Information) [43]. Information such as gene symbol, location, assembly method, gene type, and gene function were retrieved for each gene present in significant haplotype blocks for all four carcass traits (BFT, CWT, EMA, and MS).

Gene Functional Enrichment
Scheme 2. Haplotype blocks constructed based on nsnp, len and LD for each chromosome are merged to form one data set of that particular category. These data set are further used for association study.

Genome-Wide Association Study (GWAS)
Haplotype alleles determined by GHap were used for association studies using linear mixed model implemented in GEMMA (genome-wide efficient mixed model association) software [38]. As it is computationally efficient in performing large scale GWAS studies. Phenotype and genotype data for 887 animals were used to perform the studies for all four traits: BFT, CWT, EMA, and MS. The model, was used as follows where y is an n-vector of quantitative traits for n individuals and W = (w 1 ; . . . ..; w c ) is an n x c matrix of covariates (fixed effects). Slaughter age was used as fixed effect in this study, as all animals were steers and reared in same slaughterhouse. Hence, sex and slaughter place are not considered as fixed effects. α is a c vector of the corresponding coefficient, x is a n-vector of haplotype, β is the effect size of the marker, and µ is an n-vector of random effects. GEMMA obtains maximum likelihood estimate (MLE) of and outputs the corresponding p-value. These p-values were used to filter out the significant haplotype blocks for each trait. Significant haplotype blocks were identified based on a Bonferroni-corrected p-value (0.05/979942 = 5.10 × 10 −8 ).

Validation
To validate our result, the genotype and phenotype data of a further sample of 468 Hanwoo cattle, reported in a previous study, was used [11]. Similar steps were followed to construct haplotype blocks and to study association analyses as used for the original data set.

Gene Identification in Significant Haplotype Blocks
To identify candidate genes within the associated haplotype blocks, it was necessary to investigate the genomic region. Significant haplotype blocks were scanned from start to end position to identify the genes present. Ensembl database [42] (Ensembl 97: July 2019) with ARS-UCD1.2 was used as the reference genome and detailed gene information was obtained from the NCBI (National Center for Biotechnology Information) [43]. Information such as gene symbol, location, assembly method, gene type, and gene function were retrieved for each gene present in significant haplotype blocks for all four carcass traits (BFT, CWT, EMA, and MS).

Gene Functional Enrichment
After identification of genes present in significant haplotype blocks, the PANTHER classification system [44] was used to classify and cluster genes based on biological process. To understand and identify the biological pathways and the genes involved, gene ontology (GO) enrichment was also performed using DAVID (Database for Annotation, Visualization, and Integrated Discovery) v6.8 [45,46]. Significantly-enriched pathways and GO terms were identified based on enrichment scores.

Significant Haplotype Blocks
A total of 1434 haplotype blocks (BFT: 150, CWT: 1062, EMA: 161, and MS: 61) were identified as significantly-associated haplotype blocks based on this threshold (Figure 1). Detail tabulated description of each haplotype block for four carcass traits (BFT, CWT, EMA and MS) is provided in Table S1. Most significant regions were detected on BTA14 for CWT; BTA11 and BTA13 for BFT; BTA1, BTA6, BTA9 and BTA12 for MS; BTA1, and BTA10 for EMA. Manhattan plot for four carcass traits by all three method are provided in Table S2.
Genes 2020, 11, x FOR PEER REVIEW 5 of 18 BTA1, BTA6, BTA9 and BTA12 for MS; BTA1, and BTA10 for EMA. Manhattan plot for four carcass traits by all three method are provided in Table S2.

Candidate Gene Identification
Genes within the significant haplotype blocks were identified using Ensembl database for all three methods (LD, Len, and nsnp) for each trait. Venny 2.1.0 [47] revealed the number of genes which are unique or common among three methods for all carcass traits ( Figure 2, Table S3).

Candidate Gene Identification
Genes within the significant haplotype blocks were identified using Ensembl database for all three methods (LD, Len, and nsnp) for each trait. Venny 2.1.0 [47] revealed the number of genes which are unique or common among three methods for all carcass traits ( Figure 2, Table S3).
weight, EMA: eye muscle area, and MS: marbling score.

Candidate Gene Identification
Genes within the significant haplotype blocks were identified using Ensembl database for all three methods (LD, Len, and nsnp) for each trait. Venny 2.1.0 [47] revealed the number of genes which are unique or common among three methods for all carcass traits ( Figure 2, Table S3). Unique and common trait-specific genes were analyzed for each trait and are listed along with function in Table 2. It can be observed that for BFT only five genes were identified by all three methods whereas a greater number of unique genes was identified by the Len followed by the LD methods. PLCB1 and TMX4 were some of the genes identified by all three methods whereas PLCB4, is identified only by Len and LD methods and the CEL gene i.e., carboxyl ester lipase, by the LD method. For CWT, the number of common genes detected by all three methods was 45% of total genes and these included LYPLA1, TMEM68, XKR9, PREK2, KLHL38, TOX and CYP7A1, and SOX17 and PRKDC, all present on BTA14. MYC and FAM110B genes responsible for the cell cycle and cell progression were detected only by the nsnp method.   Unique and common trait-specific genes were analyzed for each trait and are listed along with function in Table 2. It can be observed that for BFT only five genes were identified by all three methods whereas a greater number of unique genes was identified by the Len followed by the LD methods. PLCB1 and TMX4 were some of the genes identified by all three methods whereas PLCB4, is identified only by Len and LD methods and the CEL gene i.e., carboxyl ester lipase, by the LD method. For CWT, the number of common genes detected by all three methods was 45% of total genes and these included LYPLA1, TMEM68, XKR9, PREK2, KLHL38, TOX and CYP7A1, and SOX17 and PRKDC, all present on BTA14. MYC and FAM110B genes responsible for the cell cycle and cell progression were detected only by the nsnp method. For EMA, SLC8A3 was detected by all three methods while COL18A1 and RORA were detected by only LD and nsnp methods, respectively. MS showed only FYN gene as common, whereas four genes were common between LD and nsnp method including PGM2. RB1 and PAPR14 emerged as unique gene detected only by Len and LD methods, respectively.

Gene Function Annotation
PANTHER was used to functionally annotate genes based on biological, cellular, and metabolic processes as in Figure 3. Further, all the genes identified by the three methods were functionally annotated via DAVID. DAVID clusters the gene dataset based on the enrichment score (Table S4,  Table S5(a-d)). The higher the enrichment score, the more significant the cluster. The result in Table 3 shows that genes such as CEL, PLCB1, PLCB4, present in significant haplotype blocks are mainly associated with the lipid catabolic process (GO: 0016042), insulin secretion (bta04911), the glucagon-signaling pathway (bta04922), insulin secretion by pancreas (bta04911), glutamatergic synapses (bta04724), and the sphingolipid-signaling pathway (bta04071) etc., which are involved in lipid metabolism and therefore can help in increasing the richness of meat. RORA, helps in regulation of fat cell differentiation (GO: 004559). Meat quality is determined by fat, protein, and cell mass. Several pathways that might play a role in increasing cell mass and size were identified by DAVID and included cell cycle (bta04110), negative regulation of G1/S transition of the mitotic cell (GO: 2000134), potassium ion trans-membrane transportation (GO: 0071805), and regulation of microtubule cytoskeleton organization (GO: 0070507). Several other related GO terms were also identified ( Table 3).

Discussion
The association study by GEMMA software revealed that BTA 13 and BTA11 are significant haplotype blocks for BFT as they contain markers which may affect candidate genes such as PLCB1, PLCB4, and CEL. Figure 2 highlights that there were numerous genes present in significant haplotype blocks, some identified commonly by all three methods, some by two methods, and others were uniquely detected by single method. BTA14 was found to have most of the significantly-associated haplotype block for CWT ( Figure 1) which agrees with previous reports [3,5]. It contains a greater number of common genes extracted by all three methods-approximately 45% of total genes, as seen in Figure 2b. Table 4 highlights the candidate genes found in significantly-associated haplotype blocks. BTA 12 was found to be more significant for EMA. BTA 12, by the LD method, BTA 9 and BTA 12, by the Len method, and BTA 3 and BTA 12, by the nsnp method, were found to be more significant for MS, as can be seen through the Manhattan plot (Table S2) and p-value generated by GEMMA software. The significant haplotype blocks identified in this study were validated in a different population of 468 Hanwoo cattle ( Table 4). The significance threshold for the validation population was 3.45 × 10 −8 . SNPs in BTA13 and BTA11 that harbored the majority of the haplotype blocks and candidate genes significant to this study were found to be associated with carcass traits in the test population (887 original data).
Functional annotation analyses with PANTHER and DAVID, showed that most of the identified genes were found to function in cellular and metabolic processes. Genes which are involved in increasing cell mass, cell number, protein formation, and stabilization and hence are mass determinant are listed in Table 3 [2]. Cellular processes such as hyperplasia and hypertrophy are responsible for increasing muscle mass [2]. Hypertrophy can be increased by increased protein synthesis. LYPLA1, TMEM68, and XKR9 which were found within the significantly-associated haplotype blocks were identified as candidate genes for improving carcass traits in a previous study [52]. These genes are responsible for lipase activity and weight gain in cattle. LYPLA1 gene has hydrolase and lipase activity. TMEM68 is expressed in bovine rumen, abomasum, intestine and adipose tissue in cattle, and likely affects lipid biosynthetic processes. XKR9 promotes phosphatidyl serine exposure on apoptotic cell surface. KLHL38 gene plays an important role in maintaining the balance between intermediate filament stability, a process in maintaining skin integrity. It is an important component of ubiquitin-protein ligase complex also formed by KLHl2 [52] responsible for proteasomal degradation. Previous studies [3,5] had reported that TOX and CYP7A1 were functional genes for CWT in Hanwoo cattle. Pausch et al. [55] identified SOX17 as positional and functional gene involved in transcriptional regulation via inhibiting Wnt signaling. RB1 [58,59] plays an important role in posttranslational modifications and hence regulates cell cycle. These genes are responsible for cell development and proliferation. Another set of gene detected were involved in muscular growth such as KLHL38 (filament stability) and KCNQ3, KCNB2, and KCNIP4 [30] which play an important role in calcium and potassium transport and in meat tenderization through their involvement in the proteolytic system responsible for postmortem tenderization and muscle contraction and these were among the candidate genes. A study by Zhao et al. [60], suggested that modulation of annexin expression in animals with high BFT may play an important role in fat formation in beef cattle. These annexins are calcium-dependent phospholipid-binding proteins involved in inhibition of phospholipase activity, exocytosis, endocytosis, signal transduction, and many other biological processes. As can be seen through the PANTHER result there is a set of lipase genes that participates in metabolic processes and plays a significant role in formation of kinases and hence metabolism of the lipid molecule. Through our study PLCB1, phospholipase C Beta 1, and PLCB4 phospholipase C Beta 4 [48,61] were identified as important phospholipases by 12 haplotype blocks and three haplotype blocks, respectively, as seen in Table 4. Haploview analysis of top haplotype blocks (based on p-value) for each candidate gene is shown in Figure 4. There were five SNPs contributing to top haplotype block, nsnp5 with block id 68104 for PLCB1 gene. There were 58 SNPs present in haplotype block Len200 with block id: 7021 affecting gene PLCB4. These phospholipases play an essential role in lipid metabolism as fat deposition plays an important role in increasing richness in meat and it has been reported that very low levels of intramuscular fat can lead to dry and less-tasty meat [62]. CEL gene i.e., carboxyl ester lipase is a bile-salt-activated lipase, detected by only one haplotype block with block id: 18600 constructed by LD0.2 method with haplotype block start position 103051196 to end 103319019 with 36 SNPs (Figure 4). These three genes are responsible for the lipid catabolic-process glucagon signaling pathway, metabolic pathways, and glutamatergic pathways, etc., as can be seen in Table 3. The CYP7B1 gene present on BTA 14 also has a minor role in bile acid synthesis [54] which plays an essential role in lipid metabolism in intestines of animal.
Phospholipase Important phospholipase C β (PLCB) isoforms were identified in this study i.e., PLCB1 and PLCB4. These genes code for 1-phosphatidylinositol 4, 5-bisphosphate phosphodiesterase proteins which are responsible for the hydrolysis of phospholipids [48]. As discussed above fat deposition is governed by the cell mass and cell size and these phospholipids play an important role in forming major components of the cell membrane and determining the physical and chemical properties of the cell [63]. When these phospholipids are acted upon by phospholipase C they form diglyceride (DG) and inositol phosphates for Ca 2+ , which is an activator of protein kinase C [64]. They are further responsible for catabolism of lipid and metabolism of membrane structural phospholipids [65] (GO: 0016042) as seen with the DAVID with p < 0.02. Their catabolic activity can reduce the fat deposition in muscle, making meat dry and lean. The genes coding for phospholipase are also involved in the sphingolipid signaling pathway (bta04071). Sphingolipids protect the cell surface from the harmful outside environment and help in making the cell mechanically stable. On the positive side, phospholipase governs the interaction of glucose, acetylcholine, and free fatty acids on insulin secretion (bta04911, bta04972: pancreatic secretion) as they increase insulin secretion leading to increased storage of glucose in body cells [66] and suppressing its counter-regulatory pathway i.e., the glucagon signaling pathway (bta04922) and other metabolic pathways (bta01100). Studies also showed that in ruminants, the adipose tissue is the principal site of de novo fatty acid (FA) synthesis [67,68] where glucose is the main substrate for lipogenic adipocytes in muscle [69].

Phospholipase
Important phospholipase C β (PLCB) isoforms were identified in this study i.e., PLCB1 and PLCB4. These genes code for 1-phosphatidylinositol 4, 5-bisphosphate phosphodiesterase proteins which are responsible for the hydrolysis of phospholipids [48]. As discussed above fat deposition is governed by the cell mass and cell size and these phospholipids play an important role in forming major components of the cell membrane and determining the physical and chemical properties of the cell [63]. When these phospholipids are acted upon by phospholipase C they form diglyceride (DG) and inositol phosphates for Ca 2+ , which is an activator of protein kinase C [64]. They are further responsible for catabolism of lipid and metabolism of membrane structural phospholipids [65] (GO: 0016042) as seen with the DAVID with p < 0.02. Their catabolic activity can reduce the fat deposition in muscle, making meat dry and lean. The genes coding for phospholipase are also involved in the sphingolipid signaling pathway (bta04071). Sphingolipids protect the cell surface from the harmful outside environment and help in making the cell mechanically stable. On the positive side, phospholipase governs the interaction of glucose, acetylcholine, and free fatty acids on insulin secretion (bta04911, bta04972: pancreatic secretion) as they increase insulin secretion leading to increased storage of glucose in body cells [66] and suppressing its counter-regulatory pathway i.e., the glucagon signaling pathway (bta04922) and other metabolic pathways (bta01100). Studies also showed that in ruminants, the adipose tissue is the principal site of de novo fatty acid (FA) synthesis [67,68] where glucose is the main substrate for lipogenic adipocytes in muscle [69].

Conclusions
PLCB1, PLCB4, and CEL genes were identified as putative candidate genes, as they regulate lipid metabolism which is responsible for increasing fat thickness and the weight of an organism, and hence more rich and marbled meat. To our knowledge, until now, no work has been reported which

Conclusions
PLCB1, PLCB4, and CEL genes were identified as putative candidate genes, as they regulate lipid metabolism which is responsible for increasing fat thickness and the weight of an organism, and hence more rich and marbled meat. To our knowledge, until now, no work has been reported which showed involvement of PLCB1, PLCB4, and CEL in increasing the carcass traits of an animal. Also, we can conclude that the LD and nsnp methods are more useful in the construction of haplotype blocks as a bigger number of unique genes are obtained by these two methods. For example, CEL was only detected by the LD method and PLCB4 only detected by the nsnp and Len methods. Alteration in expression of these lipase-forming genes can help in improving Hanwoo beef quality by increasing muscular fat thickness leading to enhancement of economic growth.

Supplementary Materials:
The following are available online at http://www.mdpi.com/2073-4425/11/5/551/s1: Table S1: Significant haploblocks identified through GEMMA; Table S2: Manhattan plot for 887 for association study through GEMMA for all methods; Figure S1: Boxplot for phenotype used in this study; Table S3: Common genes identified by all three methods: Len, LD, nsnp, Table S4: DAVID functional annotation cluster for genes and their pathways; Table S5(a-d). DAVID functional annotation chart summary of significant genes for four carcass traits identified based on p-value.
Author Contributions: S.S. performed the analysis and drafted the manuscript; D.L. was involved in conceptualization, funding acquisition, project administration, resources, and supervision; J.-H.S. was responsible for imputation of 50K and 777K genotype data; S.W. was responsible for phenotypic data collection; K.S. contributed towards network analysis; and GWAS; J.-E.P., W.P., and H.-H.C. contributed to quality control of genotype data; All authors read and agreed to the published version of the manuscript.