Gene Expression Patterns Associated with Peroxisome Proliferator-activated Receptor (PPAR) Signaling in the Longissimus dorsi of Hanwoo (Korean Cattle).

Adipose tissue deposited within muscle fibers, known as intramuscular fat (IMF or marbling), is a major determinant of meat quality and thereby affects its economic value. The biological mechanisms that determine IMF content are therefore of interest. In this study, 48 genes involved in the bovine peroxisome proliferator-activated receptor signaling pathway, which is involved in lipid metabolism, were investigated to identify candidate genes associated with IMF in the longissimus dorsi of Hanwoo (Korean cattle). Ten genes, retinoid X receptor alpha, peroxisome proliferator-activated receptor gamma (PPARG), phospholipid transfer protein, stearoyl-CoA desaturase, nuclear receptor subfamily 1 group H member 3, fatty acid binding protein 3 (FABP3), carnitine palmitoyltransferase II, acyl-Coenzyme A dehydrogenase long chain (ACADL), acyl-Coenzyme A oxidase 2 branched chain, and fatty acid binding protein 4, showed significant effects with regard to IMF and were differentially expressed between the low- and high-marbled groups (p<0.05). Analysis of the gene co-expression network based on Pearson's correlation coefficients identified 10 up-regulated genes in the high-marbled group that formed a major cluster. Among these genes, the PPARG-FABP4 gene pair exhibited the strongest correlation in the network. Glycerol kinase was found to play a role in mediating activation of the differentially expressed genes. We categorized the 10 significantly differentially expressed genes into the corresponding downstream pathways and investigated the direct interactive relationships among these genes. We suggest that fatty acid oxidation is the major downstream pathway affecting IMF content. The PPARG/RXRA complex triggers activation of target genes involved in fatty acid oxidation resulting in increased triglyceride formation by ATP production. Our findings highlight candidate genes associated with the IMF content of the loin muscle of Korean cattle and provide insight into the biological mechanisms that determine adipose deposition within muscle.

INTRODUCTION

Intramuscular fat (IMF or marbling) is deposited adipose tissue within muscle fibers in skeletal muscle. Marbling plays an important role in determining the quality of meat and therefore its economic value. Identifying the genes that regulate marbling may therefore help to optimize meat production in the livestock industry. Over the last few decades, studies have attempted to identify the genes that cause phenotypic differences in marbling using quantitative trait loci mapping, genome-wide association studies, microarrays, and proteome analyses. These results have demonstrated that marbling is controlled by multiple genes; i.e., it is a quantitative trait. Therefore, an analysis of multiple genes, especially genes related to adipogenesis or lipid metabolism, may lead to a better understanding of the effects of these processes on the phenotypes of marbled skeletal muscle. Recently, He et al. (2013) reported that 13 tag single nucleotide polymorphisms in the peroxisome proliferator-activated receptor (PPAR) signaling pathway are associated with porcine meat quality traits. On investigating fat deposition and muscle growth, seven single nucleotide polymorphisms in the reverse cholesterol transport pathway were found to be significantly associated with meat quality phenotypes in a Wagyu×Limousin reference population (Daniels et al., 2010). We recently revealed that three mitogen-activated protein kinase pathway-related genes participated in regulating the IMF content of Hanwoo, Korean native cattle (Kim et al., 2011). The skeletal muscle genes within the glycolytic pathway, which include protein kinase adenosine monophosphate-activated gamma 3 non-catalytic subunit (PRKAG3), phosphoglycerate mutase 2 (PGAM2), and pyruvate kinase muscle 2 (PKM2), have been described as markers of back fat thickness and average daily gain in Italian large white pigs (Fontanesi et al., 2008). These genes are therefore of interest with regard to understanding the factors that affect the meat quality of livestock. The PPAR signaling pathway, which is involved in lipid metabolism, is known to be an important biological pathway in determining meat quality in mammals. Peroxisome proliferator-activated receptor gamma, the primary transcription regulator in PPAR signaling, is reported to be a key factor in controlling the transcription of many genes that are involved in biological pathways based on adipogenesis (Fernyhough et al., 2007). It remains to be determined which genes of the PPAR signaling pathway play a role in livestock meat quality and more specifically, phenotypic differences in marbling.

