Hanwoo is an indigenous and popular meat-type cattle in Korea, and is particularly renowned for its rapid growth rate and quality attributes such as juiciness, tenderness, characteristic flavor, and extensive marbling of its beef (Jo et al., 2012). In recent years, both carcass and meat quality traits in Hanwoo have been extensively studied because of their economic relevance for optimizing the profitability of the beef industry. The current selection index in Hanwoo focuses on the improvement of carcass traits, such as backfat thickness (BFT), carcass weight (CW), eye muscle area (EMA), and marbling score (MS) as major selection criteria for breeding programs (Kim et al., 2017). However, other important traits such as carcass primal cuts have received inadequate attention in the Hanwoo breeding program, which affects both the quantity and quality of meat, and consequently, command premium prices. To meet consumer demand, the importance of primal cut yields is growing in the beef industry of developed countries because of its market value. Hence, cattle breeders need to address these traits, which determine selection decisions to increase carcass cut-out value and consumer acceptance of meat. Meanwhile, the existence of genetic variation and moderate to high heritability in the yield of primal cuts has been reported (Choi et al., 2015). In this sense, the improvement of primal cut yields requires knowledge of the underlying genetic background influencing these invaluable traits.

Over the last decade, with the development of high-throughput single nucleotide polymorphism (SNP) genotyping technologies, genome-wide association studies (GWAS) have become an affordable and powerful tool for detecting and localizing candidate genes and causal mutations associated with quantitative traits in different species (Matukumalli et al., 2009). Several statistical methods to conduct GWAS have been developed and applied, among which a simple regression model has been widely used, where one marker is tested at a time for significance (Meyer and Tier, 2012). However, this method was challenged by false positives and overestimation of quantitative trait loci (QTL) effects. Therefore, the marker effect models in the Bayesian approaches have been proposed for GWAS analysis (Habier et al., 2011; Wang et al., 2012; Moser et al., 2015; Wang et al., 2016) as they offer methods to overcome these challenges (Strömberg, 2009; Peters et al., 2012). One such method could have a higher power to detect SNPs with moderate effects on a trait of interest. In addition, Bayesian methods are flexible in accounting for the uncertainties of variables and parameters and allow for inferences to be made by finding their marginal posterior distributions (Yi and Shriner, 2008; Blasco and Blasco, 2017). Recently, Fernando et al. (2014; 2016) developed a class of single-step Bayesian regression methods (ssBR), which not only combines all available information as single-step genomic best linear unbiased prediction [ssGBLUP; (Misztal et al., 2009)] does, but also accommodates Bayesian models. This method can also be extended to GWAS and controls the proportion of false positives by computing the posterior probability of association of a trait with each SNP or each window of consecutive SNPs. Numerous association studies have been carried out on growth and carcass traits in beef cattle using different GWAS approaches (Peters et al., 2012; Lee et al., 2013; Magalhaes et al., 2016; Weng et al., 2016; Roberts, 2018; Bedhane et al., 2019; Naserkheil et al., 2020). However, GWAS have not yet been conducted to identify significant genomic regions for carcass primal cut traits, which are highly relevant to Hanwoo cattle breeding. Hence, the objective of this study was to perform GWAS to detect genomic regions and candidate genes associated with primal cut yields in Hanwoo cattle using the ssBR approach.

The number of records, means, minimum, maximum, standard deviations, phenotypic coefficient of variation, and heritability estimates for the 10 primal cut traits are provided in Table 1. The mean values of these traits ranged from 6.04 to 57.55 with standard deviation between 0.76 and 7.53. The coefficients of variation ranged from 11.89 to 25.72%, indicating considerable phenotypic variation of the investigated traits in the Hanwoo cattle population. Moderate to high heritability estimates were obtained for primal cut yields, which ranged from 0.32 ± 0.03 to 0.61 ± 0.04. Estimated variance components are also summarized in Supplementary Table S1.

