A transcriptomic analysis of skeletal muscle tissues reveals promising candidate genes and pathways accountable for different daily weight gain in Hanwoo cattle

Cattle traits like average daily weight gain (ADG) greatly impact profitability. Selecting based on ADG considering genetic variability can lead to economic and genetic advancements in cattle breeding. This study aimed to unravel genetic influences on ADG variation in Hanwoo cattle at the skeletal muscle transcriptomic level. RNA sequencing was conducted on longissimus dorsi (LD), semimembranosus (SB), and psoas major (PM) muscles of 14 steers assigned to same feed, grouped by low (≤ 0.71 kg) and high (≥ 0.77 kg) ADG. At P ≤ 0.05 and log2fold > 1.5, the distinct pattern of gene expression was identified with 184, 172, and 210 differentially expressed genes in LD, SB, and PM muscles, respectively. Tissue-specific responses to ADG variation were evident, with myogenesis and differentiation associated JAK-STAT signaling pathway and prolactin signaling pathways enriched in LD and SB muscles, while adipogenesis-related PPAR signaling pathways were enriched in PM muscle. Key hub genes (AXIN2, CDKN1A, MYC, PTGS2, FZD5, SPP1) were upregulated and functionally significant in muscle growth and differentiation. Notably, DPP6, CDKN1A, and FZD5 emerged as possible candidate genes linked to ADG variation. These findings enhance our understanding of genetic factors behind ADG variation in Hanwoo cattle, illuminating skeletal muscle mechanisms influencing ADG.

on small-scale experimental cohorts, with the objective of elucidating the genetic components involved in the intricate biological processes associated with complex traits and diseases [5][6][7][8] .RNA-Sequencing (RNA Seq) has emerged as a preferred method over previous techniques like microarrays due to its ability to capture the entire transcriptome, including novel transcripts, alternative splicing events, and non-coding RNAs.Moreover, it has been suggested that the patterns of gene expression show higher levels of concordance across various breeds and populations, when compared to other analysis for example genome-wide association studies (GWAS) 5,6 .
Most earlier studies employed GWAS-based methods to analyze the genetic basis of carcass traits, revealing several genomic regions containing QTLs associated with carcass weight in diverse cattle breed 7,9,10 .However, the investigation of carcass traits in Hanwoo cattle through transcriptomic studies, employing RNA sequencing for functional evaluation of pivotal candidate genes, is still in a dynamic phase and requires more comprehensive analyses.
In the present study, we conducted a transcriptome analysis on three specific skeletal muscles-longissimus dorsi muscle (LD), semimembranosus (SB), and psoas major (PM).These samples were collected from 14 steers, divided into two groups based on their ADG, in order to gain a comprehensive understanding of the genetic mechanisms underlying this complex trait and explore potential candidate genes or molecules involved.We performed functional enrichment and pathway analyses on the differentially expressed genes (DEGs) to identify significantly enriched gene ontology (GO) terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways associated with ADG.Additionally, we constructed a protein-protein interaction (PPI) network to identify hub genes based on their connectivity within the network.

Results
A cohort of 14 Hanwoo steers that were divided into two ADG groups-high ADG (≥ 0.77 kg, 7 steers) and low ADG (≤ 0.71 kg, 7 steers) were the subject of a transcriptome investigation.A full summary of our study is summarized in Fig. 1.The significant difference (P < 0.001) of ADG between two group was shown in Supplementary Fig. 1.The details regarding the animals, including their phenotype, age, and sex, can be found in Table 1 of the study.

Transcriptome analysis
In this study, RNA was extracted from specific three skeletal muscle tissue samples, including LD, SB, and PM muscle tissues.The extracted RNA was then used to construct libraries for deep sequencing.After sequencing, the obtained sequences were utilized for further analysis.Finally, we determined the total DEGs genes in each tissue, as presented in Supplementary File S1.In order to detect any underlying patterns in the samples, a principle component analysis (PCA) was carried out on DEGs without applying batch correction.The analysis effectively separated the samples of high ADG group from those the low ADG group samples (Fig. 2A-C).The results demonstrated that the first principal component (PC1) explained 30% of the variability, while PC2 explained 19% of the variability within the DEGs of LD samples, as shown in Fig. 2. For SB tissue, PC1 explained 21% of the variability, whereas PC2 explained 17% of the variability.For PM tissue, it was 25% of the variability of PC1, and 20% variability for PC2.
The expression pattern of DEGs across samples was analyzed using hierarchical clustering.The results revealed that the gene expression of the low ADG group was significantly different from that of the samples of high ADG group in all the tissues.Heatmaps were generated to visualize the expression patterns of the top 500 DEGs in each tissue of two different ADG group (Fig. 2D-F).To further investigate the differential gene expression, we generated a volcano plot (Fig. 3A).Venn diagrams were used to illustrate the overlap of upregulated and downregulated DEGs between the evaluated tissues (Fig. 3B-D).The Venn diagram shows only six genes were found to be overlapped.The correlation between the expression levels of tissues were shown in Supplementary Fig. 2.