The objective of this study was to identify genes that are differentially expressed in low- and high-marbled muscle of Korean cattle (Hanwoo) and to investigate their potential as biomarkers pertaining to the IMF content of beef.

MATERIALS AND METHODS

Animals and sample preparation

Twenty samples of low- and high-marbled m. longissimus from Hanwoo (n = 10 animals in each group) steers were obtained from the junction between the 11th and 12th lumbar vertebrae within 30 min of slaughter. These animals were bred under the same feeding conditions at the Hanwoo Experiment Station of the National Institute of Animal Science (NIAS) in Korea. The selected tissues were placed in liquid nitrogen, ground to a fine powder using a mortar, and stored at −80°C. All experimental procedures and the care of the animals were conducted in accordance with the guidelines of the Animal Care and Use Committee of the NIAS in Korea. The fat content of the m. longissimus muscle samples was analyzed using the methods established by the Association of Official Analytical Chemists (Table 1).

Table 1

Summary statistics of tissue samples for gene expression analysis

Group	Age (months)	Carcass weight (kg)	Intramuscular fat content (%)
Low	27.1 ±0.276	409.2 ±10.588	7.44 ±0.463
High	29.5 ±0.477	478.17 ±13.618**	23.54 ±0.881**

Significant differences (p<0.05 and 0.01, respectively) between low- and high-marbled groups determined using an analysis of variance model.

Quantitative real-time polymerase chain reaction and statistical analysis

The mRNA levels were analyzed by quantitative real-time polymerase chain reaction (qRT-PCR) using gene-specific primer sets (Table 2). Total RNA was prepared from each tissue sample (100 mg) using TRIzol reagent (Invitrogen Life Technologies, Carlsbad, CA, USA). For qRT-PCR analysis, 2 μg RNA was reverse transcribed in a 20 μL reaction volume using random primers (Promega, Madison, WI, USA) and reverse transcriptase (SuperScript II Reverse Transcriptase, Invitrogen Life Technologies). The reactions were incubated at 65°C for 5 min, 42°C for 50 min, and then at 70°C for 15 min to inactivate the reverse transcriptase.

Table 2

Quantitative real-time polymerase chain reaction results for the 48 genes identified by analysis of variance

Gene	Full name	Primer sequence		Expression1		p value2
		Forward	Reverse	Low	High