The identification of genomic regions related to the traits of interest was performed using a ssBR approach. Only windows with WPPA ≥0.8, were considered significant. Manhattan plots displayed the proportion of genomic variance (Figure 1), additive genetic variance (Supplementary Figure S2) explained and the posterior probability of association by each 1-Mb window for the studied traits (Figure 2). Positional candidate genes for primal cut traits were also detected within the significant windows. A summary of significant windows associated with the traits under study, such as the number of SNPs in each window, explained GV% or AGV%, and WPPA, as well as candidate genes, is shown in Table 2. A total of 16 relevant genomic regions (1-Mb window) were found to be associated with the 10 traits recorded in this study. These regions were distributed over nine different chromosomes: 3, 4, 6, 10, 11, 14, 16, 19, and 25. Of these significant windows, five genomic windows were pleiotropic QTL, meaning that the QTL had an effect on multiple traits, which were located on chromosomes 6 at 38–39, 11 at 21–22, 14 at 6–7 Mb and 26–27, and 19 at 26–27 Mb. The genomic window located on chromosome 14 at 26–27 Mb explained a large proportion of the genomic variance across the nine analyzed traits, including tenderloin, sirloin, striploin, chuck, brisket, top round, bottom round, shank, and flank. The largest QTL window was observed for shank, which explained 20.51% of the genomic variance, was located in the region of 26–27 Mb on chromosome 14. The QTL window with the smallest proportion of genomic variance (0.87%) was identified for bottom round, located on chromosome 14 at 6–7 Mb. Candidate genes responsible for the genomic variance explained by the 1-Mb window were identified using the Bos taurus genome map. A total of 154 genes were identified within the significant regions to be associated with the traits of interest. Of these, 92 genes were codified proteins, 3 were miRNA (microRNA), 6 were tRNA (RNA transporter), and 53 were pseudogenes (Table 2). Functional enrichment analysis revealed 18 biological processes, six molecular functions, seven cellular components, and seven KEGG pathways. The following biological process terms were highlighted: growth (GO:0040007), positive regulation of growth (GO:0045927), tissue development (GO:0009888), muscle structure development (GO:0061061), skeletal system development (GO:0001501), animal organ development (GO:0048513), cellular response to chemical stimulus (GO:0070887), lipid metabolic process (GO:0006629), response to lipid (GO:0033993), glycolipid biosynthetic process (GO:0009247), and cell cycle (GO:0007049) being the most significant for the traits under study. In the case of the molecular function, calcium channel activity (GO:0005262), protein binding (GO:0005515), ligase activity (GO:0016874), lipid binding (GO:0008289), catalytic activity (GO:0003824), and nucleotide binding (GO:0000166) were significant terms. The significant enrichments for cellular component, including endoplasmic reticulum membrane (GO:0005789), intracellular (GO:0005622), endoplasmic reticulum (GO:0005783), actin cytoskeleton (GO:0015629), extracellular region (GO:0005576), organelle (GO:0043226), and cytoplasmic region (GO:0099568) were obtained. Moreover, KEGG pathway analysis revealed that the identified candidate genes involved in primal cut yields were enriched in metabolic pathways (bta01100), focal adhesion (bta04510), ECM-receptor interaction (bta04512), glycolysis/gluconeogenesis (bta00010), fat digestion and absorption (bta04975), Rap1 signaling pathway (bta04015), and calcium signaling pathway (bta04020). Functional gene set annotation and enrichment pathways are presented in Table 3.