KEGG pathway and network analysis
In this study, the KEGG pathway analysis successfully annotated the DEGs to a total of 16, 14, and 10 pathways for LD, SB, PM tissues (Fig. 4).KEGG pathway analysis of predicted DEGs targets featured enrichment for bta04630:JAK-STAT signaling pathway, bta05224:breast cancer, bta05226:gastric cancer, bta04978:mineral absorption, and bta05200:pathways in cancer in LD tissue for the SB tissue, we found the bta04917:prolactin signaling pathway, bta04625:c-type lectin receptor signaling pathway, bta05145:toxoplasmosis, and bta05142:chagas disease pathways to be significantly enriched.KEGG pathway analysis of predicted DEGs of PM tissue targets featured enrichment for bta04920:adipocytokine signaling pathway, bta03320:PPAR signaling pathway, bta04923:regulation of lipolysis in adipocytes, bta04060:cytokine-cytokine receptor interaction, bta04152:AMPK signaling pathway, bta04978:mineral absorption, bta04310:WNT signaling pathway.Figure 4 showcased that genes such as IL4, FZD5, FRAT2, CDKN1A, FGF16 in LD (Fig. 4A), IL12B and MAPK13 in SB (Fig. 4B), and PLIN1, FABP4, PCK2, ADIPOQ in PM (Fig. 4C) tissue were the one which are involved in multiple pathways regulation.The expression of genes involved in pathways for each tissue is shown in heatmap (Fig. 4D-F).The KEGG graph of top pathways found in each tissue were given in Supplementary Fig. 3.This analysis helps to identify the specific pathways or signaling networks that are potentially dysregulated in relation to the observed gene expression changes.www.nature.com/scientificreports/

GO enrichment and network analysis
To evaluate the functional relationships and biological changes associated with the identified DEGs, gene enrichment analysis using predetermined gene sets and gene ranks was conducted.Specifically, GO and KEGG pathway enrichment analyses were performed on the set of 186 (LD), 175 (SB) and 211 (PM) up-regulated and downregulated DEGs.The GO enrichment analysis revealed a total of 7 significant GO terms, with four terms significantly enriched in the biological process category, two terms significantly enriched in the cellular component category, and one terms significantly enriched in the molecular function category for LD tissue (Fig. 5A).For SB tissue, a total of 12 significant GO terms was identified in GO enrichment analysis.Among these, two terms showed significant enrichment in the biological process category, six terms exhibited significant enrichment in the cellular component category, and four terms displayed significant enrichment in the molecular function category (Fig. 5B).In PM tissue, a total of 25 significant GO terms found where 15 terms significantly enriched in the biological process category, seven terms significantly enriched in the cellular component category, and three terms significantly enriched in the molecular function category for PM tissue (Fig. 5C).These enriched GO terms provide insights into the functional roles and processes that the DEGs may be involved in.Further, all the significant GO terms were shown in Fig. 2. Here, we stated the top significant GO terms that were GO:0030332 ~ cyclin binding, GO:0016020 ~ membrane, GO:0005615 ~ extracellular space found to be associated with ADG variation in LD, SB, and PM tissue, respectively.In addition, we have constructed a network of the ten top mostly enriched GO terms for each tissue (Fig. 6).www.nature.com/scientificreports/

