Identification of Biological Pathways Contributing to Marbling in Skeletal Muscle to Improve Beef Cattle Breeding

Red meat is an important dietary source that provides part of the nutritional requirements. Intramuscular fat, known as marbling, is located throughout skeletal muscle. Marbling is a trait of major economic relevance that positively influences sensory quality aspects. The aim of the present study was to identify and better understand biological pathways defining marbling in beef cattle. Pathway analysis was performed in PathVisio with publicly available transcriptomic data from semitendinosus muscle of well-marbled and lean-marbled beef. Moreover, for Bos taurus we created a gene identifier mapping database with bridgeDb and a pathway collection in WikiPathways. The regulation of marbling is possibly the result of the interplay between signaling pathways in muscle, fat, and intramuscular connective tissue. Pathway analysis revealed 17 pathways that were significantly different between well-marbled and lean-marbled beef. The MAPK signaling pathway was enriched, and the signaling pathways that play a role in tissue development were also affected. Interestingly, pathways related to immune response and insulin signaling were enriched.

Introduction

Red meat is as an important dietary source that provides part of the nutritional requirements such as proteins, minerals, B-complex vitamins, and essential fatty acids (McAfee et al., 2010). Control of meat quality is very important for meat producers and meat sellers to satisfy customer’s preferences (Bernard et al., 2007). Marbling, a trait that describes the presence intramuscular fat, is of major economic relevance for beef producing cattle that has a positive impact on sensory quality traits, such as flavor, juiciness, and tenderness of meat. Studies have shown that marbling depends on factors such as breed, genotype, age, diet, husbandry, and growth stages. Although in marbling the environmental factors play an important role, the genetic background of the animals is the major factor defining the marbling status (Yamada et al., 2006). O’Connor et al. studied the effect of breed-type on marbling, their results demonstrated that increase in meat marbling from Bos taurus cattle (Hereford, Red Angus, Angus, and Tarentaise breeds) can increase the tenderness, more than the Bos indicus cattle (Braford, Red Brangus, and Simbrah breeds) (O’Connor et al., 1997). Although, Shackelford et al., found that, the lower tenderness of meat from Bos indicus cattle is mainly because of decreased postmortem proteolysis which result from elevated calpastatin activity that the possibility existed for an interaction between breed and the influences of marbling score on tenderness (Shackelford et al., 1991). Also, Wulf et al., reported that, in Charolais and Limousin breeds, marbling correlated with calpastatin activity and shear force. He suggested that selection for increased marbling based on these genetic effects, in those two breeds, might be effective for enhancing beef tenderness (Wulf et al., 1996). There is substantial evidence from transcriptomics studies that gene expression profiles affect phenotypic variation for marbling (Cesar et al., 2015). Understanding the signaling pathways that make up the regulatory network in the marbling process can help steer the breeding process (Thaller et al., 2003). Therefore, animal breeding specialists have attempted system-oriented approaches to investigate major economic traits (Lee et al., 2010). It has been shown that marbling differences may be a function of a number of complex interactions among biological pathways. Therefore, a pathway analysis with differential gene expression patterns can result in a better understanding of muscle physiological states and their influence on beef quality and animal welfare (Hocquette et al., 2012). The public availability of transcriptomics data from beef producing cattle, provides new opportunities to explore the global gene expression in muscle to investigate physiological processes and their influence on meat sensory quality traits (Lee et al., 2010; Hocquette et al., 2012).

Within the genomic region of marbling there are several genes considered as parts of QTLs such as EDGPR1, Titin, Akirin 2, and RPL27 (Takasuga et al., 2007) which were mapped in a half-sib family of Japanese Black cattle (Yamada et al., 2006). Thus, these genes were considered as positional functional candidates for the genes responsible for marbling. This study aims at identifying genes and biological pathways regulating marbling of muscle tissue in beef cattle based on publicly available transcriptomics data obtained from a study by Sadkowski and coworkers (2014). We updated and extended the pathway collection for B. taurus at WikiPathways (Slenter et al., 2017) an online pathway repository, and a B. taurus gene product identifier mapping BridgeDb database was created to allow mapping of expression data to the gene databases identifiers used in the pathways (van Iersel et al., 2010). Sadkowski et al., 2014 measured global gene expression in skeletal muscle of three cattle breeds, i.e., Limousin, Holstein-Friesian, and Hereford, using Agilent microarray chips. Pathway and network analysis were performed to select the important biological pathways involved in marbling and their interactions.

Materials and Methods

Transcriptomics Data Set

The study by Sadkowski et al., 2014 compared gene expression in semitendinosus skeletal muscle of well-marbled beef (Holstein-Friesian and Hereford) versus lean-marbled beef (Limousin). Their publicly available microarray data set was used in the present study (NCBI GEO GSE46411). The Holstein-Friesian, Hereford, and Limousin groups consisted of four animals each. Samples for total RNA isolation were taken instantaneously after slaughter from semitendinosus muscle and were kept in liquid nitrogen for transportation and then at −80°C until analyzed. Quality of RNA samples was evaluated using Bioanalyzer 2100 (Agilent Technologies, USA). Only samples with RIN ≥ 8 were further analyzed (Sadkowski et al., 2014).

Agilent Microarray Data Analysis

Global gene expression was measured with Agilent Two-Color Mi Bovine (V2) 4 x 44K Gene Expression Microarray oligonucleotide slides (Agilent, USA). Sadkowski and coworkers checked the quality of the data and performed LOWESS normalization. The normalized transcriptomic data compared well-marbled beef Holstein-Friesian (n = 4) or Hereford (n = 4) to lean-marbled beef (n = 4). The four log10 fold change (log10FC) values for each group comparison were averaged to obtain an estimate of the 10logFC between the entire groups. Furthermore, a one-sample t-test was performed on both sets of four values, comparing those to 0 (giving a p-value indicating the significance of these values being different from 0 = no change). Bovine genes were considered to be significantly, differentially expressed with p ≤ 0.05 and an absolute FC ≥ 1.3 (Sadkowski et al., 2014).

B. taurus Pathway Collection