PPARA	Peroxisome proliferator-activated receptor alpha	ATTCCGAAACCGCGGATCTCAAGT	AAAGGCGGGTTGTTGTTGGTCTTC	0.0997	0.1261	0.308
RXRA	Retinoid X receptor, alpha	TTCACCAAGCACATCTGTGCCATC	GCTTGTCAATCAGGCAGTCCTTGT	0.0883	0.1685	0.042*
RXRB	Retinoid X receptor, beta	TGGCAAACGGCTATGTGCAATCTG	CTTACGGATGGTGCGCTTGAAGAA	0.0376	0.0325	0.784
RXRG	Retinoid X receptor, gamma	AAGAGGGAAGCTGTGCAGGAAGAA	TATCTTCATGGCCACTGTTGGCAC	0.0150	0.0251	0.153
APOA1	Apolipoprotein A-I	AGCTCAACCTGAAACTCCTGGACA	TTTCCAGGTTGTCCCAGAACTCCT	0.0042	0.0054	0.248
APOA2	Apolipoprotein A-II	GAAGCTGCTTGCATTGACTGTGCT	AAAGTAGGCCTTGGCTTGGGTTTG	0.0003	0.0004	0.158
PPARD	Peroxisome proliferator-activated receptor delta	ATGTCGCACAATGCCATCCGTTTC	TGAAGTTCTTCAGGTAGGCGCTGT	0.0762	0.0631	0.465
PPARG	Peroxisome proliferator-activated receptor gamma	CGATATCGACCAACTGAACC	ACCGGTGATTTGTCTGTCGT	0.0098	0.0511	0.021*
APOA5	Apolipoprotein A-V	AGCGCATCCAGCAGAACCTAGA	TTGAGCTGCTTCCTCCGCAT	0.0001	0.0001	0.647
APOC3	Apolipoprotein C-III	TGGACAAGATGCAGGGCTATGTGA	AGCTGAAGCTCTCAGTCATCCAGT	0.0004	0.0007	0.113
PLTP	Phospholipid transfer protein	ACGTGGCCTTCTCTGAGTTCTTCT	TCATCAGGACGATGCTCCCAAAGT	0.0205	0.0658	0.001**
SCD	Stearoyl-coa desaturase	CCAGAGGAGGTACTACAAACCTG	AGCCAGGTGACGTTGAGC	0.8597	3.0974	0.002**
CYP7A1	Similar to cholesterol 7alpha-hydroxylase	GCAATGAAAGCGGCCACTGAAGAA	TACTGTCTAGCACGGGCATGTTGT	0.0003	0.0004	0.439
CYP8B1	Cytochrome P450, family 8, subfamily B, polypeptide	ATGCGCAGAGGAAGCTAGACTTTG	AAGAGATGCTGGGATTGAGCTCCT	0.0013	0.0015	0.640
FADS2	Fatty acid desaturase 2	ACGATTACGGCCACCTCTCTGTTT	AGGCACCCTTTAAGTGACCGATGA	0.0177	0.0126	0.450
LR1H3	Nuclear receptor subfamily 1, group H, member	CTGAGTTTGCCTTGCTCATTGCCA	ATGTGTGTTGCAGCCTCTCTACCT	0.0098	0.0183	0.034
FABP1	Fatty acid binding protein 1	GGAAATCGTGCAGAATGGGAAGCA	TTCCTGCTGAACCACTGCCTTGAT	0.0025	0.0024	0.740
SLC27A4	Solute carrier family 27 (fatty acid transporter) member 1	AGCCTGGTCAAGTTCTGTTCTGGA	AGAAGAGTCGATCATCCATGCCCT	0.0417	0.0473	0.304
FABP3	Fatty acid binding protein 3	ACATGAAGTCACTCGGTGTCGGTT	TTCCTGTCATCTGCTGTGGTCTCA	1.7208	1.0294	0.023*
SLC27A1	Solute carrier family 27 (fatty acid transporter), member4	CAGGGCACCAACGACAAGAAGATT	AGCTCGTCCATCACTAGCACATCA	0.0002	0.0003	0.471
CD36	CD36 molecule (thrombospondin receptor)	TCCAGACAACCACTGTTTCTGCAC	GCATTTGCCAATGTCTAGCACACC	0.2479	0.2026	0.544
OLR1	Oxidized low density lipoprotein (lectin-like) receptor 1	ACTTCACCGCCAGAACCTGAATCT	GTGGGCATCCAAAGACAAGCAGTT	0.0002	0.0003	0.268
LPL	Lipoprotein lipase	TACCCTGCCTGAAGTTTCCAC	CCCAGTTTCAGCCAGACTTTC	0.0820	0.1676	0.226
ACSL6	Acyl-coa synthetase long-chain family member 6	ATTGTGTGTTTCACGAGTGGCACG	CAAACATGTGAGCCAAAGGCAGGA	0.0001	0.0002	0.317
ACSL5	Acyl-coa synthetase long-chain family member 5	TGACAGTCTGTGGGACAAGCTCAT	AGCCAGCAGTGCATTCTGTTTGAC	0.0008	0.0006	0.320
ACSL4	Acyl-coa synthetase long-chain family member 4	GCCTTTGGCTCATGTGCTTGAGTT	TTTGCTGGACTGGTCAGAGAGTGT	0.0024	0.0023	0.663
ACSL1	Acyl-coa synthetase long-chain family member 1	CCCAAAGGAGCAATGATCACGCAT	AGCTCCATGACAGAGCATCACACA	0.0230	0.0178	0.436
CBT1B	Carnitine palmitoyltransferase 1B	AGGCTGCTGAGAAGCACCAGAATA	ACTTGGAGACCACGTAAAGGCAGA	0.0055	0.0050	0.932
CPT1C	Carnitine palmitoyltransferase 1C	TGCACTGGGACATTCCAGACAAGA	TGTCTGAAGAGAGATGGCAGCGTT	0.0002	0.0003	0.467
CPT2	Carnitine palmitoyltransferase II	AACTCTTCACCGATGACAAGGCCA	ATGAGCCTGGATTTCAGAGGCACT	0.0020	0.0076	0.001*
ACAA1	Acetyl-Coenzyme A acyltransferase 1	TATGACATTGGCATGGCTTGTGGG	TGCTCAGCCACATTCTCTGAGGTT	0.0101	0.0075	0.184
ACADL	Acyl-Coenzyme A dehydrogenase, long chain	TCCCTGAGATGGTTGCAGGGAAAT	TGAGCCTCTCGATTTGTGATCGCT	0.0011	0.0036	0.002*
SCPX	Sterol carrier protein 2	AAGCTGCCAGAAGGTGCTATGAGA	ACACAGTCCCAGTGCTTCGTAAGT	0.0113	0.0085	0.439
ACOX1	Acyl-Coenzyme A oxidase 1, palmitoyl	ATGCTGATGAAACATGCCCAGGTG	TTCAGACTGATGCCTCACAGCACT	0.0164	0.0185	0.399
ACOX2	Acyl-Coenzyme A oxidase 2, branched chain	ATGGTGGTGGTACGTGTGGACATT	AGAGTTTCTGCTGCTGAGTCTGGT	0.0044	0.0090	0.008**
PGAR	Angiopoietin-like 4	AAATGACCTCAGATGGAGGCTGGA	TCTCCGAAGCCATCTTTGTAGGCT	0.0192	0.0172	0.905
FABP4	Fatty acid binding protein 4, adipocyte	CGTGGGCTTTGCTACCAG	TGGTTGATTTTCCATCCCAG	0.0669	0.4280	0.002**
CAP	Sorbin and SH3 domain containing 1	ACAAGAATCTGCACTGGGAAGGGA	CCATTGGCATCGGCATTGGAAGAA	0.0753	0.0810	0.488
PLIN1	Perilipin	TTGGGTGGTGTCGAGAAGATGGTA	AAGGTGTGTCTAGAGAGAGTGTTGGC	0.0572	0.0740	0.401
ADIPOQ	Adiponectin, C1Q and collagen domain containing	TGTGTATCGCTCAGCGTTCAGTGT	AGTAGAGTCCCGGAATGTTGCAGT	0.0001	0.0001	0.088
MMP1	Matrix metallopeptidase 1 (interstitial collagenase)	TCTGGAGCAATGTCACACCCTTGA	TCCGATACCTGCACCTGGTTGAAA	0.0001	0.0001	0.107
UCP1	Uncoupling protein 1 (mitochondrial, proton carrier)	TTCTGCACAACGGTTCTGTCCTCT	TTGGGCACACTTGTGTACTGTCCT	0.0003	0.0002	0.070
ILK	Integrin-linked kinase	ATCACACACTGGATGCCATACGGA	TCCAATGCAAACTTCACAGCCTGG	0.2481	0.1484	0.053
UBC	Ubiquitin C	ACGTCAAGGCCAAGATCCAAGACA	TTCACGAAGATCTGCATCCCACCT	0.1505	0.2632	0.450
PCK1	Phosphoenolpyruvate carboxykinase 1	AACCTGGCCATGATGAACCCTACT	TCTTCTGGATGGTCTTGATGGCGT	0.0233	0.0148	0.127
PCK2	Phosphoenolpyruvate carboxykinase 2	AAGAAGTGCTTTGCCCTTCGCATC	AGGTTTGTCTTGCCACAAGCACTG	0.0246	0.0372	0.406
GK	Glycerol kinase	TGGAAGGTGTGCCAATATCTGGGT	TCTCTGCCGAGTTTGTAAGCCACT	0.0145	0.0177	0.349
AQP7	Aquaporin 7	TGTCCTGGAAGAAGTTTCCCGTGT	ACAGTGTCATGTGGTCAGGAAGGT	0.3712	0.4460	0.363
β-Actin	Actin, beta	ACCTGCGAAGATGTGTCTGAACCA	TGTACACCGTCAGCTTGTTGAGGA