Hub genes, co-expression module and pathway detection by PPI network and module analysis of DEG
The PPI network was constructed using the STRING database, resulting in a network consisting of 363 nodes and 902 edges.Supplementary Fig. 4 provide visual representations of the network structure.The hub genes were selected for further investigation using PPI analysis.In order to identify hub genes for PPI analysis, all four methods within the CytoHubba plugin were utilized.The top 20 genes resulting from each method were compiled and presented in Table 3.Among these lists, a total of 12 genes that appeared in all four methods were considered as hub genes (JUN, CDKN1A, CDH1, TP53, FZD5, PPARG , LEP, IL4, PTGS2, AXIN2, CDH2, and MYC).These findings highlight the centrality of these hub genes and suggest their involvement in critical biological processes or pathways associated with growth performance in Hanwoo.
After conducting module analysis using the MCODE plugin in Cytoscape, a total of nineteen clusters were obtained.Based on their MCODE scores, the top three modules were selected as hub modules (Fig. 7).Among these modules, module 1 had the highest MCODE score (8) and consisted of 12 nodes and 23 edges.Module 2 was comprising 11 nodes and 24 edges with MCODE score 4.2.Next, the Module 3 had MCODE score of 4.182 and contained 9 nodes and 14 edges.Importantly, all genes within module 1 showed up-regulation, indicating their potential functional significance.
To gain further insights into the biological functions of the genes within each module, we performed a biological functional enrichment analysis using the DAVID database (Table 4).The genes in Module 1 showed significant enrichment in a total of 35 significant pathways.Here, we have highlighted here some of the top significant pathways such as PI3K-Akt signaling pathway, WNT signaling pathway, JAK-STAT signaling pathway, lipid and atherosclerosis, and MAPK signaling pathway.Module 2 including 11 genes was found to be associated with different biological pathways.Among them, tissue growth associated pathways were osteoclast differentiation, regulation of lipolysis in adipocytes, PPAR signaling pathway, IL-17 signaling pathway, and C-type lectin receptor signaling pathway.Module 3 consisting of 9 genes was noticed to be linked to important pathways including PI3K-Akt signaling pathway and Wnt signaling pathway.In addition, the expression level of genes found in all three modules are presented in Table 5.