The aim of this study was to identify genomic regions associated with primal cut yields using a ssBR approach in Hanwoo cattle. The marker effect model in a Bayesian framework would seem to be useful for GWAS because they account for uncertainty in parameters required to construct posterior distributions for QTL inference, thereby improving accuracy of genomic predictions and the power of QTL detection (Fernando and Garrick, 2013; Mehrban et al., 2017). In recent years, interest in exploring genomic regions that control economically important traits in beef cattle has increased due to advances in high-throughput genotyping techniques and the constant availability of molecular data, statistical methods, and ease of application of GWAS. Primal cut traits have recently been proposed as potential indicator of carcass weight and overall carcass merit (Berry et al., 2019) given that the genetic correlations between these traits and carcass weight are generally moderate to strong (Choi et al., 2015; Judge et al., 2019). Nevertheless, selection for the weight of primal cuts requires knowledge of the genetic basis for these traits, which may be useful in future genomic evaluations targeting the improvement of weight in the more valuable primal cuts, and consequently increasing the profitability of the meat production system. Our results showed that primal cut yields were moderate to highly heritable, being in accordance with those reported in Hanwoo cattle (Choi et al., 2015), Chianina cattle (Sarti et al., 2013), Simmental cattle (Zhu et al., 2019), and Irish cattle (Berry et al., 2019; Judge et al., 2019). Among all identified window regions, the QTL on chromosome 14 at position 26–27 Mb had a larger impact than any of the other QTLs and was associated with a greater number of traits. This region, which is related to tenderloin, sirloin, striploin, chuck, brisket, top round, bottom round, shank, and flank, explained between 4.26 and 20.51% of the genomic variance across all these traits. A total of 11 genes were detected on chromosome 14 at 26–27 Mb regions. Among these, FAM110B, UBXN2B, NSMAF, TOX, SDCBP, and CYP7A1 were notable. This region also seems to be most significantly associated with carcass and growth traits in beef cattle (Magalhaes et al., 2016; Roberts, 2018; Zhang et al., 2019; Grigoletto et al., 2020). Positional candidate genes of FAM110B, UBXN2B, NSMAF, CYP7A1, SDCBP, and TOX have been previously reported to be associated with carcass weight in Hanwoo cattle (Lee et al., 2013; Bhuiyan et al., 2018; Naserkheil et al., 2020). For instance, Lee et al. (2013) identified the six most significant SNPs associated with carcass weight in Hanwoo that were located in or nearby TOX, FAM110B, and SDCBP. Similarly, Srikanth et al. (2020) reported that the most significant SNPs on chromosome 14 were located in UBXN2B, CYP7A1, SDCBP, and TOX, which have been regarded as positional candidate genes for carcass weight in Hanwoo cattle. It was also reported that CYP7A1 and SDCBP are positional candidate genes for carcass weight and eye muscle area in Hanwoo cattle (Bhuiyan et al., 2018; Srivastava et al., 2020), weaning weight in Brangus cattle (Weng et al., 2016), and feed efficiency traits in Nellore cattle (Brunes et al., 2021). The TOX gene located within this region is associated with reproductive traits in Nellore (de Camargo et al., 2015), residual feed intake and mid-test metabolic weight in SimAngus (Seabury et al., 2017), and development of puberty in Brahman cattle (Fortes et al., 2012). In addition, the UBXN2B gene was found to be associated with mid-test metabolic weight in SimAngus (Seabury et al., 2017) and carcass weight, carcass fat, and carcass conformation in Simmental cattle (Purfield et al., 2019), which is known as a protein-coding gene involved in endoplasmic reticulum biogenesis. FAM110B has been reported to be associated with fat thickness in composite beef cattle (Roberts, 2018). This gene functions in the cell cycle and cell growth and might play an important role in increasing carcass weight in beef cattle by increasing cell number and size. A pleiotropic QTL on chromosome 6 at 38–39 Mb was associated with tenderloin, striploin, brisket, top round, bottom round, and shank, which explained the largest (6.46%) and smallest (2.7%) proportion of genomic variance for bottom round and tenderloin, respectively. This region harbors relevant candidate genes, including PKD2, SPP1, MEPE, IBSP, LAP3, NCAPG, and LCORL. Most of the positional genes detected on chromosome 6 have previously been associated with many economically important traits in beef and dairy cattle. In a study on Brangus beef cattle, Weng et al. (2016) reported that most positional genes associated with direct birth weight, weaning weight, and yearling weight are located on chromosome 6. Similarly, Saatchi et al. (2014) identified a large-effect pleiotropic QTL located on chromosome 6 at 37–42 Mb was associated with direct birth weight, calving ease, carcass weight, rib eye muscle area, and weaning weight across several cattle breeds. In addition, it has previously been reported that a large number of significant SNPs associated with skeleton trait in Simmental cattle that were harbored by LAP3, FAM184B, LCORL, and NCAPG genes (Xia et al., 2017). A major QTL on chromosome 6, extending from 36 to 39 Mb related to carcass weight, was also identified in Japanese Black cattle (Nishimura et al., 2012). Interestingly, the NCAPG and LCORL genes, while being associated with the skeletal type traits (Hoshiba et al., 2013; Doyle et al., 2020) have also been regarded as positional candidate genes for direct calving ease, feed intake, gain, meat, and carcass traits (Lindholm-Perry et al., 2011; Bongiorni et al., 2012; Purfield et al., 2019), as well as growth and lipid deposition (Snelling et al., 2010; Weikard et al., 2010) across several breeds. Moreover, Liu et al. (2015) identified LCORL as a positional gene related to weight and carcass composition traits in chickens based on GWAS and differentially expressed gene studies. Other notable positional candidate genes in this region including PKD2, LAP3, SPP1, MEPE, IBSP and MED28, have been associated with carcass weight and growth traits (Lindholm-Perry et al., 2011; Weng et al., 2016; Roberts, 2018; Naserkheil et al., 2020), milk production (Olsen et al., 2005), and reproductive traits and puberty (Daetwyler et al., 2008; Cánovas et al., 2014) in beef and dairy cattle.