Expression shown as the mean of the normalized expression value of each gene within the low- and high-marbled groups.

p value was calculated using analysis of variance.

Significant effects on intramuscular fat contents (p<0.05, 0.01, respectively).

The qRT-PCR was performed using the 2× Power SYBR Green PCR Master Mix (Applied Biosystems, Foster City, CA, USA) with the 7500 Real-Time PCR System (Applied Biosystems, USA) using 10 pM of each primer (Table 2). The PCR was performed for 2 min at 50°C and 10 min at 95°C, followed by 40 cycles of 95°C for 10 s, and then 60°C for 1 min. Following amplification, melting curve analysis was conducted to verify the specificity of the reactions. The end point used in the qRT-PCR quantification, Ct, was defined as the PCR threshold cycle number. The ΔCt value was determined by subtracting the β-actin Ct value for each sample from the target Ct value. The gene expression stability value of the β-actin gene was less than 0.05, which met the stability requirement to be a housekeeping gene, confirming its suitability as an internal housekeeping gene in this experiment. The qRT-PCR data were used to calculate the normalized expression values (2−ΔCt) for the statistical analysis.

To examine the association between IMF content and gene expression levels, a statistical analysis was conducted using the analysis of variance (ANOVA) model. This resulted in the following equation: Fat = μ+expression+age +residual. In this equation, Fat denotes the IMF content (%) of each animal, μ is the overall mean, Expression is a normalized gene expression value (2−ΔCt), and Age is the slaughtering age (months) as a covariate. To determine the major patterns and relationships in the gene expression data, we performed a principal component analysis (PCA) of the genes.