Discussion
Controlling factors that contribute to ADG is crucial for the Hanwoo beef industry, as a high occurrence of low weight can result in significant economic losses.Therefore, a comprehensive understanding of the genetic factors underlying key growth and economic traits of livestock is essential for identifying new genes and genetic pathways that can be leveraged through genomic selection to improve economic traits in livestock breeding programs.
Hence, the current study was carried out transcriptome analysis on these three LD, SB, and PM skeletal muscles to provide the novel molecular functions of the tissues, which may influence the ADG in Hanwoo steer.In South Korea, LD, the PM, and SB holds significant importance due to the high demand for steak.In Hanwoo beef, the LD and PM skeletal muscles, commonly referred to as top fillet and under fillet in the meat industry, stand out as highly coveted cuts of the carcass, possessing substantial commercial and culinary significance 11,12 .In our study, we selected male Hanwoo cattle with a significance difference in ADG group (high ADG ≥ 0.77 kg and low ADG ≤ 0.71 kg).The result of this study revealed that DPP6 (LD), IL22RA1 (SB), and U6 (PM) genes were the top upregulated genes.Among these genes, DPP6, also known as dipeptidyl peptidase-like protein 6, is a gene that encodes a protein involved in neuronal signaling and ion channel regulation 13 .It has been suggested as candidate associated with altered lipid profile 14 , muscle atrophy in human 13 , and the rib eye Muscle Area in Hu Sheep 13 .According to Zhao et al., it is plausible that DPP6 may have some influence on muscle physiology and potentially contribute to muscle wasting or atrophy 13 .IL22RA1 (interleukin 22 receptor subunit alpha 1) has been shown to be receptor for interleukin-22 (IL-22), a cytokine involved in various biological processes, including cell proliferation and tissue regeneration in animals and human 15 .Therefore, DPP6 gene can be suggested as one of good possible candidate gene for ADG variation in Hanwoo.While IL22RA1's role in muscle growth is not extensively studied, there is evidence suggesting its secondary involvement in muscle physiology [15][16][17] .IL22RA1 being identified as a candidate gene associated with health, adaptation and reproduction traits in cattle 18 .For PM tissue, small nuclear RNAs (U6) plays catalytic role at the core of the spliceosome 19 but straight function associated to muscle development or growth performance in cattle not been reported before.Furthermore, a number of genes that were downregulated involved in diverse function linked to cell death, metabolism, cell growth and development were screened as possible candidate genes for ADG trait.Of the top downregulated DEGs, myosin light chain 6B (MYL6B), fibroblast growth factor 6 (FGF6), acyl-CoA synthetase bubblegum family member 1(ACSBG1), two pore channel 3(TPC3), polymeric immunoglobulin receptor(PIGR), Wnt family member 16 (WNT16), alpha-2-glycoprotein 1, zinc-binding (AZGP1), and forkhead box Q1(FOXQ1) genes were found to be involve in various function associated to growth, development and metabolism.
The total upregulated and downregulated significant DEGs were used to determine the GO and KEGG pathway analysis, which was then used to identify significant enriched pathways for the three tissues.Based on the pathway analysis, JAK-STAT signaling pathway was detected as the most significant enriched pathway in LD tissue.In muscle development and regeneration, this pathway plays a vital role in regulation of adipogenesis, process of fat cell development, controlling myoblast proliferation, differentiation, and the formation of muscle fibers [20][21][22] .A previous GWAS study identified candidate genes associated with the body weight trait in Chinese sheep, revealing their involvement in the JAK-STAT pathway 23 .
Among the upregulated DEGs in JAK-STAT pathway, we observed a significant increase in the expression of CDKN1A (cyclin dependent kinase inhibitor 1A) gene and common enrichment in conjunction with other pathways in our result.This gene is associated with the regulation of the G1/S checkpoint and is known to be crucial in cellular senescence during preimplantation embryo development in humans 24 .Zhang et al. reported that CDKN1A regulates myogenesis process in muscle stem cells 25 .Recently, CDKN1A gene has been reported as strong candidate gene for high fat milk production trait in dairy cattle (Chinese Holstein population) [26][27][28] .Moreover, another vital gene, SOCS3 (suppressor of cytokine signaling family 3) known as negative regulator of JAK-STAT, was detected to be downregulated in our study 29 .It indicates that the downregulation of SOCS3 could leads to STAT3 phosphorylation which further may promotes lipid metabolism, cell cycle progression, differentiation 30,31 .Although future studies are expected to provide more comprehensive insights into the specific In SB tissue, the prolactin signaling pathway emerged as the most significant pathway.Here, our result revealed a significant upregulation of prolactin signaling pathway associated PRLR (prolactin receptor) gene by a log2FC of 3.05.The prolactin signaling pathway has significant involvement in the regulation of metabolism, immune system function, milk synthesis, cell growth, survival and pancreatic development.Activation of this signaling pathway was reported to promote cell growth and survival 32 .According to Chu et al., the PRLR gene closely related to prolactin signaling pathway is a candidate gene for prolificacy of small tail trait in han sheep 33 .In Egyptian Buffaloes, the PRLR has shown promise as a genetic marker for evaluating milk production and quality traits 34 .Additionally, another gene associated with the prolactin signaling pathway, MAPK13 (mitogenactivated protein (MAP) kinase family), exhibited downregulation.MAPK13 encodes a protein belonging to the mitogen-activated protein (MAP) kinase family, which serves as a convergence point for various biochemical signals 35 .MAP kinases are involved in numerous cellular processes including proliferation, differentiation, transcription regulation, and development 35,36 .Therefore, it can possibly suggest that PRLR and MAPK13 gene can be important possible biomarker associated to high ADG in Hanwoo.
Subsequently, the KEGG analysis identified PPAR signaling pathway as a potential candidate pathway that was the top most significant pathway found from the transcriptome analysis of PM tissue.Further, we identified the major downstream pathways associated with PPARs signaling pathway such as adipocytokine signalling pathway www.nature.com/scientificreports/and regulation of lipolysis in adipocytes.The PPAR signaling pathway has been examined in many studies with regard to its function in adipose tissue influences adipogenesis, glucose metabolism, and adipokine secretion in cattle 8,37 .Recently, He et al., described that 13 tag single nucleotide polymorphisms in the PPAR signaling pathway are associated with porcine meat quality traits 38 .A study on Hanwoo cattle investigated the gene expression patterns associated to PPAR signaling pathway.The results of this study demonstrated a significant association between PPAR signaling and genes involved in fatty acid oxidation.This association ultimately led to an increase in triglyceride formation through ATP production 8 .A muscle transcriptome study conducted earlier identified the adipocytokine signaling pathway and PPAR signaling pathway as potential candidate pathways associated with the longissimus thoracis muscles in meat-type sheep from India 39 .Hence, it can be hypothesized that the PPAR signaling pathway plays a role in controlling adipogenesis and influencing the capacity of adipose tissue to  www.nature.com/scientificreports/store lipids, thereby affecting the fat mass in Hanwoo muscles.Consequently, by regulating the size and quantity of adipocytes within the muscle, this pathway may impact the average daily gain characteristic in Hanwoo cattle.In our study, the genes involved in PPAR signaling pathway were ADIPOQ (Adiponectin, C1Q and collagen domain containing), PCK2 (Phosphoenolpyruvate carboxykinase 1), FABP4 (Fatty acid binding protein 4, adipocyte), PLIN1 (perilipin).FABP4 gene reported to be facilitates the uptake of fatty acids and their intracellular trafficking, contributing to the process of adipocyte maturation and lipid accumulation 8,40 .ADIPOQ is a well-known homeostatic factor for regulating glucose levels, and lipid metabolism and reported to be candidate gene for carcass trait 41 .PLIN1 plays a crucial role in controlling the levels of triglycerides and the size of lipid droplets in adipocytes, contributing to the maintenance of lipid balance 42 .It is also suggested as possible robust candidate gene marker for body weight in cattle breeding programs 43 .Furthermore, the SLC2A1 gene enriched in adipocytokine signaling pathway was upregulated by more than twofold here in PM tissue.SLC2A1, a glucose transporter, has also been implicated in facilitating the transport of fatty acids 44 .In summary, it is postulated that the mentioned genes may facilitate the accumulation of intramuscular fat in cattle by controlling the synthesis, transport of fatty acids, lipid balance by regulating the respective signaling pathways and potentially affecting the ADG trait.However, a more comprehensive understanding of the precise regulatory mechanisms involved requires further exploration in future studies.
Our GO analysis revealed that DEGs of LD tissues were mainly enriched in branching involved in labyrinthine layer morphogenesis, retinoic acid metabolic process, T cell chemotaxis and detoxification of copper ion biological process.For the SB tissue, we found significant biological process GO terms that are the sodium ion transmembrane transport and response to estrogen biological process.In PM tissue GO analysis, the result revealed some important biological process GO terms that are interleukin production, fat cell differentiation, carbohydrate metabolic process, ion transmembrane transport, insulin resistance and etc.Among all these GO terms, ion transmembrane transport facilitates nutrient absorption, carbohydrate metabolism provides energy, and fat cell differentiation contributes to energy storage, all of which are essential processes for weight gain in cattle.These biological processes have been reported to be have dominance effect on porcine and cattle final weight and back fat thickness 18,45,46 .Moreover, the GO network analysis revealed that several essential genes, which were commonly shared among significant GO terms, exhibited significant upregulation.For example, FZD5 (frizzled class receptor 5) involved in two different GO terms are found to be upregulated in our datasets.These genes likely to play crucial roles in the biological processes, cellular components, or molecular functions represented by the enriched GO terms.Further investigation and validation of these genes can provide insights into their specific roles and mechanisms for ADG in Hanwoo cattle.www.nature.com/scientificreports/Furthermore, we constructed a PPI network to gather valuable interaction information regarding the differentially expressed genes (DEGs).The PPI network in the current study was created with total significant DEGs (P-value of ≤ 0.05, log2FC of ≥ 1.5) collected from three experimental tissues.In present investigation, we identified 12 hub genes that are JUN (transcription factor Jun), CDKN1A, CDH1 (cadherin-1), TP53 (tumor Protein P53), FZD5, PPARG (Peroxisome proliferator-activated receptor gamma), LEP (leptin), IL4 (interleukin 4), PTGS2 (prostaglandin-endoperoxide synthase 2), AXIN2 (axin-related protein 2), CDH2 (cadherin-2) and MYC (MYC proto-oncogene).These findings highlight the centrality of these hub genes and suggest their involvement in critical biological processes or pathways associated with growth performance in Hanwoo.By integrating the DEGs from these tissues, the PPI network was generated to capture correlation patterns of gene expression across samples, the new interactions and relationships between proteins that tend to be co-regulated or functionally related 47 .
Among all the hub genes, notably the CDKN1A and FZD5 were spotted to be upregulated (> onefold) in all the three experimental tissue samples.In our study, we have identified CDKN1A as a key hub gene that works as inhibitor in cell cycle and control myogenesis process in muscle stem cells 25 .Furthermore, our module analysis indicated that CDKN1A is associated with module 1.Interestingly, in our current experiment, we observed that CDKN1A is enriched in various pathways, including those related to different types of cancers, Hepatitis B, Epstein-Barr virus infection, PI3K-Akt signaling pathway, p53 signaling pathway, and JAK-STAT signaling.CDKN1A has been described as candidate gene associated with carcass quality trait and milk production in beef cattle 48 .These findings along with our pathway analysis result collectively suggest that CDKN1A may influence muscle growth in Hanwoo cattle by regulating the muscle differentiation at the myogenin step, thereby subsequently affecting the ADG trait.However, more research is needed to confirm this theory.The other gene, FZD5 Table 5.Expression level of genes found in three modules in all the tissues.Up-regulated genes are highlighted in bold and downregulated genes in normal typeface.FC fold change, LD longissimus dorsi muscle, SB semimembranosus, PM psoas major.# The upregulated module genes identified as hub gene.www.nature.com/scientificreports/ is a protein coding gene, which functions as transmembrane receptor and mediates Wnt ligands binding in Wnt signaling pathway and plays role in the control of tissue regeneration, cell proliferation and cell differentiation 49,50 .
One recent study has proposed that the activation of FZD5 by Wnt protein could potentially play a role in regulating adipose differentiation 49 .A previous GWAS study identified FZD5 as a candidate gene associated with fertility traits in Beef Heifers 51 .They suggested that the FZD5 gene plays a role in the WNT pathway, which aligns with our findings and provides additional support for our data 51 .Therefore, it could be valuable to study further on FZD5 functions in the muscle growth and ADG trait.
In current study, the top three modules were selected for further investigations of PPI network.Our result presented that the genes in all the modules are involved in several important signaling pathways including PI3K-Akt signaling pathway, Wnt signaling pathway, JAK-STAT signaling pathway, PPAR signaling pathway, IL-17 signaling pathway, and C-type lectin receptor signaling pathway.PI3K-Akt signaling and Wnt signaling pathway were found to be common in two modules i.e. module 1 and module 3. A considerable number of studies have reported that the PI3K-Akt signaling and Wnt signaling pathways are the key pathways for regulating cell propagation, adipose proliferation, fat deposition, and muscle development in bovine thus subsequently effect the meat production and carcass traits in cattle [52][53][54][55] .Notably, the outcomes from GO/KEGG analysis were matched to the functional annotations of genes in the most top three modules of PPI network.Additionally, it is noteworthy to mention that genes enriched in module 1 (CDKN1A, AXIN2, CDH2, MYC, and TP53), module 2 (JUN, PPARG, PTGS2), and module 3 (TP53, FZD5, IL4, and MYC) were found as hub genes as well.The speculation that the identified genes, which are differentially regulated and associated with various pathways in different tissues, could be potential candidate genes affecting the ADG trait is reasonable.However, it is crucial to note that further intensive investigations are necessary to validate and confirm the significance of the identified candidate genes and pathways.Some limitations are associated with this study.Firstly, the sample size was relatively small and the age difference between low and high group reaches up to 168 days.Additionally, our research was exclusively conducted on only three skeletal muscles.One advantage of this approach is the ability to obtain results from cattle with consistent or uniform characteristics.This study is also the first to report the use of three different skeletal muscles for ADG trait in Hanwoo cattle.A larger scale transcriptome study is needed on different traditional adipose depots subcutaneous (SQF), visceral (VAT), intermuscular (INTMF), intramuscular (IM), and bone adipose for better understanding of molecular mechanism of ADG in cattle 56 .

