Skip to main content


Prediction of hub genes associated with intramuscular fat content in Nelore cattle

Article metrics



The aim of this study was to use transcriptome RNA-Seq data from longissimus thoracis muscle of uncastrated Nelore males to identify hub genes based on co-expression network obtained from differentially expressed genes (DEGs) associated with intramuscular fat content.


A total of 30 transcriptomics datasets (RNA-Seq) obtained from longissimus thoracis muscle were selected based on the phenotypic value of divergent intramuscular fat content: 15 with the highest intramuscular fat content (HIF) and 15 with the lowest intramuscular fat content (LIF). The transcriptomics datasets were aligned with a reference genome and 65 differentially expressed genes (DEGs) were identified, including 21 upregulated and 44 downregulated genes in HIF animals. The normalized count data from DEGs was then used for co-expression network construction. From the co-expression network, four modules were identified. The topological properties of the network were analyzed; those genes engaging in the most interactions (maximal clique centrality method) with other DEGs were predicted to be hub genes (PDE4D, KLHL30 and IL1RAP), which consequently may play a role in cellular and/or systemic lipid biology in Nelore cattle. Top modules screened from the gene co-expression network were identify. The two candidate modules had clear associated biological pathways related to fat development, cell adhesion, and muscle differentiation, immune system, among others. The hub genes belonged in top modules and were downregulated in HIF animals. PDE4D and IL1RAP have known effects on lipid metabolism and the immune system through the regulation of cAMP signaling. Given that cAMP is known to play a role in lipid systems, PDE4D and IL1RAP downregulation may contribute to increased levels of intracellular cAMP and thus may have effects on IF content differences in Nelore cattle. KLHL30 may have effects on muscle metabolism. Klhl protein families play a role in protein degradation. However, the downregulation of this gene and its role in lipid metabolism has not yet been clarified.


The results reported in this study indicate candidate genes and molecular mechanisms involved in IF content difference in Nelore cattle.


Intramuscular fat (IF) content is an important meat quality trait that represents the amount of fat accumulated between muscle fibers or inside muscle cells; it is the sum of phospholipids and triglycerides. IF plays a key role in meat quality, as it is mainly responsible for the palatability, tenderness, and nutritional value of meat [1,2,3]. IF content may be measured via chemical methods (which evaluate total lipids) or by marbling scores (which may be assessed after slaughter or by ultrasonography). IF content is a polygenic trait regulated by many genes involved directly or indirectly in fat metabolism [4].

Next generation sequencing applied to tissue transcriptomes (RNA-Seq) is an approach for screening the expression of functional candidate genes and consequently is possible identify important molecular mechanisms that generate variation in pathways resulting in different tissue phenotypes [5, 6]. This technology has been used to identify differentially expressed genes (DEGs) in muscle tissue of Nelore cattle divergent to different meat traits, as such: tenderness [7], fatty acid composition [8], intramuscular fat content (measured by marbling scores) [9], residual feed intake [10] and iron content [11].

Using RNA-seq (mRNA sequencing) data, Cesar et al. [9] identified differentially expressed genes (DEGs) in muscle tissue and elucidated some of the molecular mechanisms involved in the lipid metabolism in Nelore steers breeding herd (castrated males, belonging a experimental population) genetically divergent for intramuscular fat (measured by marbling scores). When castrated and uncastrated Nelore males are compared, there are differences in intramuscular fat deposition [12, 13]. Therefore, gene expression and the molecular mechanisms regulating differences in intramuscular fat deposition may be different in castrated and uncastrated animals. Most commercially Nelore cattle produced in Brazil and exported to European courtiers are uncastrated. Therefore, understanding the biological and functional mechanisms that regulate IF content in uncastrated Nelore cattle is a compelling question in meat science, since this knowledge are still unclear.

The development of IF is influenced by different genes and a complex network of genetic interactions. The biological interpretation of the results obtained from the analysis of RNA-seq, especially for DEGs, remains a challenge. Additional studies are necessary, since many genes and mechanisms that induce differences in uncastrated Nelore male IF content are still unknown, especially hubs genes (biomarkers). Hub genes are highly correlated with a large number of genes, have been shown to play key regulatory roles in gene expression networks [14,15,16], and are proposed to play an important role in the overall biology of organisms [17]. Network approaches have been used to identify complex transcriptional regulation, i.e., in the identification of hub genes. Thus, co-expression analysis may facilitate the detection of important biological pathways involved in targeted phenotypes. Consequently, the aim of this study was to use transcriptome RNA-Seq data from longissimus thoracis muscle of uncastrated Nelore males to identify hub genes based on a gene co-expression network constructed from DEGs associated with IF content.


Sample collection and phenotype

All animals (N = 80), uncastrated Nelore males belonging to the same contemporary group (i.e., remaining together from birth until slaughter), were from the Capivara Farm, which participates in the Nelore Qualitas Breeding Program. They were reared on grazing systems (Brachiaria sp. and Panicum sp. forage and free access to mineral salt) and finished in confinement for approximately 90 days. The diet was based on whole-plant silage and mix of sorghum grain, soybean meal or sunflower seeds were used as concentrate, with a concentrate/roughage ratio from 50/50 to 70/30. The animals were slaughtered at a mean age of 24 months - all on the same day and under the same conditions.

Longissimus thoracis muscle samples were collected from an area between the 12th and 13th ribs of the left half of each carcass two times: 1) at slaughter, stored in 15-mL Falcon tubes containing 5 mL RNA holder (BioAgency, São Paulo, SP, Brazil) at − 80 °C until total RNA extraction was performed for RNA sequencing analysis; and 2) 24 h after slaughter for an analysis of the IF content.

IF content was quantified for all animals, through chemical method for total lipid content, according to the methodology described by Bligh & Dyer [18]. The animals were ranked in accordance with their IF phenotype and the samples derived from the animals with the 15 highest and 15 lowest values for IF content were selected for RNA-seq analysis (Table 1). A Student’s t-test was performed to evaluate whether there were significant differences in IF between the selected groups.