Construction of the gene co-expression network and direct interaction relationships

In gene co-expression networks, we refer to nodes as genes whose degrees indicate the number of links connected to the node. We used expression values for 48 genes and evaluated pairwise correlations between the gene expression profiles of each pair of genes using Pearson’s correlation coefficients (denoted as r). The network-encoded gene co-expression data were recorded as binary information (connected = 1, unconnected = 0) using a “hard” threshold. The degree of a node is the number of connections or edges the node has to other nodes. The degree distribution of a network has a generalized power-law form of p(k) ~ k, which is the defining property of a scale-free network. We investigated hard thresholding with the power adjacency function and selected a power of tau (τ) = 0.65 using the weighted gene co-expression analysis (WGCNA) package in R (Langfelder and Horvath, 2008).

To characterize the gene co-expression network and the direct interaction relationships among the differentially expressed genes, we used the Pathway Studio software package (version 6.0, Stratagene, La Jolla, CA, USA) program (Nikitin et al., 2003). To identify the gene co-expression network, we used a p-value cut-off of 0.05, a Pearson’s correlation (r) cut-off of 0.65 according to the result of WGCNA, with the removal of 5% of the genes with the most stable expression. To characterize the overall gene co-expression network topology, we investigated the node degree (or connectivity), betweenness centrality (BC), and closeness centrality (CC) using Cytoscape (Shannon et al., 2003) The node degree is the number of connections or edges the node has to other nodes. A node with high BC has considerable influence over what flows in the network. The BC can serve as a global property, as it is a useful indicator for detecting bottlenecks in a network. After experimental validation, we also examined the relationships and the effects of regulations (positive, negative, and unknown) among the significantly differentially expressed genes based on data in the literature.

RESULTS AND DISCUSSION

Gene expression and statistical analysis using principal component analysis and analysis of variance

Ten steers with low IMF content levels (low-marbled group) and 10 steers with high IMF content levels (high-marbled group) were used in an ANOVA model to identify an association between the expression of mRNA and the IMF content. The average IMF content in m. longissimus tissue for the low-marbled group was 7.44±0.46%; this value for the high-marbled group was 23.54±0.88% (Table 1). As shown in Table 1, the IMF content (p<0.01) was significantly different between groups.

To identify the genes related to marbling, we investigated 48 genes within the PPAR signaling pathway from the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway database (bta03320). Table 2 shows the gene expression levels of 48 genes. The qRT-PCR results showed that 30 genes had higher expression levels in the high-marbled group than in the low-marbled group. The effect of 48 genes on the IMF was estimated using an ANOVA model; the least-square means were also estimated using the Student’s t-test to determine if there was a significant difference between the marbled groups. The expression levels of 10 genes: retinoid X receptor alpha (RXRA), peroxisome proliferator-activated receptor gamma (PPARG), phospholipid transfer protein (PLTP), stearoyl-CoA desaturase (SCD), nuclear receptor subfamily 1 group H member 3 (NR1H3), fatty acid binding protein 3 (FABP3), carnitine palmitoyltransferase II (CPT2), acyl-Coenzyme A dehydrogenase long chain (ACADL), acyl-Coenzyme A oxidase 2 branched chain (ACOX2), and fatty acid binding protein 4 (FABP4), showed significant differences in gene expression between the low- and high-marbled groups (p<0.05) (Table 2). Of these, nine genes showed higher expression in the high-marbled group, and one gene, FABP3, showed a low gene expression pattern in the high-marbled group.