Conclusion
To sum up, our study involved a comparative analysis of three different skeletal muscle transcriptomes in Hanwoo cattle that exhibited varying levels of average daily weight gain.A total of 200, 172, and 210 DEGs were screened from LD, SB and PM muscle tissues, respectively.Through functional enrichment, hub gene, and module analysis of DEGs in three distinct muscle tissues, we have identified several possible candidate genes, namely DPP6, CDKN1A, and FZD5, along with significant pathways, including the JAK-STAT signaling pathway, prolactin signaling pathway, and PPAR signaling pathway.These findings suggest a crucial involvement of these candidate genes and pathways in influencing the variation of average daily gain (ADG) in Hanwoo cattle.To fully comprehend their specific mechanisms related to carcass weight traits, further investigations are warranted.

Animals, phenotypes, and sample collection
The National Institute of Animal Science in the Republic of Korea's Animal Care and Use Committee gave its approval to all experimental techniques mentioned in this paper (Approval No. NIAS-2022133).The National Institute of Animal Science's (NIAS) Animal Genomics and Bioinformatics Division assembled the tissue samples of Hanwoo cattle for the dataset used in this study.Also, all the experiments in the manuscript follows the recommendations in the Animal Research Reporting in Vivo Experiments (ARRIVE) guidelines.At the slaughterhouse of the National Institute of Animal Science, the cow was stunned with a captive bolt gun and then slaughtered as per American Veterinary Medical Association (AVMA) Guidelines for the Euthanasia of Animals.The average daily weight gain (ADG, kg/d) was examined and calculated the ratio as following-In this study, 14 cattle used were at growing stage which corresponds to the early fattening stage.Additionally, these cattle were assigned to the same feed and all the cattle were weaned at the same year.Post-slaughter muscle samples including the longissimus dorsi muscle (LD), rump semimembranosus (SB), and psoas major (PM) were collected.Sterile surgical techniques were used, such as washing the muscle collection site with 70% ethanol and using a sterile scalpel to make skin incision 7 .The sample integrity was preserved through immediate preservation in liquid nitrogen and subsequent storage at − 80 °C.