The online biological pathway repository, WikiPathways (Slenter et al., 2017), contains pathways of different species, however a B. taurus collection was missing. We updated and extended the pathway collection for B. taurus. We also created a B. taurus gene identifier (ID) mapping database based on mappings present in the Ensembl-based BridgeDb framework (van Iersel et al., 2010). The newly created B. taurus ID mapping database was used to annotate genes and proteins in pathways from WikiPathways and to perform pathways analysis. A online and freely available version of the database for the Ensembl build 85 is accessible at (http://bridgedb.org/data/gene_database/archive/r85/Bt_Derby_Ensembl_85.bridge.zip). Second, the WikiPathways homology based the homology mapper which is available at GitHub (https://github.com/PathVisio/homology.mapper) was updated to improve homology coverage for gene products that were annotated with different data sources. The pathways were converted from human pathways, with a required minimum successful conversion of at least 50% of the original human genes. Third, we manually curated all converted pathways to check whether the genes were correctly annotated and pathways are relevant in B. taurus. Finally, new pathways directly derived from cow breeding literature and not present in the WikiPathways collection were designed in PathVisio (v3.2.0) (Kutmon et al., 2015), the pathway creation, visualization, and analysis tool. All pathways were uploaded in gpml format to WikiPathways using the WikiPathways plugin (https://www.pathvisio.org/plugin/wikipathways-plugin/) for PathVisio.

Pathway-Based Over-representation Analysis

To analyze and visualize the molecular changes in marbling at biological process level a pathway-based over-representation analysis was performed in PathVisio (v3.2.0) (Kutmon et al., 2015). The B. taurus WikiPathways pathway collection, containing 286 pathways (6/30/2015), and the B. taurus ID mapping database, was used in the analysis. The pathways are ranked based on a standardized difference score (Z-score) based on the expected value and standard deviation of the number of significantly (p ≤ 0.05) and differentially (absolute FC ≥ 1.3) expressed genes in a pathway. Biological pathways significantly changed when (i) Z-score > 1.96, (ii) permuted p-value < 0.05 and (iii) minimum number of changed genes is 3. Additionally, alterations in gene expression (log10FC and p value) when comparing Hereford to Limousin were visualized on the B. taurus pathways with PathVisio.

Gene Ontology Overrepresentation Analysis

To find the biological processes in which differentially expressed genes were over represented while no pathways for these processes were present in the B. taurus collection at WikiPathways we performed Gene Ontology (GO) analysis via the GO-Elite web-interface (Zambon et al., 2012). GO-Elite is a flexible tool for GO-based over-representation analysis. To identify GO processes the following settings in GO-Elite were used: (i) 2000 permutations, (ii) Z-score GO pruning algorithm, (iii) Z-score threshold >1.96, (iv) p-value threshold <0.05 and (v) minimum number of changed genes is 3 (apart from the method specific permutations those are the same criteria as used for the pathway analysis). This approach not only helps to unify the characteristics and functions of the genes but also to attain a broader perspective of the muscle physiological processes and their influence on meat quality.

Integrated Network Analysis

To visualize the pathway and GO analysis results and their interactions the network analysis and visualization tool, Cytoscape (version 3.2.0), was used (Shannon et al., 2003). First, all enriched pathways and the differentially expressed genes present in these pathways were selected. Second, all changed GO processes and the Differentially Expressed genes present in these GO classes were selected. Third, both results were combined into one network showing the interaction between pathways and GO classes based on corresponding differentially expressed genes. Finally, differences in gene expression between well-marbled and lean-marbled skeletal muscle were visualized in the network.

Results

Identification of Differentially Expressed Genes Between Well-Marbled and Lean-Marbled Skeletal Muscle

In the selected transcriptomic data set of beef marbling 42,990 microarray reporters were measured in both lean marbling beef (Limousin) and well marbling beef (Hereford and Holstein-Friesian) animals. Statistical analysis was performed on 29,677 reported genes that remained from the 42,990 reporters after quality control and annotation with Ensembl gene IDs. In the Hereford breed compared to the Limousin breed, 1,513 were higher expressed and 1,556 lower expressed (absolute log10FC >0.11 and p-value <0.05). When comparing the Holstein-Friesian breed to the Limousin breed, 1,772 genes were higher expressed and 2,458 lower expressed in the Holstein Friesian breed. The genes that met these criteria were used for further analysis.

Creating of B. taurus Pathway Collection and Pathway Design

In total, 282 human pathways were converted from human pathways to cattle pathways. All these pathways were manually checked and are available at (https://www.pathvisio.org/downloads/download-pathways/). Moreover, 4 pathways were newly created based on the bovine breeding literature: Growth hormone signaling (WP2890) (Roudbari Z and Kutmon M: Growth Hormone (GH) Signaling (B. taurus); (https://www.wikipathways.org/instance/WP2890), Growth hormone receptor signaling (WP2891) (Roudbari Z, Hanspers K, Evelo C, Kutmon M: Growth Hormone Receptor (GHR) Signaling (B. taurus) (https://www.wikipathways.org/instance/WP2891), IGF1-signaling (WP2892) (Roudbari Z, Evelo C, Willighagen E, Mélius J, Hanspers K, Kutmon M: IGF1-signaling (B. taurus) (https://www.wikipathways.org/instance/WP2892), and Gonadotropin- releasing hormone signaling (WP2901) (Roudbari Z, Kutmon M, Pico A, Willighagen E, Mélius J: Gonadotropin-releasing hormone (GNRH) signaling (B. taurus) (https://www.wikipathways.org/instance/WP2901). As an example, the newly designed GNRH signaling pathway is shown in Figure 1. The elements of this process are the key factors stimulating gonadotropin release from the pituitary, which controls the release of luteinizing hormone and follicle-stimulating hormone, and reproductive development in mammals.

FIGURE 1

Figure 1 Gonadotropin-releasing hormone (GnRH) signaling (B. taurus) (based on Widmann et al., 2013) and available at https://www.wikipathways.org/instance/WP2901.

Pathway Analysis

When comparing Hereford with Limousin breed ten biological pathways that were formerly known to be involved in marbling (Cui et al., 2012; Lim et al., 2013; Silva-Vignato et al., 2017) were found to be significantly enriched in differentially expressed genes (z-score > 1.96) (Table 1). Four biological processes such as: The Hypertrophy Model, P38 MAPK signaling, IL-1 signaling, and insulin signaling pathways, which are known to be important in marbling development are described in more detail, and the pathways are shown in Results section. Interestingly, some pathways not yet known to play a role in marbling were also found to be enriched in differentially expressed genes when comparing Hereford with Limousin breed. These included histone modifications and vitamin D metabolism pathways, in addition to Hereford breed that is used for meat production, the Holstein-Friesian breed that is a dairy cattle, was also compared with the Limousin breed. Pathway analysis revealed that some but not all of the marbling related pathways found for the comparison between Hereford breed and Limousin breed were also found when comparing Holstein-Friesian breed with Limousin. Examples of consistently affected pathways are the P38 MAPK signaling and the Hypertrophy Model pathways (Table 2).

TABLE 1

Table 1 The highest ranked pathways in skeletal muscle of Hereford compared to Limousin breed.

TABLE 2

Table 2 The highest ranked pathways in skeletal muscle of Holstein-Friesian compared to Limousin breed.

Hypertrophy Model Pathway

Muscle hypertrophy is known to increase the muscle mass, and is determined by increased protein mass per fiber which results from an increase of protein synthesis (Glass, 2005). In the B. taurus Hypertrophy model (http://www.wikipathways.org/instance/WP982) the overall gene expression was higher in skeletal muscle of Hereford and Holstein-Friesian compared to Limousin (Figure 2). The expression of the Il1a, Ifrd1, Cyr61, ATF3, and Ankrd1 genes were significantly higher in Hereford in this pathway and the Il18, Eif4ebp1, and Il1r1 genes were significantly lower in the model (p-value < 0.05). Among them was IL-1 which plays a significant role in lipid metabolism by regulating insulin levels under physiological conditions (Matsuki et al., 2003); Atf3 which works together with p38c in a common pathway in the intestine to regulate lipid metabolism and immune homeostasis (Chakrabarti et al., 2014); and Frd1, Cyr61, and Ankrd1 genes. Two of the seven genes were significantly lower expressed including Ef4ebp1 contributing to the development of obesity through increased adipogenesis and fat metabolism alterations (Le Bacquer et al., 2007) and Il18 gene. TNF and interleukin (IL)-1 may cause negative inotropic effects indirectly through activation or release of IL-18 (Mann, 2015).

FIGURE 2

Figure 2 Skeletal muscle gene expression in Hereford and Holstein-Friesian vs Limousin visualized on the Hypertrophy model. In the hypertrophy model from WikiPathways (WP982) the changes in gene expression between Hereford and Holstein-Friesian with Limousin in skeletal muscle are visualized. The logFC (Hereford and Holstein-Friesian vs Limousin) is indicated with a color gradient (blue to red over white), i.e., blue represents a negative value (= lower expressed in Hereford and Holstein-Friesian) and red a positive value (= higher expressed in Hereford and Holstein-Friesian). The p-value is colored based on a rule, i.e. p-value < = 0.05 (= significant) is shown in green and p-value > 0.05 in white.

p38 MAPK Signaling Pathway

The p38 Mitogen-activated protein kinase (p38 MAPK) signaling pathway has found to be responsible for transduction of extracellular signals to their intracellular targets in different types of cells, including skeletal muscle cells and which leads to several biological effects for example proliferation, differentiation, migration, growth, apoptosis, and more specifically to muscle cells, hypertrophy (Yu et al., 2010; Silva-Vignato et al., 2017). The p38 MAPK is one intracellular signaling pathway activated during the differentiation of myogenic cell lines and this pathway is a chief regulator of skeletal muscle development (Keren et al., 2006). The p38 MAPK signaling pathway is a well-known pathway that affects lipid metabolism (Zhang and Liu, 2002). In the B. taurus p38 MAPK signal pathway five of the seven genes present were significantly higher expressed in Hereford and Holstein-Friesian compared to Limousin (Figure 3).

FIGURE 3

Figure 3 Skeletal muscle gene expression in Hereford and Holstein-Friesian vs Limousin visualized on the p38 MAPK signaling. In the p38 MAPK signaling pathway from WikiPathways (WP1037) the changes in gene expression between Hereford and Holstein-Friesian with Limousin in skeletal muscle are visualized. The logFC (Hereford and Holstein-Friesian vs Limousin) is indicated with a color gradient (blue to red over white), i.e., blue represents a negative value (= lower expressed in Hereford and Holstein-Friesian) and red a positive value (= higher expressed in Hereford and Holstein-Friesian). The p-value is colored based on a rule, i.e., p-value < = 0.05 (= significant) is shown in green and p-value > 0.05 in white.

IL-1 Signaling Pathway

The IL-1 signal pathway is a major mediator of innate immune reactions.This pathway regulates extracellular and intracellular signaling of IL-1α or IL-1β including positive and negative-feedback mechanisms which strengthen or terminate the IL-1 response. In reply to ligand binding of the receptor, a complicated sequence of combinatorial phosphorylation and ubiquitination events lead to activation of nuclear factor kB signaling and the JNK and p38 mitogen-activated protein kinase pathways (Weber et al., 2010). The members of the B. taurus IL-1 signaling pathway (http://www.wikipathways.org/instance/WP3271), such as IL-1α, IL-1β, MAP3K1, UBE2N, MAPK14, REL, ATF2, and JUN were significantly up-regulated in Hereford breed versus Limousin breed (Figure 4). Among them IL-1α which was found to play a role as an inhibitor of the expression of peroxisome proliferator-activated receptor gamma (PPARG), a key transcriptional factor for adipocytes differentiation, (Um et al., 2011). IL-1β has been reported to inhibit adipocyte differentiation from preadipocytes and to reduce the lipid content in mature adipocytes (Simons et al., 2005). Some of aforementioned genes: MAP3K1, UBE2N, MAPK14, REL, ATF2 can directly bind to the peroxisome proliferator-activated receptor promoter and activate transcription to regulate adipocyte differentiation (Maekawa et al., 2010). The significantly down-regulated genes are: RELA, MAPK1, IKBKG, MAP2K4, SQSTM1, and IL1R1 (Figure 4). Some of them are known to participate in lipid metabolism processes; activation of p62/SQSTM1 and peroxisome proliferator-activated receptor gamma is induced by palmitate internalization, which triggers lipid metabolism and limits inflammation (Krausgruber et al., 2011).

FIGURE 4

Figure 4 Skeletal muscle gene expression in Hereford and Holstein-Friesian vs Limousin visualized on IL-1 signaling. In the IL-1 signaling pathway from WikiPathways (WP3271) the changes in gene expression between Hereford and Holstein-Friesian with Limousin in skeletal muscle are visualized. The logFC (Hereford and Holstein-Friesian vs Limousin) is indicated with a color gradient (blue to red over white), i.e., blue represents a negative value (= lower expressed in Hereford and Holstein-Friesian) and red a positive value (= higher expressed in Hereford and Holstein-Friesian). The p-value is colored based on a rule, i.e., p-value < = 0.05 (= significant) is shown in green and p-value > 0.05 in white.

Insulin Signaling Pathway

Genes engaged in the insulin signaling pathway regulate several aspects of cellular function, including most notably the regulation of cellular growth and maintaining glucose homeostasis (DeBosch and Muslin, 2008). For Hereford vs. Limousine comparison, twenty six genes present in the B. taurus Insulin signaling pathway (http://www.wikipathways.org/instance/WP966) showed significant expression differences. Fourteen genes were identified as up-regulated in Hereford breed (SOS2, PIK3R3, PIK3CA, PIK3C2A, CBLB, CBLC, SNAP25, JUN, EGR1, MAP3K1, MAPK14. ENPP1, and XBP1), and twelve were down-regulated in Hereford breed (PFKM, PFKL, ARF1, STXBP2, EIF4EBP1, PIK3CD, GAB1, IGF1R, MAPK1, MAPK13, MAP2K4, and ELK1) when compared to Limousin breed (Figure 5). The involvement of some of upregulated and downregulated genes in lipid accumulation processes were earlier confirmed, including;PIK3CA (Foukas et al., 2013), JUN (Guo et al., 2016), EGR1 (Singh et al., 2015), p38 MAPK lipid accumulation (Sun et al., 2012), NPP1 (Pan et al., 2011), XBP-1 (Zhao et al., 2012), 4E-BP1/2 (Singh et al., 2015), and IGF-1R (Freude et al., 2012).

FIGURE 5

Figure 5 Skeletal muscle gene expression in Hereford and Holstein-Friesian vs Limousin visualized on the Insulin signal pathway. In the Insulin signaling pathway from WikiPathways (WP966) the changes in gene expression between Hereford and Holstein-Friesian with Limousin in skeletal muscle are visualized. The logFC (Hereford and Holstein-Friesian vs Limousin) is indicated with a color gradient (blue to red over white), i.e., blue represents a negative value (= lower expressed in Hereford and Holstein-Friesian) and red a positive value (= higher expressed in Hereford and Holstein-Friesian). The p-value is colored based on a rule, i.e. p-value < = 0.05 (= significant) is shown in green and p-value > 0.05 in white.

Gene Ontology Analysis

Pathway analysis gave an insight in the biological processes involved in marbling. However, only 69% of all measured genes are present in the investigated pathways from the B. taurus WikiPathways collection. In order to obtain a better insight in the biological role of the differentially expressed genes not present in WikiPathways a GO analysis was performed. (Table 3). This approach not only helps to unify the characteristics and functions of the genes but also to attain a broader perspective of the muscle physiological processes and their influence on meat quality related to marbling.

TABLE 3

Table 3 The enriched processes found by GO-Analysis. A description of the process together with the positive gene number, Z Score, and P-value are given.

Integrated Network of Altered Pathways With GO-Terms

The significant pathways and GO terms were merged together and shown in Figure 6. Some of the highly connected nodes are IL1α, TGFB2, PAK1, (as Pak1 deficiency led to upregulation of reverse cholesterol transporters in ApoE−/− mice in response to Western diet feeding, it might be suggested that Pak1 exerts a negative modulatory influence on these transporters and thereby might promote lipid retention in inflamed arteries which cause atherogenesis (Singh et al., 2015), TGFB2 (that TGF-β2 might control adipocyte differentiation in bone marrow stromal cells in vivo by inducing PPARγ phosphorylation. Whether Smad activation induced by TGFβ2 might play a role with MAPK in the inhibition of adipocyte differentiation induced by TGFβ2 in vivo requires more investigation (Ahdjoudj et al., 2005), IL1β, UBE2N, IRAK2, RCAN1, LRP6 (Treatment of LRP6 knockdown-Human mesenchymal stem cells with adipogenic supplements led to the accumulation of fat vacuoles, which was demonstrated by Oil Red O staining (Peröbner et al., 2012), IKBKG, TGFB3, RELA, BCL10, FGF2 (FGF-2 treatment of human preadipocytes also resulted in increased adipocyte differentiation, suggesting that this feature might be common to members of the fibroblast growth factor family, although FGF-1 was consistently the more potent adipogenic agent, particularly in cells from subcutaneous depots (Hutley et al., 2004), IGF1R, and TLR4 (TLR4 knockdown in H9C2 cardio myocytes decreases fatty acid-induced lipid accumulation (Dong et al., 2012) which are present in at least two different pathways and GO terms.

FIGURE 6

Figure 6 Altered pathways and GO terms in skeletal muscle of Hereford vs Limousin. In the network, pathways are shown in orange rectangles and GO Terms are shown in green rectangles. Differentially expressed genes are shown in red (= higher expressed in Hereford) and blue (= lower expressed in Hereford). Linked genes common between two process pathway and GO terms (green diamonds) are depicted with a hexagon shape.

Discussion

The aim of this study was to determine transcriptional profiles of high and low marbled beef with a focus on pathways of muscle cell origin that might play a role in the regulation of marbling development. The regulation of marbling is suggested to be the result of interaction of signaling pathways in muscle, fat, and intramuscular connective tissue (Hocquette, 2010). Identifying these processes with pathway analysis can help to decipher the key processes involved in marbling development. Pathway analysis revealed 17 pathways that were significantly different (z-score > 1.96) between well-marbled and lean marbled breeds. P38 MAPK signaling pathway well known to affects lipid metabolism and muscle development, was enriched when we compared gene expression in well and low marbling breeds. In addition, the signaling pathways “Hypertrophy Model”, “MicroRNAs in cardiomyocyte hypertrophy” and “Physiological and pathological hypertrophy of the heart” that play a role in tissue development were affected. Interestingly, the analyses also demonstrated that pathways related to immune response (IL signaling, TCR signaling, and Toll-like receptor signaling pathways) and insulin signaling, mitochondrial gene expression and vitamin D metabolism were enriched and might act together with pathways related to lipid metabolism. We explored regulatory pathways that control gene expression in bovine muscle and the relationships between gene expression and the marbling trait to identify markers that effect on marbling. A similar study done by (Hong et al., 2014) investigated the biological characteristics of differentially expressed genes in high marbled muscle in pig compared to a low marbled muscle. They indicated that the differentially expressed genes were clustered to three group related to energy metabolism, protein synthesis, and immune response in high marbling pigs. These finding suggested that the genes related to energy metabolism, protein synthesis, and immune response contribute to growth performance and meat quality. Our results also showed differentially expressed genes take part in these processes. The hypertrophy model pathway was found to be enriched with the highest Z score in the present study and during muscle hypertrophy there is an equilibrium between protein synthesis and degradation that might bring about protein deposition, and hence muscle growth. Together, these processes will lead to differences in muscle and fat deposition, and for this reason animals have different proportions of ribeye area and back fat thickness (Silva-Vignato et al.,2017). Mitogen-activated protein kinase (MAPK) signals have been shown to play a significant role in intracellular signaling associated with a variety of cellular activities including cell proliferation, differentiation, survival, and death (Yu et al., 2010). In mammalian cells, three MAPK families have been characterized: classical MAPK (also known as ERK), C-Jun N-terminal kinase/stress activated protein kinase and p38 kinase pathways (Zhang and Liu, 2002). Each mammalian MAPK signaling route comprises at least three components: a MAPK kinase kinase (MAP3K), a MAPK kinase (MAP2K), and a MAPK. Activated MAPKs phosphorylate various substrate proteins including transcription factors such as ATF2 and Jun (Kim & Choi, 2010). Philip and coworkers (Philip et al., 2005) discovered that the p38 MAPK played a key role in GDF-8-induced inhibition of proliferation and upregulation of the cyclin kinase inhibitor p21. In addition, their results showed a functional link between the p38 MAPK and GDF-8-activated Smad pathways, and identify an important role for the p38 MAPK in GDF-8’s function as a negative regulator of muscle growth (Philip et al., 2005). In comparative muscle transcriptome associated with carcass traits of Nellore cattle, Silva-Vignato and colleagues indicated that MAPK signaling pathway involved in muscle and fat deposition, which are economically important carcass traits for beef production (Silva-Vignato et al., 2017). The third pathway found in the present study was IL-1, the IL-1 family of cytokines includes 11 proteins encoded by 11 different genes and gene regulation of IL-1 signal is activation of MKK4, MKK3, and MKK6 gene that activate NF-κB and p38 MAPK pathways (Weber et al., 2010). These two signaling pathway are needed to upregulate the expression of the key E3 ligases, MuRF1, which mediate the inhibition of protein synthesis (Clarke et al., 2007). Moreover, the insulin-signal transduction pathway, which was another pathway identified in the present study, is a highly conserved pathway that regulates cellular growth and when insulin binding to its cell-surface receptor, insulin receptor, activates a complex intracellular signaling network through insulin substrate proteins and the canonical PI3K and ERK cascades (Hocquette et al., 2010). Interestingly, insulin signaling is one of important factors involved in muscle development since stimulation of glucose utilization in fat and muscle cells in calves is occurring by enhancing insulin intracellular signaling (Jovanović et al., 2017).

Currently, systems biology approaches have become one of the most effective manners to accelerate the genetic improvement of beef and dairy cattle herds (Kadarmideen, 2014). It allows the selection of desired characteristics through the use of transcriptome profiles. The development of high throughput data and bioinformatics tools allow the selection of superior breeds without wasting time and money, contributing to the widespread use of transcriptome analysis in beef cattle operations. Our study shows in cattle that integration of pathway expression profiles in a systems biology approach will contribute to a better understanding of the genes and regulatory processes involved in marbling. These novel insights can be used in the future to take into account when improving the meat quality in beef cattle. The molecular mechanisms which underlie fat content in muscle can provide vital information for the production of healthier beef for human consumption.

Conclusions

The outcome of our research is the identification of biological pathways where we highlighted changed genes which are related with marbling in beef cattle. These results give a better understanding of mechanisms involving marbling in beef cattle, which is economically important carcasses trait for meat quality. Moreover, the genes involved in the highlighted pathways can potentially be utilized as an early biological marker for marbling fat content in breed-specific differences in growth performance and meat quality.

Data Availability Statement

The datasets analyzed for this study can be found in the NCBI Gene Expression Omnibus (GEO) under the accession number GSE46411 and in the wikipathways.org knowledge base under the accession numbers: WikiPathways: WP2890, WikiPathways: WP2891, WikiPathways: WP2891, and WikiPathways: WP2901.

Ethics Statement

Ethical review and approval was not required for the animal study because this study is just analysis.

Author Contributions

ZR, SC, and CE conceptualized and designed the study. ZR, SC, MK, LE, and JM analyzed the data. ZR collected and assembled the data and wrote the manuscript and SC, CE, and TS contributed to reviewing and editing the manuscript and also read and approved the final manuscript.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

The publication of the manuscript was financed by the KNOW (Leading National Research Centre) Scientific Consortium “Healthy Animal–Safe Food,” under Ministry of Science and Higher Education Decision No. 05-1/KNOW2/2015 (Poland).

References

Ahdjoudj, S., Kaabeche, K., Holy, X., Fromigue, O., Modrowski, D., Zerath, E., et al. (2005). Transforming growth factor-β inhibits CCAAT/enhancer-binding protein expression and PPARγ activity in unloaded bone marrow stromal cells. Exp. Cell Res. 303 (1), 138–147. doi: 10.1016/j.yexcr.2004.09.013

PubMed Abstract | CrossRef Full Text | Google Scholar

Bernard, C., Cassar-Malek, I., Le Cunff, M., Dubroeucq, H., Renand, G., Hocquette, J. F. (2007). New indicators of beef sensory quality revealed by expression of specific genes. J. Agric. Food Chem. 55 (13), 5229–5237. doi: 10.1021/jf063372l

PubMed Abstract | CrossRef Full Text | Google Scholar

Cesar, A. S. M., Regitano, L. C. A., Koltes, J. E., Fritz-Waters, E. R., Lanna, D. P. D., Gasparin, G., et al. (2015). Putative regulatory factors associated with intramuscular fat content. PloS One 10 (6), e0128350. doi: 10.1371/journal.pone.0128350

PubMed Abstract | CrossRef Full Text | Google Scholar

Chakrabarti, S., Poidevin, M., Lemaitre, B. (2014). The Drosophila MAPK p38c regulates oxidative stress and lipid homeostasis in the intestine. PloS Genet. 10 (9), e1004659. doi: 10.1371/journal.pgen.1004659

PubMed Abstract | CrossRef Full Text | Google Scholar

Clarke, B. A., Drujan, D., Willis, M. S., Murphy, L. O., Corpina, R. A., Burova, E., et al. (2007). The E3 Ligase MuRF1 degrades myosin heavy chain protein in dexamethasone treated skeletal muscle. Cell Metab. 6 (5), 376–385. doi: 10.1016/j.cmet.2007.09.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Cui, H.-X., Liu, R.-R., Zhao, G.-P., Zheng, M.-Q., Chen, J. L., Wen, J. (2012). Identification of differentially expressed genes and pathways for intramuscular fat deposition in pectoralis major tissues of fast-and slow-growing chickens. BMC Genomics 13 (1), 213. doi: 10.1186/1471-2164-13-213

PubMed Abstract | CrossRef Full Text | Google Scholar

DeBosch, B. J., Muslin, A. J. (2008). Insulin signaling pathways and cardiac growth. J. Mol. Cell. Cardiol. 44 (5), 855–864. doi: 10.1016/j.yjmcc.2008.03.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Dong, B., Qi, D., Yang, L., Huang, Y., Xiao, X., Tai, N., et al. (2012). TLR4 regulates cardiac lipid accumulation and diabetic heart disease in the nonobese diabetic mouse model of type 1 diabetes. Am. J. Physiol.-Heart Circulatory Physiol. 303 (6), H732–H742. doi: 10.1152/ajpheart.00948.2011

CrossRef Full Text | Google Scholar

Foukas, L. C., Bilanges, B., Bettedi, L., Pearce, W., Ali, K., Sancho, S., et al. (2013). Long-term p110α PI3K inactivation exerts a beneficial effect on metabolism. EMBO Mol. Med. 5 (4), 563–571. doi: 10.1002/emmm.201201953

PubMed Abstract | CrossRef Full Text | Google Scholar

Freude, S., Schilbach, K., Hettich, M. M., Brönneke, H. S., Zemva, J., Krone, W., et al. (2012). Neuron-specific deletion of a single copy of the insulin-like growth factor-1 receptor gene reduces fat accumulation during aging. Hormone Metab. Res. 44 (2), 99–104. doi: 10.1055/s-0031-1298018

CrossRef Full Text | Google Scholar

Glass, D. J. (2005). Skeletal muscle hypertrophy and atrophy signaling pathways. Int. J. Biochem. Cell Biol. 37 (10), 1974–1984. doi: 10.1016/j.biocel.2005.04.018

PubMed Abstract | CrossRef Full Text | Google Scholar

Guo, J., Fang, W., Sun, L., Lu, Y., Dou, L., Huang, X., et al. (2016). Reduced miR-200b and miR-200c expression contributes to abnormal hepatic lipid accumulation by stimulating JUN expression and activating the transcription of srebp1. Oncotarget 7 (24), 36207. doi: 10.18632/oncotarget.9183

PubMed Abstract | CrossRef Full Text | Google Scholar

Hocquette, J. F., Gondret, F., Baeza, E., Medale, F., Jurie, C., Pethick, D. W. (2010). Intramuscular fat content in meat-producing animals:development, genetic and nutritional control, and identification of putative markers. Animal 4 (2), 303–319. doi: 10.1017/S1751731109991091

PubMed Abstract | CrossRef Full Text | Google Scholar

Hocquette, J. F., Cassar-Malek, I., Jurie, C., Bauchart, D., Picard, B., Renand, G. (2012). Relationships between muscle growth potential, intramuscular fat content and different indicators of muscle fibre types in young Charolais bulls. Anim. Sci. J. 83 (11), 750–758. doi: 10.1111/j.1740-0929.2012.01021.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Hocquette, J. F. (2010). Endocrine and metabolic regulation of muscle growth and body composition in cattle. Animal. Int. J. Anim. Biosci. 4 (11), 1797–1809. doi: 10.1017/S1751731110001448

CrossRef Full Text | Google Scholar

Hong, X. U., Huang, Y., Li, W., Yang, M., Ge, C., Zhang, X., et al. (2014). Muscle biological characteristics of differentially expressed genes in wujin and landrace pigs. J. Integr. Agric. 13 (10), 2236–2242. doi: 10.1016/S2095-3119(13)60605-X

CrossRef Full Text | Google Scholar

Hutley, L., Shurety, W., Newell, F., McGeary, R., Pelton, N., Grant, J., et al. (2004). Fibroblast growth factor 1: a key regulator of human adipogenesis. Diabetes 53 (12), 3097–3106. doi: 10.2337/diabetes.53.12.3097

PubMed Abstract | CrossRef Full Text | Google Scholar

Jovanović, L., Pantelić, M., Prodanović, R., Vujanac, I., Đurić, M., Tepavčević, S., et al. (2017). Effect of peroral administration of chromium on insulin signaling pathway in skeletal muscle tissue of Holstein calves. Biol. Trace Element Res. 180 (2), 223–232. doi: 10.1007/s12011-017-1007-1

CrossRef Full Text | Google Scholar

Kadarmideen, H. N. (2014). Genomics to systems biology in animal and veterinary sciences: Progress, lessons and opportunities. Livestock Sci. 166, 232–248. doi: 10.1016/j.livsci.2014.04.028

CrossRef Full Text | Google Scholar

Keren, A., Tamir, Y., Bengal, E. (2006). The p38 MAPK signaling pathway: a major regulator of skeletal muscle development. Mol. Cell. Endocrinol. 252 (1–2), 224–230. doi: 10.1016/j.mce.2006.03.017

PubMed Abstract | CrossRef Full Text | Google Scholar

Kim, E. K., Choi, E.-J. (2010). Pathological roles of MAPK signaling pathways in human diseases. Biochim. Biophys. Acta (BBA)-Mol. Basis Dis. 1802 (4), 396–405. doi: 10.1016/j.bbadis.2009.12.009

CrossRef Full Text | Google Scholar

Krausgruber, T., Blazek, K., Smallie, T., Alzabin, S., Lockstone, H., Sahgal, N., et al. (2011). IRF5 promotes inflammatory macrophage polarization and T H 1-T H 17 responses. Nat. Immunol. 12 (3), 231. doi: 10.1038/ni.1990

PubMed Abstract | CrossRef Full Text | Google Scholar

Kutmon, M., van Iersel, M. P., Bohler, A., Kelder, T., Nunes, N., Pico, A. R., et al. (2015). PathVisio 3: an extendable pathway analysis toolbox. PloS Comput. Biol. 11 (2), e1004085. doi: 10.1371/journal.pcbi.1004085

PubMed Abstract | CrossRef Full Text | Google Scholar

Le Bacquer, O., Petroulakis, E., Paglialunga, S., Poulin, F., Richard, D., Cianflone, K., et al. (2007). Elevated sensitivity to diet-induced obesity and insulin resistance in mice lacking 4E-BP1 and 4E-BP2. J. Clin. Invest. 117 (2), 387–396. doi: 10.1172/JCI29528

PubMed Abstract | CrossRef Full Text | Google Scholar

Lee, S.-H., Gondro, C., van der Werf, J., Kim, N.-K., Lim, D., Park, E.-W., et al. (2010). Use of a bovine genome array to identify new biological pathways for beef marbling in Hanwoo (Korean Cattle). BMC Genomics 11 (1), 1–11. doi: 10.1186/1471-2164-11-623

PubMed Abstract | CrossRef Full Text | Google Scholar

Lim, D., Lee, S. H., Kim, N. K., Cho, Y. M., Chai, H. H., Seong, H. H., et al. (2013). Gene co-expression analysis to characterize genes related to marbling trait in Hanwoo (Korean) cattle. Asian-Australasian J. Anim. Sci. 26 (1), 19–29. doi: 10.5713/ajas.2012.12375

CrossRef Full Text | Google Scholar

Maekawa, T., Jin, W., Ishii, S.. (2010). The role of ATF-2 family transcription factors in adipocyte differentiation: antiobesity effects of p38 inhibitors. Mol. Cel. Biol. 30 (3), 613C625. doi: 10.1128/MCB.00685-09

CrossRef Full Text | Google Scholar

Mann, D. L. (2015). Innate immunity and the failing heart: the cytokine hypothesis revisited. Circ. Res. 116 (7), 1254–1268. doi: 10.1161/CIRCRESAHA.116.302317

PubMed Abstract | CrossRef Full Text | Google Scholar

Matsuki, T., Horai, R., Sudo, K., Iwakura, Y. (2003). IL-1 plays an important role in lipid metabolism by regulating insulin levels under physiological conditions. J. Exp. Med. 198 (6), 877–888. doi: 10.1084/jem.20030299

PubMed Abstract | CrossRef Full Text | Google Scholar

McAfee, A. J., McSorley, E. M., Cuskelly, G. J., Moss, B. W., Wallace, J. M. W., Bonham, M. P., et al. (2010). Red meat consumption: an overview of the risks and benefits. Meat Sci. 84 (1), 1–13. doi: 10.1016/j.meatsci.2009.08.029

PubMed Abstract | CrossRef Full Text | Google Scholar

O’Connor, S. F., Tatum, J. D., Wulf, D. M., Green, R. D., Smith, G. C. (1997). Genetic effects on beef tenderness in Bos indicus composite and Bos taurus cattle. J. Anim. Sci. 75 (7), 1822–1830. doi: 10.2527/1997.7571822x

PubMed Abstract | CrossRef Full Text | Google Scholar

Pan, W., Ciociola, E., Saraf, M., Tumurbaatar, B., Tuvdendorj, D., Prasad, S., et al. (2011). Metabolic consequences of ENPP1 overexpression in adipose tissue. Am. J. Physiol.-Endocrinol. Metab. 301 (5), E901–E911. doi: 10.1152/ajpendo.00087.2011

PubMed Abstract | CrossRef Full Text | Google Scholar

Peröbner, I., Karow, M., Jochum, M., Neth, P. (2012). LRP6 mediates Wnt/β-catenin signaling and regulates adipogenic differentiation in human mesenchymal stem cells. Int. J. Biochem. Cell Biol. 44 (11), 1970–1982. doi: 10.1016/j.biocel.2012.07.025

PubMed Abstract | CrossRef Full Text | Google Scholar

Philip, B., Lu, Z., Gao, Y. (2005). Regulation of GDF-8 signaling by the p38 MAPK. Cell. Signalling 17 (3), 365–375. doi: 10.1016/j.cellsig.2004.08.003

CrossRef Full Text | Google Scholar

Sadkowski, T., Ciecierska, A., Majewska, A., Oprządek, J., Dasiewicz, K., Ollik, M., et al. (2014). Transcriptional background of beef marbling—novel genes implicated in intramuscular fat deposition. Meat Sci. 97 (1), 32–41. doi: 10.1016/j.meatsci.2013.12.017

PubMed Abstract | CrossRef Full Text | Google Scholar

Shackelford, S. D., Koohmaraie, M., Miller, M. F., Crouse, J. D., Reagan, J. O. (1991). An evaluation of tenderness of the longissimus muscle of angus by Hereford versus brahman crossbred heifers. J. Anim. Sci. 69 (1), 171–177. doi: 10.2527/1991.691171x

PubMed Abstract | CrossRef Full Text | Google Scholar

Shannon, P., Markiel, A., Ozier, O., Baliga, N. S., Wang, J. T., Ramage, D., et al. (2003). Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 13 (11), 2498–2504. doi: 10.1101/gr.1239303.metabolite

PubMed Abstract | CrossRef Full Text | Google Scholar

Silva-Vignato, B., Coutinho, L. L., Cesar, A. S. M., Poleti, M. D., Regitano, L. C. A., Balieiro, J. C. C. (2017). Comparative muscle transcriptome associated with carcass traits of nellore cattle. BMC Genomics 18 (1), 506. doi: 10.1186/s12864-017-3897-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Simons, P. J., van den Pangaart, P. S., van Roomen, C. P. A. A., Aerts, J. M. F. G., Boon, L. (2005). Cytokine-mediated modulation of leptin and adiponectin secretion during in vitro adipogenesis: evidence that tumor necrosis factor-α-and interleukin-1β-treated human preadipocytes are potent leptin producers. Cytokine 32 (2), 94–103. doi: 10.1016/j.cyto.2005.08.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Singh, M., Shin, Y.-K., Yang, X., Zehr, B., Chakrabarti, P., Kandror, K. V. (2015). 4E-BPs control fat storage by regulating the expression of Egr1 and ATGL. J. Biol. Chem. 290 (28), 17331–17338. doi: 10.1074/jbc.M114.631895

PubMed Abstract | CrossRef Full Text | Google Scholar

Singh, N. K., Kotla, S., Dyukova, E., Traylor, J. G., Jr., Orr, A. W., Chernoff, J., et al. (2015). Disruption of p21-activated kinase 1 gene diminishes atherosclerosis in apolipoprotein E-deficient mice. Nat. Commun. 6, 7450. doi: 10.1038/ncomms8450

PubMed Abstract | CrossRef Full Text | Google Scholar

Slenter, D. N., Kutmon, M., Hanspers, K., Riutta, A., Windsor, J., Nunes, N., et al. (2017). WikiPathways: a multifaceted pathway database bridging metabolomics to other omics research. Nucleic Acids Res. 46 (D1), D661–D667. doi: 10.1093/nar/gkx1064

CrossRef Full Text | Google Scholar

Sun, C., Qi, R., Wang, L., Yan, J., Wang, Y. (2012). p38 MAPK regulates calcium signal-mediated lipid accumulation through changing VDR expression in primary preadipocytes of mice. Mol. Biol. Rep. 39 (3), 3179–3184. doi: 10.1007/s11033-011-1084-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Takasuga, A., Watanabe, T., Mizoguchi, Y., Hirano, T., Ihara, N., Takano, A., et al. (2007). Identification of bovine QTL for growth and carcass traits in japanese black cattle by replication and identical-by-descent mapping. Mamm. Genome 18 (2), 125–136. doi: 10.1007/s00335-006-0096-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Thaller, G., Kühn, C., Winter, A., Ewald, G., Bellmann, O., Wegner, J., et al. (2003). DGAT1, a new positional and functional candidate gene for intramuscular fat deposition in cattle. Anim. Genet. 34 (5), 354–357. doi: 10.1046/j.1365-2052.2003.01011.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Um, J.-Y., Rim, H.-K., Kim, S.-J., Kim, H.-L., Hong, S.-H. (2011). Functional polymorphism of IL-1 alpha and its potential role in obesity in humans and mice. PloS One 6 (12), e29524. doi: 10.1371/journal.pone.0029524

PubMed Abstract | CrossRef Full Text | Google Scholar

van Iersel, M. P., Pico, A. R., Kelder, T., Gao, J., Ho, I., Hanspers, K., et al. (2010). The BridgeDb framework: standardized access to gene, protein and metabolite identifier mapping services. BMC Bioinf. 11 (1), 5. doi: 10.1186/1471-2105-11-5

CrossRef Full Text | Google Scholar

Weber, A., Wasiliew, P., Kracht, M. (2010). Interleukin-1 (IL-1) pathway. Sci. Signal. 3 (105), cm1–cm1. doi: 10.1126/scisignal.3105cm1

PubMed Abstract | CrossRef Full Text | Google Scholar

Widmann, P., Reverter, A., Fortes, M. R. S., Weikard, R., Suhre, K., Hammon, H., et al. (2013). A systems biology approach using metabolomic data reveals genes and pathways interacting to modulate divergent growth in cattle. BMC Genomics, 14 (1), 798. doi: 10.1186/1471-2164-14-798

PubMed Abstract | CrossRef Full Text | Google Scholar

Wulf, D. M., Tatum, J. D., Green, R. D., Morgan, J. B., Golden, B. L., Smith, G. C. (1996). Genetic influences on beef longissimus palatability in charolais- and limousin-sired steers and heifers. J. Anim. Sci. 74 (10), 2394–2405. doi: 10.2527/1996.74102394x

PubMed Abstract | CrossRef Full Text | Google Scholar

Yamada, T., Taniguchi, Y., Nishimura, S., Yoshioka, S., Takasuga, A., Sugimoto, Y., et al. (2006). Radiation hybrid mapping of genes showing intramuscular fat deposition-associated expression changes in bovine musculus longissimus muscle. Anim. Genet. 37 (2), 184–185. doi: 10.1111/j.1365-2052.2006.01426.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Yu, W., Chen, C., Fu, Y., Wang, X., Wang, W. (2010). Insulin signaling: a possible pathogenesis of cardiac hypertrophy. Cardiovasc. Ther. 28 (2), 101–105. doi: 10.1111/j.1755-5922.2009.00120.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Zambon, A. C., Gaj, S., Ho, I., Hanspers, K., Vranizan, K., Evelo, C. T., et al. (2012). GO-Elite: a flexible solution for pathway and ontology over-representation. Bioinformatics 28 (16), 2209–2210. doi: 10.1093/bioinformatics/bts366

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, W., Liu, H. T. (2002). MAPK signal pathways in the regulation of cell proliferation in mammalian cells. Cell Res. 12 (1), 9. doi: 10.1038/sj.cr.7290105

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhao, S., Zhu, L., Duan, H., Liu, S., Liu, Q., Liu, W., et al. (2012). PI3K/Akt pathway mediates high glucose-induced lipid accumulation in human renal proximal tubular cells via spliced XBP-1. J. Cell. Biochem. 113 (10), 3288–3298. doi: 10.1002/jcb.24207

PubMed Abstract | CrossRef Full Text | Google Scholar