FABP3 is involved in the uptake, transport, and metabolism of fatty acids (FAs) and is mainly expressed in cardiac and skeletal muscle. Ovilo et al. (2002) reported that the FABP3 gene had important effects on the IMF content. Other studies have also revealed that FABP3 is positively associated with the IMF content (Li et al., 2007; Zhao et al., 2009). However, Cui et al. (2012) reported that the expression level of the FABP3 gene in two chicken breeds was lower in later growth stages. Higher IMF and triglyceride (TG) levels are deposited later in growth compared with earlier growth stages in the pectoralis major tissue. Recently, Ye et al. (2014) also reported that the FABP3 gene was down-regulated in higher IMF content muscle. In the present study, the FABP3 gene expression pattern was lower in the high-marbled group than in the low-marbled group of Hanwoo cattle. We therefore suggest that lower expression of the FABP3 gene in the longissimus muscle is not a result of fat deposition and is most likely regulated by other biological factors. However, further investigations are needed to confirm the relationship between the FABP3 gene and marbling traits in skeletal muscle and adipocytes.

The PCA is a useful tool for data simplification and for the visualization of relationships. Therefore, we applied PCA to the gene expression dataset. The first two principal components (PCs) explained approximately 73.3% of the total variance, allowing most of the information to be visualized in two dimensions. The PCA indicated that the most important pattern of gene expression (PC1, accounting for 60% of the variance in the data) was associated with differences in the IMF levels. Individual samples were clearly partitioned into two separate groups: high- and low-marbled groups based on PC1. In this analysis, the first PC illustrated a link between the nine differentially expressed genes, which had a positive relationship according to PC1 (Figure 1). Whereas, the FABP3 gene had a negative relationship according to PC1. Our experimental results suggested that these genes warrant further investigation as metabolic indicators of marbling in the bovine model.

Figure 1

Analysis of gene expression data by principal component analysis. Biplot of the first two principal components. L and H represent low-and high-marbled samples, respectively.

Construction of a gene co-expression network based on gene expression patterns

The nodes represent 48 genes within the PPAR signaling pathway, and the links between the nodes represent the associations between expression profiles across all samples. The absolute value of Pearson’s correlation coefficient was calculated for all pairwise comparisons. This correlation matrix was transformed into a matrix of adjacency using a “hard” threshold (τ). When the hard threshold is 0.65, the plot between log10 (p(k)) and log10 (k) shows that the network mostly follows a scale-free topology (R2 = 0.54, slope = −0.45). To find clusters of highly correlated genes, we used the Pathway Studio program. We found 18 correlated genes and identified four distinct gene co-expression sub-clusters (Figure 2). The overall trend observed in the large co-expression pattern was also in agreement with the PCA result showing 10 significantly differentially expressed genes. From the PCA result, the significantly up-regulated genes according to PC1 could be divided into two clusters. The first cluster contained RXRA, PLTP, FABP4, and PPARG. The second contained ACOX2, CPT2, SCD, NR1H3, and ACADL. The gene co-expression network also showed a similar pattern. The genes in the PPARG-FABP4-PLTP group were interrelated. The ACOX2-CPT2 group also showed connections. This suggested that GK plays a role in connecting the PPARG, FABP4, PLTP, ACOX2, and CPT2 groups. The main cluster was large with a mean degree of 2.8, and it consisted of 10 genes, all of which were up-regulated in the high-marbled group. This cluster was connected by nine target genes, including one transcription factor, PPARG. The nine target genes were acyl-CoA synthetase long-chain family member 6 (ACSL6), apolipoprotein A-I (APOA1), PLTP, GK, CPT2, aquaporin 7 (AQP7), acyl-Coenzyme A oxidase 1 palmitoyl (ACOX1), ACOX2, and FABP4. From the gene correlations, the FABP4-PPARG pair was highly correlated in the cluster (r = 0.92). PPARG and FABP4 are reportedly key regulators in the PPAR signaling pathway associated with adipogenesis and marbling traits in beef cattle. Therefore, the expression levels of these genes in muscle are reported to be positively correlated with marbling (Michal et al., 2006; Ahn et al., 2014). This is consistent with the fact that PPARG is directly linked to an increase in FABP4 among its targets or coactivators.