RNA isolation and sequencing
The total RNA was isolated from each homogenized tissue sample (LD, SB, and PM) using the Qiagen RNeasy kit (Qiagen, Hilden, Germany) following the manufacturer's instructions.The RNA's purity and integrity were assessed using a nano photometer (IMPLEN, CA, USA) and the RNA Nano 6000 Assay kit on the Bioanalyzer 2100 equipment (Agilent Technologies, CA, USA).Samples with RNA integrity number (RIN) values exceeding 8 were selected for library preparation.The NEBNext UltraTM RNA Library Prep Kit for Illumina (NEB, Iswich, MA, USA) was employed to prepare the library, involving random fragmentation, adapter ligation, and tagmentation as per the manufacturer's guidelines.The sequencing was performed using the Illumina NovaSeq 6000 platform.The adapter-ligated fragments were amplified via PCR, followed by gel purification.Cluster generation was achieved using the TruSeq PE cluster kit v4-cBot-HS (Illumina) following the manufacturer's instructions.For sequencing, the library was loaded into the Illumina NovaSeq 6000 System, which produced paired-end reads.Base calling on the raw images generated by the Illumina sequencer was performed using real-time analysis (RTA) software.Finally, the generated BCL image files were converted to FASTQ raw readings using the Illumina BCL to fastq program, which were then used for post-analysis.

