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 . The skeletal muscle genes within the glycolytic pathway, which include protein kinase adenosine monophosphateactivated 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.

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).

Quantitative real-time polymerase chain reaction and statistical analysis
The mRNA levels were analyzed by quantitative realtime polymerase chain reaction (qRT-PCR) using genespecific 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.
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  ). The network-encoded gene coexpression 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 -r , 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 coexpression 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.

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 (highmarbled 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 lowmarbled 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.

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 (R 2 = 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 upregulated 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.
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.
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 , 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).
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 acidbinding proteins that belongs to a superfamily of lipidbinding 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 highand 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, 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). 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.