Figure 2

Gene co-expression network of the 48 genes in the peroxisome proliferator-activated receptor pathway based on their expression values. To assess correlations, we calculated the Pearson’s correlation coefficient (r) between all gene pairs. Edges join genes with co-expression correlation (Pearson’s r) ≥0.65. The red members indicate significantly differentially expressed genes from the analysis of variance test.

From the network topology analysis, the GK gene has also been shown to mediate the expression levels of genes involved in the network with the largest BC and CC values (BC = 0.35, CC = 0.69). The BC value is an indicator of a global central node. The effect of removing nodes with large BC values is similar to that of removing hub nodes because large BC nodes are correlated with hub nodes. The BC also implies that a site is relatively central with regard to all other sites. Thus, such sites are advantageously located to act as intermediaries. Therefore, we investigated the communication between genes and confirmed that hub and large-BC genes are located in the topological center of the network by calculating the CC in the network. The GK is directly connected to the following genes: ACSL6, ACOX1, ACOX2, AQP7, and PPARG. It was found that GK is functional in muscle cells and contributes to intramuscular TG synthesis for fat storage. Rahib et al. (2009) suggested that GK affects transcription factors that are important in muscle differentiation from a transcriptomic analysis of GK in mouse skeletal muscle. This gene showed no significant differences in expression between the high- and low-marbled groups. However, it may operate at the interface between fat and carbohydrate metabolism. Moreover, its deletion leads to alterations in the expression levels of genes involved in lipid and carbohydrate metabolism, PPARG, AQP7, or fatty acid synthase (Rahib et al., 2007). Therefore, GK may contribute to the cooperation between various genes and PPARG as a master regulator of fat and glycerol metabolism within the PPAR pathway.

The functional relationships between the significantly differentially expressed genes

As a candidate pathway involved in fat metabolism and development, KEGG analysis indicated that the PPARs of the PPAR signaling pathway in cattle trigger the induction of several downstream biological pathways, such as the pathways for lipid metabolism, adipocyte differentiation, and gluconeogenesis. Our results showed that eight genes were involved in the downstream pathways as transcription regulators, with the exception of PPARG and RXRA. Therefore, we categorized the significantly differentially expressed genes into the downstream pathways and investigated the relationships between these genes.

Lipid metabolism

The following genes were categorized into the lipid metabolism pathway: PLTP (lipid transport); SCD (lipogenesis); liver X receptor α, nuclear receptor subfamily 1, group H, member (LR1H3) (cholesterol metabolism); FABP3 (fatty acid transport); CPT2, ACADL, and ACOX2 (fatty acid oxidation); and FABP4 (adipocyte differentiation). These genes may be responsible for the high IMF content of Hanwoo.

The liver X receptor α (LR1H3 or LXRα) is a member of the nuclear receptor family of transcription factors and is involved in controlling cholesterol and fatty acid metabolism. This process is governed by sterol regulatory element-binding protein-1c gene (SREBP-1c), which regulates the acetyl-CoA carboxylase, fatty acid synthase, and SCD genes involved in the pathway of fatty acid biosynthesis (lipogenesis). The synergy between LR1H3 and SREBP-1c increases the synthesis of TG, which accumulates as cytoplasmic fat droplets and thus increases their secretion by PLTP into plasma in very-low-density lipoproteins (VLDL) (Repa et al., 2000). VLDL-TG production is directly related to free fatty acids, which are taken up by hepatocytes or myocytes for re-esterification or metabolism through β-oxidation occurring in the mitochondria or in peroxisomes. Recent reports have also shown that activated LR1H3 stimulates adipocyte differentiation through the induction of PPARG expression (Seo et al., 2004). These reports demonstrated that the pathway of the liver X receptor α (LxRα) gene through the induction of PPARG was related to fat metabolism processes such as TG formation and fatty acid oxidation. Earlier reports have also shown that the LR1H3 (Yu et al., 2006), SCD (Cho et al., 2011), and PLTP (Seong et al., 2013) genes are related to fat metabolism and IMF. Interestingly, these findings are also in agreement with our results, in which three genes were found to be up-regulated in the high-marbled group as opposed to the low-marbled group, and expression levels were positively related to the IMF content in Hanwoo (Table 2). Moreover, we found direct regulatory relationships among these genes (Figure 3). The LR1H3-SCD and LR1H3-PLTP pairs showed a direct interactive relationship with the significantly differentially expressed genes identified in our study. Liver X receptor (LXR) ligands induce the expression of the PLTP gene and then regulate plasma lipoprotein metabolism (Laffitte et al., 2003). The lipogenic gene SCD is regulated partly by the insulin-sensitive transcription factor LxRa (Hebbachi et al., 2008). Other relationships, such as those between the PPARG-RXRA, PPARG-SCD, and RXRA-FABP4 pairs, have been thoroughly investigated in many studies (Chan et al., 2009).