RNA sequencing data processing and analysis
For the data analysis in this study, the input was derived from FASTQ files.The quality of the raw RNA-seq data was assessed using FastQC v0.11.Trimmomatic v0.39 was utilized to perform trimming of N bases, removal of adapters from the 5' and 3' ends, and filtering of reads with a quality score below Q20 for quality control during the generation of raw reads with the paired-end option.The filtered sequences were mapped to the Bos taurus reference genome (version ARS-UCD 1.2) using the Hisat2 v2.2.1 mapping program with the corresponding index file.Subsequently, the mapped reads were used to generate a count matrix using FeatureCounts from the Subread package (v2.0.1).A PCA plot without any batch correction was then generated to visualize the data.Finally, the differential expression analysis of genes (DEGs) was performed in R (v3.4.4) using each read counts matrix file.

Figure 1 .
Figure 1.The diagram illustrates the experimental design.DEGs differentially expressed genes, GO Gene Ontology, KEGG Kyoto Encyclopedia of Genes and Genomes, PPI protein-protein interaction.

Figure 2 .
Figure 2. (A) Data preprocessing and differential expression gene (DEG) analysis on the LD, SB, and PM muscle tissue samples from low and high average daily weight gain Hanwoo cattle.Principal component analysis indicating the overall profiles of three muscles datasets (A) LD, (B) SB, and (C) PM.Heatmap of the top 500 DEGs of (D) LD, (E) SB, and (F) PM muscles at P value < 0.05 and log2FC > 1.5 and column represents the gene expression pattern of each sample.Up-and down-regulated transcripts are shown in red and blue color reFigurespectively.LADG low average daily weight gain, HADG high average daily weight gain, LD longissimus dorsi, SB semimembranosus, PM psoas major.

Figure 3 .
Figure 3. (A) Volcano plots visualizing significant DEGs with P ≤ 0.05 and log2FC ≥ 1.5 between low and high average daily weight gain (ADG) Hanwoo cattle muscle tissue samples.Red points symbolize up-regulation, while blue points designate down-regulation; gray points signify non-significant expression.The number of up-and down-regulated genes presented as red arrow and blue arrow symbol respectively.A Venn plot to identify common (B) up-regulated, (C) down-regulated, (D) total genes among the muscle tissue samples.LD longissimus dorsi, SB semimembranosus, PM psoas major, FC fold change.

Figure 4 .
Figure 4. Function enrichment analysis of differently expressed genes.Sankey plot representing the each enriched KEGG pathway (P ≤ 0.05) for (A) LD, (B) SB, (C) PM tissues.Heatmap represents the expression of each gene associated with respectively enriched pathway in (D) LD, (E) SB, (F) PM muscle tissues.LD longissimus dorsi, SB semimembranosus, PM psoas major, KEGG Kyoto Encyclopedia of Genes and Genomes.

Figure 6 .
Figure 6.The network analysis of gene ontology of (A) LD, (B) SB, (C) PM muscle tissue samples showing genes involved in different gene ontology terms.LD longissimus dorsi, SB semimembranosus, PM psoas major.

Figure 7 .
Figure 7. Module analysis from protein-protein interaction (PPI) network.Top three significant modules including (A) module 1, (B) module 2, and (C) module 3 were obtained from PPI network of differentially expressed genes with MCODE.

