In recent years, various molecular genetics studies have revealed that livestock breeding is quite often dependent on individual phenotypical values or pedigree values (Dang et al., 2011). Developments in molecular genetics in the 2000s have led to the discovery of useful genes and improved research on domestic animal breeding while confirming the relationship between polymorphisms of these major genes and meat traits (Seong et al., 2014). Meat traits and genetic polymorphisms in livestock are mainly analyzed using single-nucleotide polymorphisms (SNPs) (Gill et al., 2009). Genetic polymorphisms in the genome include duplication, deletion, inversion, translocation, and insertion. Copy number variation (CNV) is a genomic structural variation, the size of which is not clearly defined; it appears because of duplications or deletions ranging from a size of 50 bp to several Mbps, affecting the phenotype of the individual.

CNV detection is now possible using next-generation sequencing (NGS) (Xie and Tammi, 2009) and microarray (Sebat et al., 2004). Some domestic animals, such as cattle (Liu and Bickhart, 2012), sheep (Yang et al., 2018), and pigs (Ramayo-Caldas et al., 2010) have been studied for CNV research. In particular, the relationship between CNV and expressive traits has been analyzed in cattle. Upadhyay et al. (2017) identified the CNV pattern across European bovine populations and confirmed the presence of CNV polymorphisms among various European bovine populations. da Silva et al. (2016) analyzed the quantitative trait loci (QTL) present in Nellore cattle CNVs and 34 CNVs associated with production in Holstein milk. Hanwoo cattle is an inherent species raised in the Republic of Korea, and in the past, this species was often used for transport and cultivational activities. Since the 1960s, agricultural machinery has taken over most of the work, and Hanwoo has been devoted to meat cattle. In particular, while producing Korean beef, meat has important traits, such as marbling (Lee et al., 2010), longissimus muscle area (Lee and Kim, 2009), backfat thickness (Seong et al., 2014), carcass weight (Lee et al., 2013), and the corresponding genes. Studies on SNPs in Hanwoo have been conducted and have undergone many improvements using molecular genetics (Lee et al., 2012). In addition, there have been studies on CNV in Hanwoo using NGS (Choi et al., 2013) or SNP chips (Bae et al., 2010), but an analysis of the relationship between CNV in Hanwoo and the meat traits is lacking.

The main purpose of this study was to identify CNV and evaluated the relationship between CNVR in Hanwoo steers using Hanwoo SNP 50K chip. In conclusion, our results certified that the CNVs in Hanwoo are essential to their meat traits. These results contributed to the greater understanding of CNV in Hanwoo and its role in genetic variation among cattle livestock.

Hanwoo cattle information