Figure 3

The direct relationships between the 10 significantly differentially expressed genes. Seven of the genes had direct relationships based on the literature or our expression studies. Genes were classified according to function (regulation, promoter binding, and expression).

In mammals, the oxidation of fatty acids occurs in two major systems: the mitochondrial and the peroxisomal β-oxidation systems. The first step of the β-oxidation cycle is catalyzed by acyl-coenzyme A dehydrogenases and acyl-CoA oxidases in the two respective systems. Prior to oxidation, fatty acids are transported into the mitochondria by carnitine palmitoyl transferases, they then enter the tricarboxylic acid cycle, where they are used by the electron transport chain to produce adenosine triphosphate (ATP). TG is the major component of marbled fat in skeletal muscle and is found mostly within intramuscular adipocytes. During the TG formation process, the first activation reaction is the consumption of two molar equivalents of ATP in a process that demands a high level of energy. Interestingly, the chemical inhibition of the mitochondrial respiratory chain, which produces ATP from products of the citric acid cycle and fatty acid oxidation, leads to reduced fat accumulation in the adipocytes of Caenorhabditis elegans, indicating that the mitochondria play a key role in fat synthesis (McKay et al., 2003). Acyl-coenzyme A dehydrogenase encodes one of the adipocyte fatty acid-binding proteins that belongs to a superfamily of lipid-binding proteins that are involved in the development of fatness traits. Recently, Jiang et al. (2014) found that ACADL genes were positively correlated with the adipocyte volume. We also found that the expression levels of the CPT2, ACADL, and ACOX2 genes were increased in the high-marbled group in the present study. Moreover, our previous reports found that gene expression within the mitochondrial oxidative phosphorylation pathway, involving NADH dehydrogenase 1, NADH dehydrogenase 2, cytochrome oxidase III, and ATP synthase F0 subunit 6 (ATP6), strongly correlated with IMF content (Kim et al., 2008). Taken together, these findings indicate that fatty acid oxidative metabolism in skeletal muscle is positively correlated with IMF content and fatty acid use for TG formation by providing the ATP.

Adipocyte differentiation

FABP4 is involved in adipocyte differentiation in the PPAR signaling pathway. FABP4 has been examined in many studies with regard to its function in fatness traits. Jurie et al. (2007) reported that TG content is also strongly correlated with FABP4 activity and oxidative enzyme activities at different IMF levels in cattle. Intracellular fatty acids in intramuscular adipocytes bind to FABP4, facilitating the transport of fatty acids from the plasma membrane to sites of fatty acid oxidation or esterification into the TG. In the present study, we found that the FABP4 gene showed a higher level of gene expression in the high-marbled sample, with the level differing significantly among the groups (Table 2). According to these results, we suggest that more marbled tissue exhibits more FABP4 content and uses more fatty acid for TG formation.

In conclusion, the 48 gene expression profiles examined in this study revealed 10 genes in the PPAR signaling pathway that were differently expressed between the high- and low-marbled groups, indicating a significant association with marbling. From the gene co-expression relationships and a network topology analysis, the main network was identified as 10 genes with a gene connectivity value of 2.8, all of which were up-regulated in the high-marbled group. Finally, we investigated the functional relationships according to the downstream pathways. We found that the major downstream pathway associated with the differently expressed genes was fatty acid oxidation. Activation of PPARG promotes terminal differentiation through the induction of a range of genes important for TG formation through fatty acid oxidation among downstream pathways involving CPT2, ACADL, and ACOX2.