The other large-effect pleiotropic QTL identified in multiple traits (tenderloin, sirloin, brisket, top round, bottom round, shank, and flank) was located on chromosome 19 at position 26–27 Mb. This region explained between 1 and 2.07% of the genomic variance across the traits of interest, and harbors candidate genes, including PITPNM3, WSCD1, MIS12, and RABEP1. A recent study on linear type traits conducted by Doyle et al. (2020b) identified PITPNM3 as a potential gene for chest depth in Angus cattle. The WSCD1 gene encodes a protein with sulfotransferase activity involved in glucose metabolism, and is associated with udder depth in dairy cattle (Li et al., 2021), feed efficiency and feeding behaviors in pigs (Guo et al., 2015), and body size in beef cattle (An et al., 2020). The MIS12 and RABEP1 genes have also been recently linked to milk production traits in dairy cattle (Cai et al., 2020; Li et al., 2021). Two other pleiotropic genomic windows were located on chromosome 11 at 21–22 Mb, and on chromosome 14 at position 6–7 Mb, which the proportion of genomic variance explained by these windows ranged from 1.13 to 1.69 and from 0.87 to 1.9, respectively. The region on chromosome 14 that was found to be associated with sirloin, brisket, bottom round, and shank had no positional candidate genes in cattle, whereas its orthologous region on human chromosome 8 containing KHDRBS3 gene that may play a role as a negative regulator of cell growth and inhibition of cell proliferation (Ma et al., 2017). For QTL region on chromosome 11, a total of 15 genes annotated were found to be related to brisket, top round, bottom round, and shank. Overall, 11 chromosomal regions identified in this study were trait-specific QTLs for six traits. Three trait-specific QTLs were identified for tenderloin, which are distributed on chromosomes 4 at 77–78 Mb, 6 at 37–38 Mb, and 10 at 32–33 Mb. The trait-specific QTL on chromosome 6 was responsible for 2.74% of the genomic variance in tenderloin and harbors the positional candidate genes PPM1K, ABCG2, and PIGY. The PPM1K gene is involved in cellular survival, phosphorus metabolic process, amino acid dephosphorylation, and development by regulating mitochondrial permeability transition pore function, which have been shown to be associated with increased carcass weight, mid-point metabolic weight, reduction of residual feed intake, feed efficiency conversion ratio, and marbling score in crossbred beef cattle (McClure et al., 2010). The ABCG2 gene is related to body weight (Weng et al., 2016) and milk yield and composition traits (Cohen-Zinder et al., 2005), which is involved in iron transport and metabolism. PIGY is a member of the PIG gene family, which encodes the glycosylphosphatidylinositol-N-acetylglucosaminyltransferase (GPI-GnT) complex and plays an important role in cell–cell interactions. A previous study reported that there were significant effects of copy number variation of the PIGY gene on growth traits in three Chinese sheep breeds (Feng et al., 2020). Brisket also had three trait-specific QTLs located on chromosomes 3 at 98–99, 10 at 49–50, and 25 at 32–33 Mb and explained 1.27, 0.92, and 0.92% of the genomic variance for this trait, respectively. The top round had two trait-specific QTLs. One was located on chromosome 4 at 6–7 Mb (GV% = 1.46), and the other on chromosome 25 at 4–5 Mb (GV% = 0.93). A QTL for the bottom round was identified on chromosome 16 at 25–26 Mb, which accounted for 1.13% of the genomic variance and harbors the positional candidate gene DUSP10. Previous studies have shown that the DUSP10 gene is associated with growth traits (Ribeiro et al., 2021) and carcass weight (Chang et al., 2018) in beef cattle. The QTL located on chromosome 6 at position 39–40 Mb was associated with shank and had a greater proportion of genomic variance (4.89%) than other trait-specific QTLs. Only a QTL on chromosome 4 at 8–9 Mb, which was responsible for 3.23% of the genomic variance for the rib, was identified. The CDK14 gene is located in this region and is associated with fatty acid profile (C18:1 trans-9) in the intramuscular fat of Longissimus thoracis muscle of Nellore cattle (Lemos et al., 2016). According to these results, most primal cut traits are probably controlled by several QTLs with large effects. Among these, genomic regions located on chromosomes 6 and 14 were considered as hot spots for several causal variants related to many economically important traits in beef cattle. A similar conclusion was given by Mehrban et al. (2021), who identified major QTLs on chromosomes 6 and 14 for EMA, yearling weight, and particularly for CW using the weighted single-step GWAS in Hanwoo cattle.