Table 1 Descriptive statistics for ribeye muscle area and intramuscular fat content of Nelore cattle

RNA-seq library construction

Total RNA was isolated from the longissimus thoracis samples (an average of 50 mg each) using the RNeasy Lipid Tissue Mini Kit (Qiagen, Valencia, CA, USA) according to the manufacturer’s protocol. The following three methods were used for RNA quantification and qualification: RNA purity was determined by evaluating absorbance using a NanoDrop 1000 spectrophotometer (Thermo Fisher Scientific, Santa Clara, CA, USA; 2007); RNA concentration was measured using a Qubit® 2.0 Fluorometer (Invitrogen, Carlsbad, CA, USA; 2010), and RNA integrity was assessed using an RNA Nano 6000 Assay Kit with the Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA; 2009).

Sequencing libraries were prepared using the TruSeq RNA Sample Preparation Kit® (Illumina, San Diego, CA) following the manufacturer’s protocol. Libraries were pooled to enable multiplexed sequencing and generated an average of approximately 25 M reads per sample. RNA sequencing (RNA-Seq) was carried out on a HiSeq 2500 System (Illumina®) that generated 100 bp paired-end reads.

Data filtering and alignment of reads

Trimmed data (trimmed reads) were obtained by removing low quality reads (adapter sequence and reads containing poly-N) from raw data using Trimmomatic v.0.36 [19]. All downstream analysis was based on the trimmed data with high quality reads. HISAT2 v.2.0.5 [20] was used to align the paired-end trimmed reads to the bovine reference genome (UMD3.1.1 Bos taurus) and chromosome Y (Btau 4.6.1), both deposited in National Center for Biotechnology Information (NCBI) (

Differentially expressed genes (DEGs) and enrichment analysis

The Cufflinks2 v.2.1.1 suite of tools [21] was used for transcriptome assembly and differential expression analysis. Cufflinks2 assembles transcriptomes from RNA-Seq data and quantifies their expression in fragments per kilobase of transcript per million reads mapped (FPKM). After assembly, the Cufflinks2 output per sample was merged together by Cuffmerge2. A uniform set of transcripts was obtained for all samples and then Cuffdiff2 was used to test for genes that were differentially expressed between the IF groups. The Cuffdiff2 uses the t-test to calculate the p-values for differentially expressed genes, based on the normalized FPKM values (log2FPKM) between two conditions (HIF and LIF groups). The method assumes that the data present normal distribution [22]. False discovery rates (FDR) were controlled using the Benjamini-Hochberg procedure in which we considered a gene to be differentially expressed if it had a q-value ≤ 0.05 [23].

The enrichment and pathway analyses of DEG sets was performed using the Database for Annotation, Visualization, and Integrated Discovery (DAVID 6.8) [24]. The DAVID Pathway was used to map the enriched pathways from the Kyoto Encyclopedia of Genes and Genomes (KEGG) database [25]. FDR were controlled by the Benjamini-Hochberg method [23] considering a p-value of less than 10% (EASE score ≤ 0.1).

Co-expression network analysis and prediction of hub genes

All analyses were performed with plugins or applications from Cytoscape v.3.7, a free software package for visualizing, modeling, and analyzing the integration of biomolecular interaction networks with high-throughput expression data and other molecular states [26]. The co-expression network was constructed from the ExpressionCorrelation plugin [27]. The similarity matrix of count data from DEGs (normalized by FPKM) was computed using a Pearson correlation. A histogram tool was used for the screening criteria at node score cut-offs > 0.75 and < − 0.75 and employed to identify statistical significance of the pairwise correlations. We established a co-expression network between the significantly co-expressed DEGs.

To identify network modules, co-expression network was analyzed using Reactome FI plugin (cluster FI Network function) [28]. The modules within the network were mapped onto pathways and biological processes using the same plugin with analyze module function [28]. FDR were controlled by the Benjamini-Hochberg method [23] considering a p-value of less than 10%.

The cytoHubba plugin [29] was used to explore co-expression network nodes (hub genes). This plugin provides a user-friendly interface to explore important nodes in biological networks and computes using eleven methods (degree, edge percolated component, maximum neighborhood component, density of maximum neighborhood component, maximal clique centrality -MCC and six centralities - bottleneck, EcCentricity, closeness, radiality, betweenness, and stress). All these parameters can indicate the robustness of the analysis [29, 30]. MCC method has a better performance on the precision of predicting essential genes from the co-expression network [29] and to generate a subnetwork [30]. The top three hub genes were chosen to define top modules screened from the gene co-expression network. The top modules could be functionally relevant and therefore is used to assess biological function [30]. The Reactome FI Plugin [28] was then used to determine top modules- distinct groups of genes interactions [31] and to analyze the functional enrichment of biological pathways within the of each top module.

Results and discussion

RNA-seq data alignment

After quality control of the raw reads (approximately 25 million reads), the mean number of reads per sample (paired-end) was approximately 21.9 million. In the library, 88% of the clean reads were uniquely mapped to bovine reference genome UMD3.1.1 Bos taurus and chrY of Btau 4.6.1. The mean number of reads mapped in pairs using Hisat2 software was approximately 20.7 million (96.44%) with 45× sequencing coverage (coverage for all transcripts of all samples). The box plot containing the transformed FPKM values (log10) for each group and the plot of principal component analysis (PCA) were showed in Additional file 2: Figure S1. These plots were constructed using the cummerRbund package [32]. The distribution of quartiles, on box plot, was consistent between groups, indicating high quality of the data. In addition, the medians were similar in the two groups and close to − 1, indicating that the level of sequencing coverage permitted the identification of low-expressed genes [10]. PCA showed the formation of different groups (highest and lowest IF), indicating differences in the expression of genes between the HIF and LIF groups.

Differentially expressed genes (DEGs) and enrichment analysis

A total of 65 DEGs (q ≤ 0.05) were found, including 21 upregulated genes and 44 downregulated genes (Table 2) in the HIF group. The genes were distributed across almost all bovine chromosomes in different proportions, except chromosomes 11, 14, 24, mitochondrial (MT), and Y. Important genes that play a role in lipid metabolism, muscle development and immune systems were differentially expressed. DEGs associated with lipid metabolism included: CDP-diacylglycerol synthase 1 (CDS1) and solute carrier superfamily genes (SLC22A4 or OCTN1 and SLC2A3 or GLUT3), all of which were upregulated in HIF animals. Also were found regulatory genes: mediator complex subunit 25 (MED25) and forkhead box O1 (FOXO1), both downregulated in HIF animals.

Table 2 Differentially expressed genes in the longissimus thoracis muscle of Nelore cattle divergent for intramuscular fat content phenotypes

CDS1, plays an important role in the synthesis of triglycerides. Qi et al. [33] elucidated the role of CDS, a key enzyme in phospholipid metabolism, in cellular lipid storage as well as in the differentiation of adipocytes. They demonstrated that CDS1 may regulate the expansion of lipid droplets essential for adipogenesis. The CDS1 up-regulation in HIF animals could contribute in the adipocytes differentiation of these animals.

SLC22A4 (OCTN1) is a plasma membrane carnitine transporter (an organic cation transporter) which is an important component of the carnitine system. Adequate carnitine levels are required for normal lipid metabolism and are important for energy metabolism [34, 35]. SLC2A3 (GLUT3), was classically defined as a neuronal glucose transporter due to its high level of expression and initial characterization in nervous tissue [36]. However, SLC2A3 is also highly expressed in other tissues with high energy needs, including muscle [37, 38]. The animals in the HIF group had higher expression levels of genes associated with lipid metabolism, indicating that these genes may contribute to adipogenesis and the maintenance of lipid metabolism in HIF animals.

A few notable examples of transcription factors that regulate lipid metabolism include FoxOs proteins and hepatic nuclear factors (HNFs) [39]. The FOXO1 gene (belonging to the FoxOs proteins family), is one of the most important transcriptional effectors, promoting the expression of genes involved in gluconeogenesis [40, 41]. FOXO1 plays a role in hepatic lipid and lipoprotein pathways, potentially considered the central link for the regulation and coordination of insulin action on carbohydrate and lipid metabolism [42, 43]. Previous studies showed that lipogenesis is modulated by FOXO1 through sterol regulatory element binding protein 1c (SREBP-1c) [42, 43].

Functional annotation of DEGs revealed biological processes GO terms such as translation regulatory and transporter activities. The regulatory genes, such as MED25, we suggest that these genes could be involved in modulation of others gene expression and thus influencing in intramuscular fat deposition. MED25 is a transcription factor that belongs to a variable functional complex that controls the constitutive expression of genes. Rana et al. [44] showed that MED25 is a cofactor of HNF4α (nuclear factor), i.e., MED25 plays a vital role in the modulation of the transcriptional activity of HNF4α. HNF4α regulates genes that are responsible for lipid and drug metabolism, such as cytochrome P-450. These researchers also showed that down regulation of MED25 impairs a specific set of HNF4α target genes, suggesting a role for MED25 in metabolism and lipid homeostasis.

Gene: Symbol of the differentially expressed gene; Locus: location of the gene in the Bos taurus genome; HIF (μ): mean normalized counts from highest intramuscular fat content; LIF: mean normalized counts from lowest intramuscular fat content; Fold change (log2): FPKM values obtained for HIF and LIF; q-value: p-value adjusted.

GO terms significantly enriched in the biological processes and pathways categories for all DEGs found this study are shown in Additional file 1: Table A1. The GO terms were closely associated with cell regulation and may have a putative association with lipid and muscle metabolism, for example: Cellular response to fibroblast growth factor stimulus (GO:0044344). Two genes were related to this GO term: C-C motif chemokine ligand 2 (CCL2) and delta-like canonical notch ligand 4 (DLL4), both of which were upregulated in HIF animals. CCL2 was associated with intramuscular adipocyte differentiation in Japanese Black cattle [45] and DLL4 was associated with beef tenderness induced by acute stress in Angus cattle [46].

GO terms and pathways related to the immune and hormonal systems were also overrepresented, for example: antigen processing and presentation of peptide antigen via MHC class I (GO:0002474), immune response (GO:0006955) and insulin resistance (bta04931). Bola family member (BOLA), upregulated in HIF animals, was found in most GO terms significantly enriched in the biological processes category and pathways associated with the immune system (see Additional file 1: Table A1). BOLA has been associated with marbling in Hanwoo (Korean Cattle) [47], tenderness [7], and reproductive performance [48] in Nelore cattle.

Other gene found in the most GO terms related immune and hormonal systems was NFKB inhibitor alpha (NFKBIA), downregulated in HIF animals. Immune cells such as macrophages depend of the recognition of lipid ligands by membrane proteins, surface/extracellular, and intracellular immune receptors. Involvement of the lipid receptors triggers an immune response. Oxidized lipids activate nuclear receptors, which play a role in lipid homoeostasis and also regulate, for example, the immune response directed by NFKB. Activation of macrophages promotes the production of cytokines and the induction of acute phase response, accompanied by systemic lipid changes [49].

NFKBIA downregulation may contribute to intracellular the increase in 3′–5′-cyclic adenosine monophosphate (cAMP) binding levels. cAMP is known to have a significant role in adipogenesis [50] (other details will be discussed afterward) and thus may have effects on IF content differences in Nelore cattle, if cAMP plays a role in lipid systems. The NFKBIA gene was found in two GO terms and six pathways. Insulin resistance (bta04931) was among the pathways found. Insulin is synthesized by beta cells of islets of Langerhans in the pancreas and is the most important hormone in the regulation of energy metabolism. This hormone is essential for carbohydrate intake, protein synthesis, and fat storage [51]. Insulin resistance causes some disorders of lipid metabolism: increased triglyceride and decreased HDL levels [52]. This gene was previously reported to be associated with gain in crossbred beef steers [53].

Gene co-expression network analysis and prediction of hub genes

The expression levels of transcripts changed dynamically between HIF and LIF animals as shown in the heat map (Fig. 1). We computed the correlation matrix using DEG expression profiles (Fig. 2). Most DEGs have a moderately high correlation, positive or negative. We used the correlation matrix to construct the co-expression network (Fig. 3). The network has 274 significantly correlated gene pairs (edges) that were discovered for 54 genes (nodes). From the co-expression network, module analysis of the components of the associated pathways could be carried out [31]. The parameters that indicate the robustness of the analysis [29, 30] are shown in Additional file 1: Table A2. The modules (same colour code in Fig. 3a; listed in Additional file 1: Table A3). were mapped onto biological processes and pathways. Four modules were identified within the Network, which had clear associated biological pathways. These include biological processes and pathways related to fat development, muscle contraction, cell adhesion, and muscle differentiation, among others (Fig. 3a, for more details see Additional file 1: Table A4).

Fig. 1

Heat map displays differentially expressed genes found in Nelore cattle divergent for intramuscular fat content. The differing colors represent differing levels of expression of those genes

Fig. 2

Correlation matrix plot. The differing colors represent differing levels of correlation between genes

Fig. 3

Gene co-expression network construction and module analysis. a Gene co-expression network. The nodes represent the differentially expressed genes involved in four modules (different colour code), the lines represent the interaction between nodes and summary of enrichment analysis of the modules (same colour code). b Top three hub genes with a maximal clique centrality -MCC. The more forward ranking is represented by a redder color. c Top modules screened from the gene co-expression network. The top three hub genes were chosen to define the subnetwork (top modules), as this represents the most functional elements of the network

Since we are looking for hub genes within the co-expression network, the first three genes in the MCC method, predicted by CytoHubba plugin (Fig. 3b), were chosen to define the subnetwork top modules (Fig. 3c, for more details see Additional file 1: Table A5). Two modules were identified within the subnetwork, which had clear associated biological pathways (for more details see Additional file 1: Table A6). The two candidate modules to contained considerable amount of genes and represent high intensity of the connectedness within the module in order to maintain high precision in their predictive capability. This analysis allowed the identification of genes that are related to lipid metabolism. The genes found this study, may therefore could to play a role as biomarkers of bovine IF.

Three top hub genes were predicted in this study were: phosphodiesterase 4D (PDE4D), kelch-Like family member 30 (KLHL30), and interleukin 1 receptor accessory protein (IL1RAP) (Fig. 3b). Hub genes, highly interconnected with nodes in a network, have been shown to be functionally significant [54], that is, the hub genes found may have a putative effect on lipid metabolism in Nelore cattle. According to our results, the PDE4D gene exhibited the highest centrality, indicating that it has influence on the other DEGs identified in this study, as such NFKBIA. PDE4D was downregulated in HIF animals and was present in module 1 of the subnetwork (Additional file 1: Table A5). This gene participated in more than 20 GO terms biological process and pathways in accordance to reactome FI analysis (FDR < 0.1), as for example: cAMP signaling pathway, regulation of cAMP metabolic process and positive regulation of interleukin-5 production (Additional file 1: Table A6), revealing that this gene also has effects on the immune system of the animals studied.

PDE4D belongs to a family of four PDE4 genes, all encoding phosphodiesterase that hydrolyze the second messenger cAMP binding [55, 56]. The cAMP signaling pathways are among the well-characterized mechanisms controlling adipocyte differentiation [50] and immune system [57]. PDE inhibitors and synthetic cAMP analogs are commonly employed to switch on the adipogenic program in vitro [58]. PDE4 genes can be regulated by transcriptional and/or by post-translational mechanisms [59,60,61]. This reflects a way of negative feedback for cAMP signaling, phosphorylation and activation of PDE4 genes by the cAMP-dependent protein kinase (PKA) have been reported [59,60,61], that is, where PDE4 is inhibited, cAMP levels are increased [57, 62]. This mechanism could putatively allow the rapid regulation of PDE4D activity in selected regions within cells and the impact on cellular functions requiring local PDE4D mediated control of cAMP levels [59]. Elevation of cellular cAMP concentration could, for example, influence early adipocyte differentiation [63].

As mentioned earlier, to maintain homeostasis, cyclic nucleotide levels are regulated by PDEs, with PDE4s predominantly responsible for degradation of cAMP, this regulation could have effects in inhibition activity of other promoters such as NFKBIA. NFKBIA was differentially expressed this study (downregulated in HIF animals) and was present in subnetwork module 0 (Fig. 3c).

Such effects on NFKBIA gene cause increased expression of anti-inflammatory signals and decreased mRNA expression of cytokines and other inflammatory mediators, such as interleukins. Therefore, cAMP signaling supports the maintain immune homeostasis by modulating the production of pro-inflammatory and anti-inflammatory mediators [57]. When intracellular cAMP concentrations are heightened, inflammatory signaling is reduced [57]. The fact that PDE4D was downregulated in HIF animals may increase the intracellular cAMP through the regulation of other genes and thus may contribute to an increase in IF content and maintain immune homeostasis of Nelore cattle. PDE4D was previously related with marbling score in a commercial Hanwoo cattle population from genome-wide association study (GWAS) [64].

KLHL30 (downregulated in HIF animals) was identified as a hub gene and was present in module 1 of the subnetwork (Additional file 1: Table A5). This gene was not present in GO terms biological processes or pathways overrepresented this study. KLHL30 belonging to the kelch (Klhl) superfamily that consists of a large number of structurally and functionally diverse proteins characterized by the presence of a kelch-repeat domain [65].

Klhl proteins are involved in a number of cellular and molecular processes such as cell migration, cytoskeletal arrangement, regulation of cell morphology, protein degradation (ubiquitination process), and gene expression [66]. The ubiquitin-proteasome system degrades the major proteins of contractile skeletal muscle and plays an important role in muscle loss [67]. Although KLHL30 was not present in pathways and biological processes related to fat and the exact function that this gene has on production animal metabolism is still unknown, Piórkowska et al. [68] showed that KLHL30 was differentially expressed in muscle tissue of chickens selected for differential shear force. This result, together with ours, shows the putative pleiotropic activities of KLHL30 on qualitative traits in animals.

IL1RAP, downregulated in HIF, was identified as a hub gene encoding the interleukin 1 receptor accessory protein [69]. It is a necessary part of the interleukin 1 receptor complex which initiates signaling events that result in the activation of interleukin 1-responsive genes. Interleukins are a group of cytokines which play a role, primarily, in the immune system and in glucose/lipid metabolism [70]. Interleukins are expressed in a wide range of cells of the immune, neural and endocrine systems, reflecting the pleiotropic activities of this molecule [71]. The immune system can influence lipids and lipoprotein levels [72]. Lipids, besides being structural components of cellular membranes and serving as fuel stores, also play roles as effectors and second messengers associated with the immune system [73]. Interleukins can directly modulate lipid metabolism by increased of the activity of lipoprotein lipase, the enzyme that hydrolyzes triglycerides in lipoproteins. For example, the Interleukin-1 may inhibits adipocyte maturation and the synthesis of fatty acid transport proteins in adipose tissue in vitro [74].

IL1RAP was present in module 0 of the subnetwork (Additional file 1: Table A5). This gene participated in more than ten GO terms biological process and seven pathways in accordance to reactome FI analysis (FDR < 0.1), as for example: positive regulation of interleukins production, positive regulation of synapse assembly, protein complex assembly and inflammatory mediator regulation of TRP channels pathway (Additional file 1: Table A6). The elevation of cAMP can cause gene-specific inhibition of interleukin expressions and consequently of their receptors. The downregulation of IL1RAP in HIF animals may contribute to increase the intracellular cAMP and thus may to an increase in IF content and maintain immune homeostasis of Nelore cattle, as discussed before.

Of the hub genes identified from the coexpression analysis, only IL1RAP had colocalization with known QTLs controlling mainly reproductive traits of bovine (17 traits were associated, with 18 QTL/association found). This result, shows the putative pleiotropic activities of IL1RAP in producing qualitative traits in bovine.

We identify hub genes that play a role in lipid and muscle metabolism and immune system. These genes were the main regulators of others DEGs found this study and this promotes a better understanding of the essential biological mechanisms involved in IF deposition of Nelore cattle.


Our results show that muscle cells of Nelore cattle phenotypically divergent for IF expressed genes related to lipid and muscle metabolism, as well as genes related to the immune system. Through gene co-expression network analysis, we identify modules and hub genes, which are the main regulators of other DEGs identified and consequently may play a role in Nelore cattle’s cellular or systemic lipid biology. We show that the fact that PDE4D was downregulated in HIF animals may increase the intracellular cAMP, and this, could influence adipocyte differentiation. The cAMP regulation could have effects in inhibition activity of other promoters such as NFKBIA., and thus may contribute to an increase in IF content and maintain immune homeostasis of Nelore cattle. Such effects on NFKBIA gene cause increased expression of anti-inflammatory signals and decreased mRNA expression of cytokines and other inflammatory mediators, such as interleukins (IL1RAP was a hub gene and was downregulated in HIF, this gene was member of interleukins complex). Since, the interleukins can directly modulate lipid metabolism by increased of the activity of lipoprotein lipase, the downregulation of IL1RAP in HIF animals may decrease of the activity of lipoprotein lipase and thus may contribute to an increase in IF content in animals. This study identified potential biomarkers and molecular mechanisms involved in IF content difference in Nelore cattle.

Availability of data and materials

The dataset utilized in this study belongs to a Qualitas Nelore breeding program company, and could be available on request. The author does not have authorization to share the data.



BARX homeobox 2


BolA family member

CDS1 :

CDP-diacylglycerol synthase 1


CLOCK interacting pacemaker

COL18A1 :

Collagen type XVIII alpha 1 chain


Differentially expressed genes

EGR1 :

Early growth response 1


False discovery rates


Forkhead box O1


Forkhead O box proteins


Fragments per kilobase of transcript per million reads mapped


GlcNAc-1-phosphotransferase subunits alpha/beta


Gene ontology


Genome-wide association study


Highest intramuscular fat content

IER2 :

Immediate early response 2


Intramuscular fat


Interleukin 1 receptor accessory protein

KLHL30 :

Kelch-Like family member 30


Lowest intramuscular fat content


Maximal clique centrality

MED25 :

Mediator complex subunit 25


Myosin binding protein H


Nuclear factor and interleukin 3 regulated


NFKB inhibitor alpha


Phosphodiesterase 4D

SLC22A4 and SLC2A3 :

solute carrier superfamily genes


  1. 1.

    Wood JD, Enser M, Fisher AV, Nute GR, Sheard PR, Richardson RI, et al. Fat deposition, fatty acid composition and meat quality: a review. Meat Sci. 2008;78(4):343–58.

  2. 2.

    Pickworth CL, Loerch SC, Velleman SG, Pate JL, Poole DH, Fluharty FL. Adipogenic differentiation state-specific gene expression as related to bovine carcass adiposity. J Ani Sci. 2010;89(2):355–66.

  3. 3.

    Costa ASH, Pires VMR, Fontes CMGA, Prates JAM. Expression of genes controlling fat deposition in two genetically diverse beef cattle breeds fed high or low silage diets. BMC Vet Res. 2013;9(118):1–16.

  4. 4.

    Michal JJ, Zhang ZW, Gaskins CT, Jiang Z. The bovine fatty acid binding protein 4 gene is significantly associated with marbling and subcutaneous fat depth in wagyu x Limousin F2 crosses. Anim Genet. 2006;37(4):400–2.

  5. 5.

    Trapnell C, Williams BA, Pertea G, Mortazavi A, Kwan G, Van Baren MJ, et al. Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation. Nature Biotech. 2010;28(5):511–5.

  6. 6.

    He H, Liu X. Characterization of transcriptional complexity during longissimus muscle development in bovines using high-throughput sequencing. PLoS One. 2013;8:e64356.

  7. 7.

    Fonseca LFS, Gimenez DFJ, Silva DBS, Barthelson R, Baldi F, Ferro JA, et al. Differences in global gene expression in muscle tissue of Nellore cattle with divergent meat tenderness. BMC Genomics. 2017;18(945):1–12.

  8. 8.

    Berton MP, Fonseca LFS, Gimenez DFJ, Utembergue BL, Cesar ASM, Coutinho LL, et al. Gene expression profile of intramuscular muscle in Nellore cattle with extreme values of fatty acid. BMC Genomics. 2016;17:972.

  9. 9.

    Cesar AS, Regitano LC, Koltes JE, Fritz-Waters ER, Lanna DP, Gasparin G, et al. Putative regulatory factors associated with intramuscular fat content. PLoS One. 2015;10(6):1–21.

  10. 10.

    Tizioto PC, Coutinho LL, Oliveira PSN, Cesar AS, Diniz WJ, Lima AO, et al. Gene expression differences in longissimus muscle of Nelore steers genetically divergent for residual feed intake. Sci Rep. 2016;6:39493.

  11. 11.

    Diniz WJ, Coutinho LL, Tizioto PC, Cesar AS, Gromboni CF, Nogueira AR, et al. Iron content affects Lipogenic gene expression in the muscle of Nelore beef cattle. PLoS One. 2016;11(8):e0161160.

  12. 12.

    Vaz FN, Restle J, Feijó GLD, Brondani IL, Rosa JRP, Santos AP. Qualidade e composição química da carne de bovinos de corte inteiros ou castrados de diferentes grupos genéticos Charolês x Nelore. Rev Bras Zootec. 2001;30(2):518–25.

  13. 13.

    Ribeiro EL, Hernandez JA, Zanella EL, Shimokomaki M, Prudêncio-Ferreira SH, Youssef E, et al. Growth and carcass characteristics of pasture fed LHRH immunocastrated, castrated and intact Bos indicus bulls. Meat Sci. 2004;68:285–90.

  14. 14.

    Filteau M, Pavey SA, St-Cyr J, Bernatchez L. Gene Coexpression networks reveal key drivers of phenotypic divergence in Lake whitefish. Mol Biol Evol. 2013;30(6):1384–96.

  15. 15.

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

  16. 16.

    Kogelman LJA, Cirera S, Zhernakova DV, Fredholm M, Franke L, Kadarmideen HN. Identification of co-expression gene networks, regulatory genes and pathways for obesity based on adipose tissue RNA sequencing in a porcine model. BMC Med Genet. 2014;7(1):57.

  17. 17.

    Langfelder P, Mischel PS, Horvath S. When is hub gene selection better than standard meta-analysis? PLoS One. 2013;8:e61505.

  18. 18.

    Bligh EG, Dyer WJ. A rapid method of total lipid extraction and purification. Can J Biochem Phys. 1959;37(8):911–7.

  19. 19.

    Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics. 2014;30:2114–20.

  20. 20.

    Kim D, Langmead B, Salzberg SL. HISAT: a fast spliced aligner with low memory requirements. Nat Methods. 2015;12:357–60.

  21. 21.

    Trapnell C, Roberts A, Goff L, Pertea G, Kim D, Kelley DR, et al. Differential gene and transcript expression analysis of RNA-Seq experiments with TopHat and cufflinks. Nat Protocols. 2012;(3):562–78.

  22. 22.

    Trapnell C, Hendrickson DG, Sauvageau M, Goff L, Rinn JL, Pachter L. Differential analysis of gene regulation at transcript resolution with RNA-seq. Nat Biotechnol. 2013;31:46–53.

  23. 23.

    Benjamini Y, Hocheberg Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J R Stat Soc. 1995;57(1):1–12.

  24. 24.

    Huang W, Sherman BT, Lempick RA. Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat Protocols. 2009;4:44–57.

  25. 25.

    Kanehisa M, Goto SKEGG. Kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 2000;28(1):27–30.

  26. 26.

    Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, et al. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 2003;13(11):2498–504.

  27. 27.

    Hui S, Sander C, Potylitsine E, Whitaker W, Bader G, Morrison L, et al. ExpressionCorrelation. 2015. Available online at:

  28. 28.

    Cline MS, Smoot M, Cerami E, Kuchinsky A, Landys N, Workman C, et al. Integration of biological networks and gene expression data using Cytoscape. Nat Protoc. 2007;2:2366–82.

  29. 29.

    Chin CH, Chen SH, Wu HH, Ho CW, Ko MT, Lin CY. cytoHubba: identifying hub objects and sub-networks from complex interactome. BMC Syst Biol. 2014;8(4):11.

  30. 30.

    Stevens A, De Leonibus C, Hanson D, Dowsey AW, Whatmore A, Meyer S, et al. Network analysis: a new approach to study endocrine disorders. J Mol Endocrinol. 2013;52:R79–93.

  31. 31.

    Newman ME. Modularity and community structure in networks. Proc Natl Acad Sci U S A. 2006;103:8577–82.

  32. 32.

    Goff L, Trapnell C, Kelley D. cummeRbund: Analysis, exploration, manipulation, and visualization of cufflinks high-throughput sequencing data. R package version 2.24.0. 2018.

  33. 33.

    Qi Y, Kapterian TS, Du X, Ma Q, Fei W, Zhang Y, et al. CDP-diacylglycerol synthases regulate the growth of lipid droplets and adipocyte development. J Lipid Res. 2016;57(5):767–80.

  34. 34.

    Akisu M, Kultursay N, Coker I, Huseyinov A. Myocardial and hepatic free carnitine concentrations in pups of diabetic female rats. Ann Nutr Metab. 2002;46:45–8.

  35. 35.

    Santiago JL, Martínez A, de la Calle H, Fernández-Arquero M, Figueredo MÁ, de la Concha EG, et al. Evidence for the association of the SLC22A4 and SLC22A5 genes with type 1 diabetes: a case control study. BMC Medical Genetics. 2006;7:54.

  36. 36.

    Vannucci SJ, Maher F, Simpson IA. Glucose transporter proteins in brain: delivery of glucose to neurons and glia. Glia. 1997;21:2–21.

  37. 37.

    Zhao FQ, Glimm DR, Kennelly JJ. Distribution of mammalian facilitative glucose transporter messenger RNA in bovine tissue. Int J BioChemiPhysics. 1993;25(12):1897–903.

  38. 38.

    Shimamoto S, Ijiri D, Kawaguchi M, Nakashima K, Ohtsuka A. Gene expression pattern of glucose transporters in the skeletal muscles of newly hatched chicks. Biosci Biotechnol Biochem. 2016;80(7):1382–5.

  39. 39.

    Brown MS, Goldstein JL. Selective versus total insulin resistance: a pathogenic paradox. Cell Metab. 2008;7(2):95–6.

  40. 40.

    Matsumoto M, Pocai A, Rossetti L, Depinho RA, Accili D. Impaired regulation of hepatic glucose production in mice lacking the forkhead transcription factor Foxo1 n liver. Cell Metab. 2007;6:208–16.

  41. 41.

    Munekata K, Sakamoto K. Forkhead transcription factor Foxo1 is essential for adipocyte differentiation. In Vitro Cell Dev Biol Anim. 2009;45:642–51.

  42. 42.

    Sparks JD, Dong HH. FoxO1 and hepatic lipid metabolism. Curr Opin Lipidol. 2009;20(3):217–26.

  43. 43.

    Li Y, Ma Z, Jiang S, Hu W, Li T, Di S, et al. A global perspective on FOXO1 in lipid metabolism and lipid-related diseases. Prog Lipid Res. 2017;66:42–9.

  44. 44.

    Rana R, Surapureddi S, Kam W, Ferguson S, Goldstein JA. Med25 is required for RNA polymerase ii recruitment to specific promoters, thus regulating xenobiotic and lipid metabolism in human liver. Mol Cell Biol. 2011;31(3):466–81.

  45. 45.

    Mizoguchi Y, Hirano T, Itoh T, Aso H, Takasuga A, Sugimoto Y, et al. Differentially expressed genes during bovine intramuscular adipocyte differentiation profiled by serial analysis of gene expression. Anim Genet. 2010;41(4):436–41.

  46. 46.

    Zhao C, Tian F, Yu Y, Luo J, Mitra A, Zhan F, et al. Functional genomic analysis of variation on beef tenderness induced by acute stress in angus cattle. Comparative and Functional Genomics. 2012;2012:756284.

  47. 47.

    Lee SH, Gondro C, Van der Werf J, Kim NK, Lim DJ, Park EW, et al. Use of a bovine genome array to identify new biological pathways for beef marbling in Hanwoo (Korean cattle). BMC Genomics. 2010;11:623.

  48. 48.

    Melo TP, Camargo GMF, Albuquerque LG, Carvalheiro R. Genome-wide association study provides strong evidence of genes affecting the reproductive performance of Nellore beef cows. PLoS One. 2017;12(5):e0178551.

  49. 49.

    Neyen CD, Gordon S. Macrófagos em lipid and immune homeostasis. eLS. John Wiley & Sons Ltd. 2014;1, 1.

  50. 50.

    Ravnskjaer K, Madiraju A, Montminy M. Role of the cAMP pathway in glucose and lipid metabolism. Handb Exp Pharmacol. 2016;233:29–49.

  51. 51.

    Dimitriadis G, Mitrou P, Lambadiari V, Maratou E, Raptis SA. Insulin effects in muscle and adipose tissue. Diabetes Res Clin Pract. 2011;93(Suppl 1):S52–9.

  52. 52.

    Ikeoka D, Krusinova E. Insulin resistance and lipid metabolism. Rev Assoc Med Bras. 2009;55(3):1.

  53. 53.

    Lindholm-Perry AK, Keel BN, Zarek CM, Keele JW, Kuehn LA, Snelling WM, et al. 325 genes in skeletal muscle associated with gain and intake identified in a multiseason study of crossbred beef steers. J Ani Sci. 2017;95(4):161.

  54. 54.

    Casci T. Network fundamentals, via hub genes. Nat Rev Genet. 2006;7:664–5.

  55. 55.

    Kim JH, Ovilo C, Park EW, Fernndez A, Lee JH, Jeon JT, et al. Minimizing a QTL region for intramuscular fat content by characterizing the porcine phosphodiesterase 4B (PDE4B) gene. BMB Rep. 2008;41(6):466–71.

  56. 56.

    Lee KT, Byun MJ, Kang KS, Park EW, Lee SH, Cho S. Et al. neuronal genes for subcutaneous fat thickness in human and pig are identified by local genomic sequencing and combined SNP association study. PLoS One. 2011;6(2):e16356.

  57. 57.

    Schafer PH, Parton A, Capone L, Cedzik D, Brady H, Evans JF, et al. Apremilast is a selective PDE4 inhibitor with regulatory effects on innate immunity. Cell Signal. 2016;26(9):2016–29.

  58. 58.

    Russell TR, Ho R. Conversion of 3T3 fibroblasts into adipose cells: triggering of differentiation by prostaglandin F2alpha and 1-methyl-3-isobutyl xanthine. Proc Natl Acad Sci U S A. 1976;73(12):4516–20.

  59. 59.

    Houslay MD, Sullivan M, Bolger GB. The multienzyme PDE4 cyclic adenosine monophosphate-specific phosphodiesterase family: intracellular targeting, regulation, and selective inhibition by compounds exerting anti-inflammatory and antidepressant actions. Adv Pharmacol. 1998;44:225–342.

  60. 60.

    Conti M, Jin SL. The molecular biology of cyclic nucleotide Phosphodiesterases. Prog Nucleic Acids Res Mol Biol. 1999;63:1–38.

  61. 61.

    Liu H, Palmer D, Jimmo SL, Tilley DG, Dunkerley HA, Pang SC, et al. Expression of phosphodiesterase 4D (PDE4D) is regulated by both the cyclic AMP-dependent protein kinase and mitogen-activated protein kinase signaling pathways. A potential mechanism allowing for the coordinated regulation of PDE4D activity and expression in cells. J Biol Chem. 2000;275(34):26615–24.

  62. 62.

    Xu Z, Huo J, Ding X, Yang M, Li L, Dai J, Hosoe K, et al. Coenzyme Q10 improves lipid metabolism and ameliorates obesity by regulating CaMKII-mediated PDE4 inhibition. Sci Rep. 2017;7(1):8253.

  63. 63.

    Petersen RK, Madsen L, Pedersen LM, et al. Cyclic AMP (cAMP)-mediated stimulation of adipocyte differentiation requires the synergistic action of Epac- and cAMP-dependent protein kinase-dependent processes. Mol Cell Biol. 2008;28(11):3804–16.

  64. 64.

    Sudrajad P, Sharma A, Dang CG, Kim JJ, Kim KS, Lee JH, et al. Validation of single nucleotide polymorphisms associated with carcass traits in a commercial Hanwoo population. Asian-Australas J Anim Sci. 2016;29(11):1541–6.

  65. 65.

    Dhanoa BS, Cogliati T, Satish AG, Bruford EA, Friedman JS. Update on the Kelch-like (KLHL) gene family. Hum Genomics. 2013;7:13.

  66. 66.

    Gupta VA, Beggs AH. Kelch proteins: emerging roles in skeletal muscle development and diseases. Skelet Muscle. 2014;4:11.

  67. 67.

    Attaix D, Ventadour S, Codran A, Béchet D, Taillandier D, Combaret L. The ubiquitin-proteasome system and skeletal muscle wasting. Essays Biochem. 2005;41:173–86.

  68. 68.

    Piórkowska K, Żukowski K, Nowak J, Połtowicz K, Ropka-Molik K, Gurgul A. Genome-wide RNA-Seq analysis of breast muscles of two broiler chicken groups differing in shear force. Anim Genet. 2016;47(1):68–80.

  69. 69.

    Cullinan EB, Kwee L, Nunes P, Shuster DJ, Ju G, McIntyre KW, Chizzonite RA, Labow MA. IL-1 receptor accessory protein is an essential component of the IL-1 receptor. J Immunol. 1998;161(10):5614–20.

  70. 70.

    Tsao CH, Shiau MY, Chuang PH, Chang YH, Hwang J. Interleukin-4 regulates lipid metabolism by inhibiting adipogenesis and promoting lipolysis. J Lipid Res. 2014;55(3):385–97.

  71. 71.

    Tocci MJ, Schmidt JA. In: Remick DG, Friedland JS, editors. Interleukin-1: structure and function. Cytokines in health and disease, second edition. New York: Marcel Dekker, Inc; 1997. p. 1–27.

  72. 72.

    Getz GS, Reardon CA. The mutual interplay of lipid metabolism and the cells of the immune system in relation to atherosclerosis. Clin Lipidol. 2014;9(6):657–71.

  73. 73.

    Dowds CM, Kornell SC, Blumberg RS, Zeissig S. Lipid antigens in immunity. Biol Chem. 2014;395(1):61–81.

  74. 74.

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

Download references


The authors thank the Qualitas Nelore breeding program company for providing the tissue samples and database used in this study.


This study was financed in part by the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior – Brazil (CAPES) – Finance Code 001 and São Paulo Research Foundation - FAPESP (#grant 2009/16118-5 and #grant 2015/16850-9).

Author information

DBSS, LFSF, FB, JAF, LALC and LGA developed and planned the design of the study. DBSS, LFSF, MMMM, AFBM and DGP led performance of the experiments and interpreted the results. DBSS, LFSF, AFBM, DGP and LGA drafted and revised the manuscript. All authors read and approved the final manuscript.

Correspondence to Danielly Beraldo dos Santos Silva.

Ethics declarations

Ethics approval and consent to participate

All experimental procedures were approved by Ethics Committee of the School of Agricultural and Veterinarian Sciences, São Paulo State University (UNESP), Jaboticabal, SP, Brazil (protocol number 18.340/16). The animals were provided by Qualitas Nelore breeding program company and they were slaughtered in commercial slaughterhouses. These slaughterhouses have animal welfare departments staffed by professionals trained by WAG (World Animal Protection) to ensure that the animals are killed humanely using a captive bolt pistol for stunning.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Additional files

Additional file 1:

Table A1. Biological Process GO terms and pathways obtained with the DAVID software for differentially expressed genes in Nelore catlle muscle divergent for intramuscular fat content phenotypes. Table A2. Parameters from CytoHubba Plugin from Cytoscape software and hub genes. Table A3. Intramuscular fat associated modules as from co-expression network analysis. Table A4. Biological Process GO terms and pathways obtained with the Reactome FI Plugin from Cytoscape software for co-expressed genes network modules. Table A5. Intramuscular fat associated modules as from co-expression subnetwork analysis. (XLSX 75 kb)

Additional file 2:

Figure S1. Global statistics and quality control. a Box plot of FPKM distributions for individual conditions. b PCA plot for gene-level features. (PDF 65 kb)

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (, which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark


  • Bovine
  • Co-expression
  • Genes
  • Lipids
  • Transcriptome