Figure 8 .
Figure 8. Verification of expression level of selected DEGs in the (A) LD, (B) SB, (C) PM muscle tissue samples by quantitative reverse-transcription-PCR (qRT-PCR).Ten genes randomly selected from both up-regulated and down-regulated of each tissue were analyzed.The fold changes obtained from qRT-PCR was compared with the fold changes of the RNA-seq result.DEG differentially expressed gene, LD longissimus dorsi, SB semimembranosus, PM psoas major.

( 1 ) 9 Table 6 .
ADG kg/d = Body weight kg at Slaughter day − Initial Born weight (kg) Age days Vol:.(1234567890)Scientific Reports | (2024) 14:315 | https://doi.org/10.1038/s41598-023-51037-List of primer used in this study for qRT-PCR verification of RNA-seq results.LD longissimus dorsi muscle, SB semimembranosus muscle, PM psoas major muscle.ACC TAC CTC CT TCG GTT GTC AAG TCC ATG GA Up-genes PIPOX CTG GCA AGA GAA GGT TCC TG CAT TGT TGC CGT GGT GAT AG FRMPD3 GAG GAG GAC GTG AGT GAA GC ATC GGA CCT TCA CAG GAT TG BHMT2 TCT GGA TAG TGG GGA GGT TG TAC AGC TTC CCA CAG GCT TT ERICH6 GCG GAG AAA ACA AGA ACG AC TCC CGT CTG GAA ATG AAG TC ZDHHC22 CTG GGC AAC TAT GTC CTG GT GCG AGG TGT AGA GGC AGA AC Down-genes DRC1 GGA TGG GAA GTC ACA GGA GA TGG CAG GGA TAA CCG TAG TC CDKL4 GGC AAA CAC TCC AAG CTC TC TAT CCC CCA CAA GAA GTT CG MFRP GTG CCA AGT TTT CAG GGT GT GGC TGA AGT TGT GGA AGA GC PLPPR3 SB tissue IL22RA1 CCA CAG ACA AGG GTC CAA GT CCT CGG AAC AGA GAG TCC AG Up-genes CLDN7 TGC CTT GAT AGC TTG CTC CT CGG TAT CCA GCT TTG CTC TC EHF ATG AGT CTG CAG GAG TTC AC GCT ACC ATA GTT GGT GTC GT SLC34A2 AGC CTG AGA AAG CTA AGG AG CTC TCT CTG ACC ACT TGA GC PITX1 CAG CTC CAT CTC CTC CAT AGT GCT GCT TAG ACT TGA GC KCNG2 GAG TGC TCC CCT AAG TGT CG ACA CGT AGA AGG GCA GGA TG Down-genes STMN2 ATC TGC TCT TGC TTC TAC CC GCC TCC TGA GAC TTT CTT CT OTOGL CAC AGG CTT CCA CAC TTT GA CAT TGC CCC TTC ATC AGT CT MYL6B TCC AGG GTT GAG CTA CCA TC TAT CCC AAG ACC CTG AGC AC FGF6 CCC AGC TTC CAA GAA GAG TG CTA GGG AGG AAG TGC GTG AC PM tissue SRARP GAA CCT GGA GAC CAG CTC AG GAG CCT CCA GAG TCA GGT TG Up-genes CSTB CAG CTG GAG GAG AAG GAG AA GGT CTG GTA GCT GGT CAA GG PADI4 CCA GGA GTG GTT GTG AAG GT GCA CTC AGG GAG ATT TCG AC ISG20 ACT ATA GGA CCC CGG TCA GC GGC GTA GTC GCT CAT GTT CT SPP1 GAT GGC CGA GGT GAT AGT GT GTT GGT TTC CTT GCT GTG GT MPTX1 CCT GAA AGC CTT CAC AGA CC TCC CAG CTG ACA CAG AGA TG Down-genes MGAM TGC CTT GAC ATA CCG TAC CA CCA TAG CGA CAC AGC TGA AA PIGR ACC ATC AAC TGC CCT TTC AC TGA GCT TGA CTC GGT TGA TG WNT16 CAT GAC CGA GTG TTC CTG TG CTG CCT TCC AGC TTC ATT GT GNLY AGG AGA AGA GCT GGG CCT AC TTT CCA GCT GTG ATG TCT GC

Table 1 .
Summary of phenotypic data of animals used in this study (n = 14).ADG Average daily weight gain.***p < 0.001 indicate significant difference of ADG between Low-ADG and High-ADG group.

Table 3 .
The top 20 hub genes identified from protein-protein interaction network using the cytohubba plugin in cytoscape.Degree node connected degree, EPC edge percolated component, MCC maximal clique centrality, MNC maximum neighborhood component.

Table 4 .
The biological pathway enrichment analysis of genes from the three modules of protein-protein interaction network.