Gene ontology and pathway enrichment analyses were carried out to gain insight into the genes identified within QTL windows using g:Profiler and DAVID functional classification clustering tools (Table 3). Our analyses revealed the significant GO terms classified into the biological processes, cellular components, molecular functions, and seven KEGG pathways were enriched for the studied traits. It is interesting to note that the majority of the common genes identified for the primal cut traits are involved in growth-related processes: growth, positive regulation of growth, muscle structure development, regulation of cell development, muscle organ development, tissue development, skeletal system development, biomineral tissue development, and animal organ development. Among them, six genes, namely HLX, SOS1, SDCBP, ANXA2, FZD1, and MEPE, were highlighted as the main candidates for the traits under study in at least three biological pathways. The HLX gene, located on chromosome 16 at approximately 25–26 Mb, is a protein-coding gene that is involved in embryogenesis and hematopoiesis. It has also been shown to be associated with intramuscular fat in composite beef cattle (Roberts, 2018). The SDCBP gene is located on a conserved region on chromosome 14. It encodes a protein that binds to a variety of transmembrane proteins and plays a crucial role in carcass and meat quality traits. ANXA2 gene is known to encode a member of a widely distributed, phospholipid-binding, calcium-regulated, peripheral membrane protein family known as annexins. This gene is involved in molecular functions related to calcium channel activity and lipid binding. It is not surprising then, that the role of calcium in meat tenderness and muscle contraction, and is a key regulator of muscle growth in beef cattle (Sadkowski et al., 2009). ANXA2 knockout affects white adipose tissue hypotrophy due to reduced fatty acid uptake in mice (Salameh et al., 2016). Furthermore, this gene is associated with feed conversion efficiency in beef cattle (Al-Husseini et al., 2014). The GO terms related to the lipid metabolic process, cellular responses to the chemical stimulus, cell cycle, localization, and regulation of fat cell differentiation were also significantly represented in the pleiotropic QTL windows (Table 3). Pathway enrichment revealed that eight genes from six window regions (located on chromosomes 4, 6, 11, 14, and 25) were significantly associated with the metabolic pathway (bta01100). Among the genes harbored in this pathway, CYP7A1 is involved in the transport, synthesis, and secretion of cholesterol, steroids, and other lipids (Zhao et al., 2013), which play a crucial role in digestion and absorption of dietary fat and contribute in maintaining the balance of cholesterol and lipid metabolism within the body (Monte et al., 2009). Interestingly, the extracellular matrix (ECM)-receptor interaction (bta04512) and focal adhesion (bta04510) pathways were enriched in the genes SPP1 and IBSP from a pleiotropic QTL located on chromosome 6 at 38–39 Mb, which functions as a positive regulator in skeletal muscle cells and is involved in bone mineralization processes. The ECM-receptor interaction pathway plays a key role in tissue, organ morphogenesis, cell maintenance, and tissue structure and function (Mariman and Wang, 2010). It was previously reported that the ECM-receptor is upregulated in subcutaneous fat and intramuscular fat and appears to be involved in adipogenesis and meat tenderness (Taye et al., 2018). Focal adhesion also participates in important biological processes and serves as a mechanical link to ECM receptors and other molecules. These results will improve our understanding of enriched molecular processes, pathways, and genes associated with the primal cut traits and shed some light on how different pathways control these traits.

In the current study, 16 genomic regions (SNP windows) were found to be associated with 10 primal cut traits in Hanwoo cattle using a single-step Bayesian regression GWAS. Within these windows, five QTLs had pleiotropic effects, with the most significant region located on chromosome 14 at position 26–27 Mb. Several candidate genes with potential functions in tissue development, regulation of growth, and lipid metabolism for the related traits were highlighted, among which SPP1, IBSP, PKD2, SDCBP, PIGY, CYP7A1, and MEPE were well-known. Moreover, our findings can contribute to a better understanding of the genetic basis and biological processes underlying the traits of interest; consequently, information on QTL regions can be used to search for causal mutations and marker-assisted or genomic selection in Hanwoo breeding schemes.