We collected carcass information (phenotypic value) from the Korea Institute for Animal Products Quality Evaluation website (https://www.ekape.or.kr/) for 473 Hanwoo steers born between 2015 to 2018.

DNA extraction

DNA was extracted from the tissue of the tail hair root. gDNA was extracted using the MagMAX^TM DNA Multi-Sample Ultra 2.0 Kit (ThermoFisher, US). The quality of the extracted gDNA was confirmed using Epoch (Biotek, US).

Genome-wide SNP genotyping using Hanwoo SNP 50k chip

We used the Illumina HanwooSNP50 ver1 BeadChip (Illumina, US) to identify 53,866 SNP. HanwooSNP50 ver1 is based on bovineSNP50 ver3. All signal intensity log R ratio (LRR) and allelic intensity B allele frequency (BAF) ratios of the samples were reported using Illumina GenomeStudio software (Illumina, US).

CNV and CNVR calling

We used PennCNV (v1.0.5) (Wang et al., 2007), based on a hidden Markov model, to determine the number of copies and genotypes of each CNV. First, the analysis required a PFB file which was created via “compile_pfb.pl”. Further, the LRR and BAF values for each sample were calculated. Raw CNV data were used to calculate these values using the ‘Detect_cnv.pl’ option ‘-lastchr 29’. After CNV calling, samples that were filtered out had a high-intensity noise (LRR SD > 0.3), BAF drift < 0.01, and waviness factor > 0.05. We identified CNVRs using CNVRuler (v1.3.3.2) (Kim et al., 2012). We chose the CNVR method and the threshold was set at 0.5, and the CNVRs were sorted into gain, loss, and mixed. In addition, drawing CNVRs in chromosomes used a Rldeogram of the R package.

Gene annotation

Genes overlapping with CNVRs were retrieved from the Ensembl Gene 94 database (http://asia.ensembl.org/index.html), using the BioMart tool. HanwooSNP50K ver. 1, with UMD3.1 (genome sequence assembly), was used as the reference genome and the genes in CNVR with different frequencies were verified. Additionally, DAVID Bioinformatic Database (https://david.ncifcrf.gov/) was used to annotate the significantly associated regions. Gene ontology (GO) (http://geneontology.org/) that used three GO categories (biological processes, molecular function, and cellular component) and KEGG pathway (https://www.genome.jp/kegg/) analyses were performed.

QTL analysis

To investigate the potential associations between CNVRs and economically important traits in Hanwoo, the animal QTL database (https://www.animalgenome.org/cgi-bin/QTLdb/BT/index) was used on UMD 3.1. We downloaded the bed file of the cattle QTL database and mapped the CNVRs of this research with at least one bp overlapping.

Association with meat traits

We analyzed the p-value of CNVs using CNVRuler and graphically represented with the statistical software R-4.1.2 (https://cran.r-project.org/). We analyzed the association between the minor allele threshold 0.05 CNVR and four meat traits: carcass weight, longissimus muscle area, backfat thickness, and marbling score, using the CNVRuler program linear regression method. The following model was used.

Y_ij is each meat trait (carcass weight, backfat thickness, longissimus muscle area, and marbling score), μ is the mean, Ag_ei is the carcass month, CNV_j is the effect CNVs of each sample, and e_ij is the random error. The CNVR was determined by the significant effect of each of the four meat characteristics if the p-value < 0.05.

Identification of CNVs

We identified 2,575 autosomal CNVs in the 473 Hanwoo, out of which 933 and 1,642 CNVs were found to be losses and gains, respectively. The size of the identified CNVs ranged from 2,183 bp to 983,333 bp, with an average length of 132,848 bp. The distribution of CNV sizes with 22.0% being ~50K, 30.9% being 50K-100K, 30.1% being 100K-200K, 12.8% being 200K-300K, 2.5% being 300K-400K, and 1.7% 400K~ were relatively rare (Table 1).

Table 1. Size distributions of copy number variations (CNVs).

Size (bp)	Count	Percentage (%)
~50K	567	22.0
50K-100K	795	30.9
100K-200K	775	30.1
200K-300K	330	12.8
300K-400K	64	2.5
400K~	44	1.7
Total	2,575	100

Identification of CNVRs

A total of 416 CNVRs were obtained by merging the overlapping CNVs using the CNVRuler software which contained 289 gains (29,629,042 bp), 111 losses (13,370,341 bp), and 16 mixed events (1,937,662 bp). The chromosome with the highest number of CNVRs is chr1 and chr4, both having 27 CNVRs. On the other hand, the smallest number of CNVRs is chr25, which has 2 gains and 1 loss. The chromosomes with the highest number of gain CNVRs are chr1 and chr5, each having 17 gain CNVRs. Furthermore, these Hanwoo cattle CNVRs ranged in size from 10,004 bp to 381,836 bp meaning 108,022 bp (Table 2 and 3). Fig. 1 summarizes the map of CNVRs illustrating the distribution of autosomal chromosomes constructed using the Hanwoo genome and the Hanwoo SNP 50K bead chip.

Table 2. Distribution of copy number variable regions (CNVRs) in the Hanwoo genome.

Chr	Gain	Loss	Mixed	Total	Total length (bp)
1	17	10	0	27	3,397,875
2	15	4	0	19	1,736,068
3	15	3	0	18	1,805,976
4	16	8	3	27	3,003,020
5	17	4	0	21	1,949,611
6	14	9	1	24	2,568,299
7	14	3	2	19	2,225,824
8	7	2	0	9	948,864
9	5	4	0	9	1,129,687
10	13	1	0	14	1,181,146
11	12	5	2	19	1,925,950
12	5	9	2	16	2,135,892
13	7	3	0	10	952,130
14	12	4	0	16	1,573,673
15	11	6	0	17	1,336,528
16	11	2	0	13	1,406,768
17	12	5	1	18	2,241,899
18	12	4	0	16	2,193,728
19	14	3	0	17	1,660,771
20	5	6	0	11	1,360,317
21	11	4	0	15	1,453,179
22	7	2	1	10	807,742
23	7	2	1	10	1,187,854
24	7	1	1	9	783,710
25	2	1	0	3	419,047
26	6	1	0	7	1,002,764
27	2	2	1	5	801,027
28	7	1	0	8	832,252
29	6	2	1	9	915,444
Total	289	111	16	416	44,937,045

Table 3. Size distribution of CNVRs.

Size (bp)	Count	Percentage (%)
~50K	71	17.1
50K-100K	164	39.4
100K-200K	138	33.2
200K-300K	35	8.4
300K-400K	8	1.9
Total	416	100

Figure 1.CNVR in Hanwoo chromosome. Green, purple, and orange represent loss, gain, and mixed events, respectively.

Association of CNVRs with meat traits

The proportion of polymorphisms was analyzed for the association between meat traits and 5% or more CNVR area. CNVR_138, CNVR_185, CNVR_322, CNVR_388, CNVR_406, and CNVR_88 were associated with carcass weight, and CNVR_138 and CNVR_322 were associated with longissimus muscle area. CNVR_138, CNVR_185, and CNVR_406 were associated with backfat thickness, and CNVR_357 and CNVR_88 were associated with the marbling score (Table 4).

Table 4. Association of CNVRs with traits.

Trait	CNVR	chr	p-value
Carcass weight	CNVR_138	7	2.58E-06
	CNVR_322	19	0.00186
	CNVR_185	10	0.005047
	CNVR_406	28	0.010267
	CNVR_388	26	0.011912
	CNVR_88	4	0.026506
Longissimus muscle area	CNVR_138	7	0.008473
	CNVR_322	19	0.043607
Backfat thickness	CNVR_138	7	0.004984
	CNVR_185	19	0.01287
	CNVR_406	28	0.025993
Marbling score	CNVR_357	22	0.011328
	CNVR_88	4	0.039865

Gene annotation

BioMart tools were employed to map genes to the 416 identified CNVRs, and 347 genes were found to overlap with 347 CNVRs. These included 309 protein-coding genes, 10 snRNA genes, 9 miRNA genes, 9 pseudogenes, 5 snoRNA genes, 4 rRNA genes, and 1 processed_pseudogene. Furthermore, we analyzed genes with overlapping CNVRs. Findings overlapping were found between Obscurin (OBSCN) gene and CNVR_138, solute carrier family 8 member A3 (SLC8A3) gene and CNVR_185, mannose receptor C type 2 (MRC2) gene and CNVR_322. Lastly, proteasome assembly chaperone 3 (PSMG3), transmembrane protein 184A (TMEM184A), integrator complex subunit 1 (INTS1), and MICAL like 2 (MICALL2) genes overlapped with CNVR_388 as well (Table 5). In addition, to search for CNVs including gene functions, we analyzed enrichment protein-coding genes. The GO analysis showed 20 molecular functions, 11 cellular components, and 15 biological processes (Table 6). As well as, pathway analysis revealed that the number of genes associated with hepatitis C, metabolic pathways, ether lipid metabolism, and bile secretion was 7, 26, 4, and 4, respectively (Table 7).

Table 5. CNVRs are associated with meat traits, overlapping genes, and QTLs.

CNVR	Overlapped QTL ID	Overlapped gene	Associated traits
CNVR_88	125340		Carcass weight, marbling score
CNVR_138	106634	OBSCN	Carcass weight, back fat thickness, longissimus muscle area
CNVR_185	35235, 44964, 44965, 44966	SLC8A3	Carcass weight, back fat thickness
CNVR_322	56361-56368, 56636-56640, 20292, 49981-49991	MRC2	Carcass weight, longissimus muscle area
CNVR_357	56580, 56581, 127091, 122114		Marbling score
CNVR_388	24678, 32337	PSMG3, TMEM184A, INTS1, MICALL2	Carcass weight
CNVR_406	20829		Carcass weight, back fat thickness

Table 6. Gene ontology terms of genes associated with the identified CNVRs.

Term	Count	Percentage (%) of genes	p-value
Molecular function
GO:0016787~hydrolase activity	50	18.45	0.002
GO:0004722~protein serine/threonine phosphatase activity	5	1.85	0.005
GO:0016788~hydrolase activity, acting on ester bonds	19	7.01	0.009
GO:0005216~ion channel activity	13	4.80	0.013
GO:0022838~substrate-specific channel activity	13	4.80	0.017
GO:0022839~ion gated channel activity	4	1.48	0.024
GO:0005261~cation channel activity	10	3.69	0.025
GO:0015267~channel activity	13	4.80	0.025
GO:0022803~passive transmembrane transporter activity	13	4.80	0.025
GO:0016646~oxidoreductase activity, acting on the CH-NH group of donors, NAD or NADP as acceptor	3	1.11	0.028
GO:0015269~calcium-activated potassium channel activity	3	1.11	0.031
GO:0022836~gated channel activity	10	3.69	0.032
GO:0098811~transcriptional repressor activity, RNA polymerase II activating transcription factor binding	4	1.48	0.032
GO:0001190~transcriptional activator activity, RNA polymerase II transcription factor binding	4	1.48	0.032
GO:0036094~small molecule binding	45	16.61	0.032
GO:0052689~carboxylic ester hydrolase activity	6	2.21	0.033
GO:0000166~nucleotide binding	42	15.50	0.040
GO:1901265~nucleoside phosphate binding	42	15.50	0.040
GO:0046965~retinoid X receptor binding	2	0.74	0.046
GO:0042578~phosphoric ester hydrolase activity	10	3.69	0.050
Cellular component
GO:0005737~cytoplasm	141	52.03	0.005
GO:0005773~vacuole	22	8.12	0.008
GO:0005769~early endosome	8	2.95	0.018
GO:0005774~vacuolar membrane	12	4.43	0.022
GO:0005622~intracellular	183	67.53	0.023
GO:0005768~endosome	15	5.54	0.024
GO:0044424~intracellular part	176	64.94	0.026
GO:0044437~vacuolar part	12	4.43	0.029
GO:0098588~bounding membrane of organelle	23	8.49	0.043
GO:0015630~microtubule cytoskeleton	21	7.75	0.048
GO:0005815~microtubule organizing center	14	5.17	0.049
Biological process
GO:0007033~vacuole organization	8	2.95	0.010
GO:1901700~response to oxygen-containing compound	21	7.75	0.014
GO:0010508~positive regulation of autophagy	5	1.85	0.017
GO:0010738~regulation of protein kinase A signaling	3	1.11	0.023
GO:0050728~negative regulation of inflammatory response	5	1.85	0.030
GO:0055086~nucleobase-containing small molecule metabolic process	14	5.17	0.035
GO:0044281~small molecule metabolic process	29	10.70	0.037
GO:0010506~regulation of autophagy	7	2.58	0.042
GO:0006182~cGMP biosynthetic process	3	1.11	0.042
GO:0009152~purine ribonucleotide biosynthetic process	7	2.58	0.044
GO:0043149~stress fiber assembly	4	1.48	0.046
GO:0030038~contractile actin filament bundle assembly	4	1.48	0.046
GO:0006164~purine nucleotide biosynthetic process	7	2.58	0.048
GO:0009165~nucleotide biosynthetic process	8	2.95	0.049
GO:0009409~response to cold	3	1.11	0.049

Table 7. KEGG pathways of genes associated with identified CNVRs.

Term	Count	Percentage (%) of genes	p-value
bta05160: Hepatitis C	7	2.58	0.006
bta01100: Metabolic pathways	26	9.59	0.006
bta00565: Ether lipid metabolism	4	1.48	0.019
bta04976: Bile secretion	4	1.48	0.050

QTL analysis

In addition, we identified 2,550 QTLs that overlap with CNVRs, including 1,026 milk traits, 461 production traits, 447 reproduction traits, 240 exterior traits, 221 meat and carcass traits, and 155 health traits (Table 8). The main factors were QTL, shear force, intramuscular fat, subcutaneous fat, and tridecylic acid content in the meat and carcass traits (Table 9). The QTL in CNVR revealed an association with meat traits (Table 5): body weight in CNVR_138, the efficiency of gain, milk fat percentage, calving ease and stillbirth in CNVR_185, fat thickness at the 12th rib in CNVR_322, marbling score and milk conjugated linoleic acid percentage in CNVR_388, and shear force in CNVR_406 overlapped. In particular, CNVR_138 has a body weight QTL that is considered to be associated with carcass weight. CNVR_322 was associated with carcass weight and longissimus muscle area but overlapped thickness at the 12th rib QTL. CNVR_388 had a carcass weight trait association and marbling score QTL. We assumed that CNVs are directly or indirectly associated with quantitative traits.

Table 8. QTLs found to overlap with CNVRs identified in this study.

QTL	Count
Reproduction trait	447
Production trait	461
Meat and carcass traits	221
Health traits	155
Exterior trait	240
Milk traits	1,026
Total	2,550

Table 9. Meat and carcass traits overlapped with the CNVRs identified in this study.

Meat and carcass traits	Count of CNVR
Meat color L	3
Shear force	49
Lean meat yield	5
Myristic acid content	6
Subcutaneous rump fat thickness	6
Pelvic area	8
Myristoleic acid content	5
Tenderness score	1
Intramuscular fat	13
Longissimus muscle area	9
Fat thickness at the 12th rib	9
Marbling score	5
Docosahexaenoic acid content	5
Lignoceric acid content	8
Omega-3 unsaturated fatty acid content	5
Carcass weight	4
Iron content	3
Stearic acid content	1
Oleic acid content	6
Linolenic acid content	1
Monounsaturated fatty acid content	3
Intermuscular fat percentage	1
Elaidic acid content	1
Yield grade	6
Subcutaneous fat	14
Tridecylic acid content	13
Dihomo-gamma-linolenic acid content	4
Eicosapentaenoic acid content	2
Trans-6/9-C18:1 fatty acid content	1
Margaric acid content	9
Palmitoleic acid content	5
Medium-chain fatty acid content	2
Palmitic acid content	1
Long-chain fatty acid content	1
Saturated fatty acid content	1
Atherogenic index	1
Myristic and palmitic acid ratio	1
Muscle protein percentage	2
cis-Vaccenic acid content	1

Two studies analyzed the CNVs in Hanwoo cattle using an SNP chip, of which Bae et al. (2010) used the bovine 50K SNP chip to analyze 265 cattle and confirmed 368 CNVRs. Among our results, we found that 89 CNVRs overlapped with Bae et al. (2010) results. However, there are differences in Bae et al. (2010) studies and methods of CNV calling and how to configure CNVR, so it is somewhat impossible to identify Hanwoo-specific CNVs. In addition, Bae et al. (2010) selected a Btau 4.0 reference for their analysis, which may introduce differences when compared to results obtained using LiftOver and UMD 3.1. On the other hand, Shin and Oh (2016) analyzed 571 Hanwoo cattle were analyzed using a BovineHD chip, resulting in the identification of 1,659 CNVRs. It was confirmed CNVR of up to 50 kb in the whole CNVR, is approximately 52%. The number of CNVRs is relatively different in this study compared to the research by Bae et al. (2010) and Shin and Oh (2016). There is enough possibility that different results will be obtained even if the same livestock group is analyzed because the limit of the range that determines the size of CNV is the algorithm for estimating CNV and the SNP density of the platform used (Henrichsen et al., 2009). CNV analysis using SNP chips has been studied in various standards, but the distance between the SNP markers and the basic limits, such as density, is a problem yet to be overcome. Generally, the marker with the Hardy-Weinberg equilibrium can be excluded from the SNP chip, which distinguishes the potential CNV using the SNP chip that is otherwise difficult to identify. CNV identification and analysis require more studies to reduce these deviations by sample and reduce these deviations.

In addition, we determined genes with overlapping CNVRs. CNVR_138 overlapped in the OBSCN gene, which is known to play a role in adipogenesis and lipid metabolism (Jager et al., 2013) and previous reports have highlighted the occurrence of imprinted fat accumulation resulting from the early weaning of calves in the Hanwoo breed (Reddy et al., 2017). Also, the OBSCN gene is expected to affect fat synthesis, which is considered to have an indirect effect on carcass results. Furthermore, CNVR_185 is a candidate in the SLC8A3 gene; it is a candidate gene that is estimated to affect phenotypes, such as meat quality in Hanwoo and Qinchuan (Lim et al., 2013). This suggests that SLC8A3 has been identified in the path of the endoplasmic reticulum membrane and has a cooking loss, softness, juiciness, and coupled tissue relevance (Leal-Gutiérrez et al., 2019). Also, CNVR_322 is located in the MRC2 gene. The MRC2 gene affects tail syndrome, but in Belgian Blue cattle, this gene has a major effect on the muscle meat phenotype (Druet et al., 2014). Of the CNVR_388 genes, the relationship between CNV in MICALL2 and growth traits was analyzed (Xu et al., 2013).

As well as, GO analysis revealed primarily associated with the cytoplasm and intracellular vacuoles, of which the cytoplasm is reported to be useful in assessing hygiene and freshness in food production (Kalyuzhnaya et al., 2020). The cytoplasm appears to affect meat quality. In addition, substances, such as intracellular calcium signaling, apoptosis, and cathepsin, can be affected by the extracellular matrix, which has been reported as one of the potential mechanisms affected by fleshy development (Xing et al., 2019). It is also known that intracellular Ca²⁺ concentration in pigs and cattle affects meat quality (Küchenmeister and Kuhn, 2003). Therefore, the result suggested that cytoplasm and vacuoles are significantly associated with meat traits.

In conclusion, we identified 2,575 CNVs and 416 CNVRs using the Hanwoo SNP 50K chip in 473 Hanwoo steers. In addition, CNV is associated with four meat traits: carcass weight, longissimus muscle area, backfat thickness, and marbling score were analyzed and verified. It was confirmed that there are six CNVRs (CNVR_138, CNVR_322, CNVR_185, CNVR_406, CNVR_388, CNVR_88) in carcass weight, two CNVRs (CNVR_138, CNV_322) in longissimus muscle area, three CNVRs (CNVR_138, CNVR_185, CNVR_406) in backfat thickness, and two CNVRs (CNVR_357 and CNVR_88) which were significantly associated with marbling score. We expect that our results can serve as a valuable asset for future explorations into CNVs and contribute to a better understanding of the relationship between CNVs and meat’s important traits in cattle. Among our results, three CNVRs (CNVR_138, CNVR_185, and CNVR_322) were determined to have a gene that affects meat quality. These results suggest that the CNVs in Hanwoo are associated with their meat traits. Further research on the same would prove to be beneficial in the future.

None.

Conceptualization, H.S.K.; methodology, C.M.B., K.T., G.H.L., and G.J.S.; investigation, C.M.B., and G.J.S.; data curation, C.M.B., K.T., G.H.L., and G.J.S.; writing^—original draft preparation, C.M.B., K.T., G.H.L., and G.J.S.; writing^—review and editing, C.M.B., K.T., G.H.L., G.J.S., and H.S.K.; supervision, H.S.K.; project administration, H.S.K.; funding acquisition, H.S.K.

None.

The study was approved by the Hankyong National University Animal Ethics Committee (No. 2021-3).

Not applicable.

Not applicable.

Not applicable.

No potential conflict of interest relevant to this article was reported.