Skip to main content
Download PDF
  • Research article
  • Open access
  • Published:

Gene co-expression networks associated with carcass traits reveal new pathways for muscle and fat deposition in Nelore cattle

BMC Genomics volume 20, Article number: 32 (2019) Cite this article

Abstract

Background

Positively correlated with carcass weight and animal growth, the ribeye area (REA) and the backfat thickness (BFT) are economic important carcass traits, which impact directly on producer’s payment. The selection of these traits has not been satisfactory since they are expressed later in the animal’s life and multigene regulated. So, next-generation technologies have been applied in this area to improve animal’s selection and better understand the molecular mechanisms involved in the development of these traits. Correlation network analysis, performed by tools like WGCNA (Weighted Correlation Network Analysis), has been used to explore gene-gene interactions and gene-phenotype correlations. Thus, this study aimed to identify putative candidate genes and metabolic pathways that regulate REA and BFT by constructing a gene co-expression network using WGCNA and RNA sequencing data, to better understand genetic and molecular variations behind these complex traits in Nelore cattle.

Results

The gene co-expression network analysis, using WGCNA, were built using RNA-sequencing data normalized by transcript per million (TPM) from 43 Nelore steers. Forty-six gene clusters were constructed, between them, three were positively correlated (p-value< 0.1) to the BFT (Green Yellow, Ivory, and Light Yellow modules) and, one cluster was negatively correlated (p-value< 0.1) with REA (Salmon module). The enrichment analysis performed by DAVID and WebGestalt (FDR 5%) identified eight Gene Ontology (GO) terms and three KEGG pathways in the Green Yellow module, mostly associated with immune response and inflammatory mechanisms. The enrichment of the Salmon module demonstrated 19 GO terms and 21 KEGG pathways, related to muscle energy metabolism, lipid metabolism, muscle degradation, and oxidative stress diseases. The Ivory and Light yellow modules have not shown significant results in the enrichment analysis.

Conclusion

With this study, we verified that inflammation and immune response pathways modulate the BFT trait. Energy and lipid metabolism pathways, highlighting fatty acid metabolism, were the central pathways associated with REA. Some genes, as RSAD2, EIF2AK2, ACAT1, and ACSL1 were considered as putative candidate related to these traits. Altogether these results allow us to a better comprehension of the molecular mechanisms that lead to muscle and fat deposition in bovine.

Background

Heavy carcass weight is critical for meat producers since it is used as the primary parameter for payment in the slaughterhouses [1]. Positively correlated with carcass weight and animal growth, the ribeye area (REA) can be used as an indicator of muscularity, prime cuts, and the edible mass of carcass. Inversely proportional to REA, the backfat thickness (BFT) is related to the percentage of fat in the carcass. BFT is essential to protect the carcass during cooling, avoiding problems such as cold shortening, drip loss and dark cutting [2,3,4].

Despite the economic relevance for producers and to affect the final weight of the animals, the selection of these traits has not been satisfactory, since their maximum potential is expressed later in animal’s life [1]. Therefore, a better understanding of the biological processes that regulate these characteristics could help to elucidate the mechanisms of genetic inheritance, and consequently will increase animal selection accuracy. In this context, next-generation sequencing technologies have revolutionized genome and transcriptome analysis of complex organisms, generating large datasets and allowing the identification of new genes, metabolic pathways and biological processes that influence the phenotype [5,6,7]. Moreover, novel approaches have been developed to more reliably analyze large and multivariate datasets by associating them with traits of interest. Correlation network analysis has been widely used to analyze large datasets as a promissory method from systems biology, able to represent the complexity of a cellular transcription network [8,9,10].

The Weighted Correlation Network Analysis (WGCNA) is a system biology method that is used to explore the correlation patterns among genes in transcriptomic studies, providing a unique insight into the structure and behavior of molecular interactions [8, 9]. The WGCNA describes and permits the visualization of networks derived from large datasets. WGCNA can be used to explore the structure of modules within a co-expression network, to measure the relationship between genes and modules (module membership), explore the relationship between modules, or even to rank genes or modules associated to the studied traits [8].

Several works utilized gene co-expression networks analysis, in particular, the WGCNA tool, to evaluate complex traits in different species such as mice [11], humans [12], pork [13, 14], lamb [10] and, cattle [15,16,17,18] demonstrating the correlation between genes and phenotype successfully.

Kong et al. [16] utilizing the WGCNA tool to understand the molecular differences between efficient and inefficient cattle concerning residual feed intake in a Hereford x Angus population, identified a significant module with 764 genes negatively correlated with the trait of interest. With the results, the authors could infer that efficient animals probably have an increased energy production and better absorption of food nutrients compared to inefficient ones. Sabino et al. [10] used the WGCNA to do a nutrigenomics investigation in lambs, considering diet and sex differences in muscle and liver tissue, revealing a sex-dependent dietary effect on the transcriptome of the studied animals. In previous work from our group, Oliveira et al. [18] employed WGCNA tool to do an integrative miRNA-mRNA (microRNA – messenger RNA) study associated with intramuscular fat deposition in Nelore cattle, revealing potential regulatory mechanisms of gene signaling networks involved in fat deposition in bovine.

The present study aimed to gain molecular insights into economic important carcass traits of Nelore cattle and to identify putative candidate genes and metabolic pathways that regulate REA and BFT, by constructing a gene co-expression network using WGCNA and RNA sequencing data. Modules correlated with REA and BFT were identified, and the genes within each significant module were extracted to perform functional enrichment analysis, permitting us to better understand gene interaction, biological processes, and the metabolic pathways behind these complex traits.

Results

For this study, we used RNA-sequencing normalized data by transcript per million (TPM) from 14,529 genes of 43 animals with contrasting genomic estimated breeding values (GEBV). Table 1 shows animal identification, phenotypic values, GEBVs, number of raw reads, mapped reads and percentage of mapped reads. The heritability values for REA and BFT were h2 = 0.27 and h2 = 0.21, respectively [19]. Additional file 1: Table S1 shows the contrasting GEBV groups for both traits, demonstrating that mean values in the low group were statistically different (p-value< 0.001) from the high group for REA and BFT. For REA, mean GEBV were − 2.94 in the low group (standard deviation, SD = 0.59) and, 3.14 in the high group (SD = 0.56); for BFT, mean GEBV were − 0.91 (SD = 0.12) and 1.24 (SD = 0.28) in the low and high group, respectively. The correlation analysis between the phenotypic (REA mean and SD values = 59.75 ± 9.50 cm2; BFT mean and SD values = 7.00 ± 3.41 mm) and GEBV values both for REA and BFT showed a strong correlation, with r = 0.79 and r = 0.84 for REA and BFT respectively, justifying our selection using genomic values. On the other hand, the correlation analysis between the GEBV of REA x BFT demonstrated that these traits were independent of one another in this particular study, with r = − 0.06.

Table 1 Phenotypic and genomic values, number of raw-reads, number and percentage of mapped reads from the selected 43 Nelore cattle
Full size table

The gene co-expression network analysis, performed using WGCNA program, resulted in the identification of 46 module eigengenes (ME) (Fig. 1). Figure 2 illustrates the hierarchical clustering tree (dendrogram) of all genes and modules colors and, in the supplementary material, there is a heatmap plot of the gene network (Additional file 2: Figure S1). Among the 46 identified ME, the Ivory, Light Yellow, Green Yellow, and Salmon modules showed a significant module-trait association (p-value< 0.1) with at least one of the studied phenotypes, demonstrating positive correlations of r = 0.3 with BFT (Ivory, Light Yellow and Green Yellow) and, a negative correlation of r = − 0.3 with REA (Salmon) (Fig. 1). The entire lists of genes in each of these modules were further analyzed using DAVID (Database for Annotation, Visualization and Integrated Discovery) version 6.8 and, WebGestalt 2017 (WEB-based GEne SeT AnaLysis Toolkit).

Fig. 1
figure 1

Module-trait associations between the module eigengenes (ME) and the studied traits, ribeye area (REA) and backfat thickness (BFT). Each row corresponds to a module eigengene, column to a trait. Each cell contains the Pearson’s correlation coefficients (numbers outside parentheses) and, the p-values of the correlation (numbers within parentheses). The graphic is color-coded by correlation according to the color legend, red represents a positive correlation and blue represents a negative one

Full size image
Fig. 2
figure 2

Cluster dendrogram of all genes from the selected Nelore steers. Cluster dendrogram of all genes, with dissimilarity based on topological overlap. The different colors in the bottom represent gene modules

Full size image

The Green Yellow module, positively correlated to BFT (p-value< 0.1), presented 146 co-expressed genes assigned for the enrichment analysis. Eight Gene Ontology (GO) terms divided into four Biological Processes (BP), four Molecular Functions (MF), and two KEGG (Kyoto Encyclopedia of Genes and Genomes) pathways were identified in the analysis performed using DAVID (FDR 5%) (Fig. 3, Additional file 3: Table S2).

Fig. 3
figure 3

Functional enrichment analysis from the gene list of the Green Yellow module, performed by DAVID v6.8 (FDR < 0.05). Gene Ontology terms: BP – Biological Process; MF – Molecular Function

Full size image

The two KEGG pathways found were Influenza A (bta05164) and Herpes simplex infection (bta05168). The biological processes were defense response to virus (GO:0051607), innate immune response (GO:0045087), negative regulation of viral genome replication (GO:0045071) and, ISG15-protein conjugation (GO:0032020). The MF were double-stranded RNA binding (GO:0003725), ubiquitin-protein transferase activity (GO:0004842), NAD+ ADP-ribosyltransferase activity (GO:0003950) and, GTPase activity (GO:0003924). Most of them related to inflammation mechanisms and the immune system.

The functional enrichment analysis by WebGestalt revealed three KEGG pathways (Table 2), the different pathway identified in this analysis was the NOD-like receptor signaling pathway (bta04621), also related to immune response.

Table 2 KEGG Pathways (FDR < 0.05) identified by WebGestalt 2017 from the gene list of the Green Yellow module
Full size table

The Salmon module, negatively correlated with REA (p-value< 0.1), was constituted by 136 genes. Figure 4 and Additional file 4: Table S3 demonstrates the results of the functional enrichment analysis performed by DAVID (FDR 5%), we found five BP, eleven Cellular Components (CC) and, three MF terms.

Fig. 4
figure 4

Functional enrichment analysis from the gene list of the Salmon module, performed by DAVID v6.8 (FDR < 0.05). Gene Ontology terms: BP – Biological Process; CC – Cellular Component; MF – Molecular Function

Full size image

The biological processes were mostly related to energy metabolism, like ATP synthesis coupled proton transport (GO:0015986), tricarboxylic acid cycle (GO:0006099) and, fatty acid beta-oxidation (GO:0006635). From the eleven CC identified, seven were mitochondrial constituents, such as mitochondrial inner membrane (GO:0005743), mitochondrial respiratory chain complex I (GO:0005747), mitochondrial matrix (GO:0005759) and, mitochondrial proton-transporting ATP synthase complex (GO:0005753). The MF were proton-transporting ATP synthase activity (GO:0046933), NADH dehydrogenase (ubiquinone) activity (GO:0008137) and, electron carrier activity (GO:0009055), also associated to the muscle energy metabolism.

Seventeen KEGG pathways were identified by DAVID (Fig. 4, Additional file 4: Table S3) contrasting with 21 showed in WebGestalt enrichment analysis (FDR 5%, Additional file 5: Table S4). All 17 pathways found by DAVID also appeared in WebGestalt analysis. The four different KEGG pathways identified by the second program (Table 3) were Arginine biosynthesis (bta00220), Phenylalanine, tyrosine and tryptophan biosynthesis (bta00400), Butanoate metabolism (bta00650) and, Lysine degradation (bta00310). Among the common pathways, we can highlight Fatty acid metabolism (bta01212), Oxidative phosphorylation (bta00190) and Citrate cycle (TCA cycle) (bta00020) associated to the muscle energy metabolism; and PPAR signaling pathway (bta03320), associated to the muscle lipid metabolism. And some pathways related to muscle degradation and oxidative stress diseases, such as Parkinson’s disease (bta05012) and Alzheimer’s disease (bta05010).

Table 3 KEGG Pathways (FDR < 0.05) unique identified by WebGestalt 2017 from the gene list of the Salmon module
Full size table

At last, the Ivory and Light Yellow modules, positively correlated to the BFT (p-value< 0.1), presented 17 co-expressed genes and, 87 co-expressed genes respectively (Additional files 6 and 7: Tables S5 and S6, respectively). The genes within these modules were not significantly enriched by DAVID (FDR 5%) nor by WebGestalt (FDR 5%).

Discussion

The commercial value of the bovine carcass is determined by the adequate development of muscle and adipose tissue. Therefore, an increase in muscle and fat masses in growing cattle is an important issue for farmers and beef industry. Studies involving Bos indicus animals have shown a lower propensity for subcutaneous fat deposition and Longissimus muscle area when compared to Bos taurus breeds [2, 3]. Compared to other Nelore studies, our phenotypic values for REA (mean value = 59.75 cm2) can be considered on the average, taking into account that studies in Nelore cattle show amounts varying from 40 to 100 cm2 [1, 2, 20, 21]. On the other hand, the BFT phenotypic values (mean value = 7.00 mm) were higher than that found in the literature for Nelore, with mean values ranging from 1.93 to 4.84 mm and, even higher than presented in some Bos indicus x Bos taurus crossbreed studies [1, 3, 20,21,22].

The variation in ribeye area and backfat thickness in Nelore cattle is harmful both for producers and the beef industry since these traits influence carcass yield and consequently, producer’s payment [3]. It’s well known that both environment and genetics contribute to the phenotypic variations within a population or in the same breed. So, to better understand the genetic and molecular influences behind these economically important traits, our goal was to construct co-expressed gene networks and identify putative candidate genes and metabolic pathways that regulate these traits, using the WGCNA tool and RNA sequencing data.

Forty-six gene clusters were constructed, between them, three were positively correlated (p-value< 0.1) to the BFT (Green Yellow, Ivory, and Light Yellow modules) and, one cluster was negatively correlated (p-value< 0.1) with REA (Salmon module). From the three positive correlated modules assigned to BFT, just the genes within the Green Yellow were significantly (FDR 5%) enriched by DAVID and WebGestalt.

The functional enrichment analysis for the co-expressed genes from the Green Yellow module demonstrated that they are mainly involved in inflammatory mechanisms and immune response. In previous work, Oliveira et al. [18], studying intramuscular fat deposition in this same Nelore cattle population also found co-expressed gene modules enriched for inflammatory response and immune system GO terms. Tao et al. [23] studying the muscle and adipose transcriptomic profile of Indigenous x Western Chinese pig breeds associated with growth performance and quality carcass traits – intramuscular fat, marbling, and loin muscle area –, identified genes related to immune response and inflammation. These authors identified highly expressed and up-regulated genes in the muscle tissue of Chinese Indigenous breeds (higher BFT content) enriched for immune response, mitochondria, Herpes simplex infection, Parkinson’s disease, and apoptosis GO terms, corroborating our findings.

The backfat, not only, is important in the industrialization process acting as a thermal insulator during carcass cooling but is an indispensable source of energy in the animal’s body, carrying fat-soluble vitamins and essential fatty acids [2, 3, 24, 25]. This trait also is representative of the total body fat in the carcass. According to Schröder and Staufenbiel [26], an increase of just 1 mm in the BFT reflects in approximately 5 kg of the total body fat content in bovine.

The increase in total body fat is defined by an adipose tissue expansion together with adipocyte hypertrophy [27]. The fat deposition is a consequence of the balance between energy intake and expenditure. The excess of nutrients and energy can lead to an increase in fat accumulation – also called obesity in human and mouse studies – and, consequentially activate inflammatory and stress responses. So, this inflammatory process created by the increase in fat can disrupt systemic metabolic homeostasis and inhibit insulin receptor signaling, involving immune cells and immune response pathways. Although the adipose tissue is the primary source of inflammation induced by excess fat accumulation, it is sufficient to activate inflammatory signaling pathways and increase the amount of pro-inflammatory immune cells in other tissues, like skeletal muscle, liver, pancreas, and brain [28,29,30,31,32]. The molecular mechanisms behind fat deposition in bovine are still unclear, so studying the lipid metabolism from another mammalian species can clarify our knowledge about the differences in backfat in this cattle population.

In our preview study, with the same Nelore population [64], we identified biological processes related to the immune response in differentially expressed genes assigned to the BFT trait, similarly to the results found herein. Interestingly, the RSAD2 gene (Radical domain of S-adenosyl methionine containing 2) was down-regulated in the group with lower GEBV values for BFT in our previous study and, herein it was identified in Green Yellow module, which was positively correlated to BFT (p-value< 0.1).

RSAD2 is an interferon-regulated gene associated with innate immune response during viral infections [33, 34]. Additionally, this gene participates in lipid biosynthesis and the modulation of lipid droplet contents [35, 36]. Dogan et al. [36], working with obese-induced mouse detected higher expression of RSAD2 in animals with lower fat amounts, showing that this gene controls lipid droplets formation, but it is the RSAD2 impairment that drives fat accumulation. These authors further associated this gene with endoplasmic reticulum (ER) stress and the activation of inflammatory mechanisms in obese animals. According to Warfel et al. [31] and Ramos-Lopez et al. [37], one of the main contributors to the activation of inflammatory pathways and immune response during obesity is the metabolic stress of organelles, such as ER and mitochondria.

Like RSAD2, the EIF2AK2 (Eukaryotic translation initiation factor 2 alpha kinase 2; Alias PRKR, PKR) is a key inducer of inflammation, responsive to interferons (IFN) and, associated with ER stress and fat accumulation [37, 38]. In Green Yellow module, EIF2AK2 was associated with inflammatory pathways, and immune response GO terms. This gene encodes a protein called PKR, a serine/threonine kinase, activated by double-stranded RNA, cytokines, stress signals, and IFN. The PKR protein can be activated by lipids and, participate in major inflammatory signaling events, such as JNK (c-Jun N-terminal kinase) activation during lipids exposition and ER stress [28, 38]. Nakamura et al. [28] hypothesized that the increase in PKR activity in obese states could be caused by an excess of energy and nutrient supply, representing an adaptive attempt to interact with synthetic pathways that would further accumulate energy. In our study, higher expression of the EIF2AK2 gene was positively correlated with BFT, supporting the association of this gene with fat accumulation.

Another gene family identified in the Green Yellow module was poly (ADP-ribose) polymerases (PARP) represented by PARP9, PARP10, PARP12, and PARP14. PARPs activity is stimulated by excess fat, high-fat diet, aging, oxidative stress, DNA damage and, inflammation states. These enzymes are involved in lipid metabolism, mostly by controlling redox balance and NAD+ homeostasis in mitochondrial metabolism [39, 40]. According to Jokinen et al. [39], obesity can be characterized by low levels of NAD+ in the adipose tissue that can stimulate PARP activity. Mohamed et al. [41] examined the effect of high-fat diet in mouse skeletal muscle cultured cells (C2C12), and they found that oversupplied obese animals had their levels of PARP2 increased together with reduced mitochondrial functions and NAD+ levels in the cultured muscle cells.

The NOD-like receptor signaling pathway (bta04621), the only one enriched by WebGestalt (Table 2), participates in the innate immune response and, can be activated by increased levels of glucose and free fatty acids, as well as reactive oxygen species derived from the inflammation state during obesity [27, 42]. Yin et al. [27] studying adipocytes metabolism from the adipose tissue of obese x non-obese women, verified an up-regulation of the NOD-like receptor signaling pathway in the adipocytes from the obese group. So, genes and metabolic pathways presented in the Green Yellow module, related to immune response and inflammatory processes, have an important paper in the lipid metabolism in mammals, demonstrating that there is a correlation between the increase in the expression of genes involved in these pathways and the BFT in cattle.

The other significantly enriched gene module, Salmon, showed a negative correlation (r = − 0.3, p-value< 0.1) with the ribeye area. REA is an important quality carcass trait related to the amount of meat, that is used as an indicator of cuts yield, the percentage of muscle, animal growth and carcass weight [1,2,3, 43].

Some of some GO terms and KEGG pathways identified here, like the fatty acid beta-oxidation (GO:0006635), mitochondrion (GO:0005739), mitochondrial proton-transporting ATP synthase complex (GO:0005753), Citrate cycle (TCA cycle) (bta00020) and, Oxidative phosphorylation (bta00190) are part of the complex cascade of events that occur in the skeletal muscle for generate energy. Muscles have an essential role in energy metabolism, the regulation of skeletal muscle metabolism involves multiple pathways and different molecules committed in the uptake and storage of energy. The glucose and fatty acid metabolism are the major sources of energy in this tissue. Their energy demand is mainly fulfilled by phosphocreatine and ATP produced during glucose and fatty acid oxidation [44]. In the glucose metabolism, glucose delivered by blood enters the myocyte, and its oxidation generates energy by phosphorylation. In the fatty acid metabolism, the primary source of energy for the muscle is the non-esterified fatty acids (NEFA) derived from circulation and, from lipolysis of triacylglycerols (TG) located mostly in the adipose tissue or, accumulated in the muscle (intramuscular TG). Once in the cytosol, NEFA are esterified to long-chain acyl CoA that is destined for mitochondrial beta-oxidation and, subsequently enter the TCA cycle generating ATP [45, 46]. Cesar et al. [47] and Oliveira et al. [18], studying fat-related traits in this Nelore population, also identified energy metabolism pathways in the skeletal muscle of the animals.

Marrades et al. [48], investigating two groups of subjects (lean versus obese) high feed diet consuming, identified ACADM gene (Medium chain acyl-CoA dehydrogenase) in the Fatty acid β-oxidation pathway and, SUCLG2 (β subunit succinate-CoA ligase, GDP-forming) in the TCA cycle pathway, down-regulated in the subcutaneous adipose tissue of obese subjects. Likewise, Jeong et al. [49] identified the TCA Cycle and Fatty Acid Oxidation genes in the Longissimus muscle of Korean bulls, following castration. Their enrichment analysis also demonstrated CPT1B (Carnitine palmitoyltransferase IB) and Hydroxyacyl-CoA dehydrogenase (HADH) genes assigned to Fatty oxidation pathway; and, Succinate-CoA ligase [GDP-forming] (SUCLG), and Isocitrate dehydrogenase (IDH) gene families in the TCA Cycle pathway. Also, they found some ATP (Adenosine triphosphatase), NDUF (NADH dehydrogenase [ubiquinone]) and COX (Cytochrome c oxidase) genes present in the Oxidative phosphorylation pathway, similar to our findings.

Three HADH genes appeared in Salmon module enrichment analysis, HADH (Hydroxyacyl-CoA dehydrogenase), HADHA (Hydroxyacyl-CoA Dehydrogenase Trifunctional Multienzyme Complex Subunit Alpha) and, HADHB (Hydroxyacyl-CoA Dehydrogenase Trifunctional Multienzyme Complex Subunit Beta). Together they were identified in ten KEGG pathways, including two of the WebGestalt analysis (Table 3 and Additional file 4: Table S3). Costa et al. [50], studying the fatty acid profile of the Longissimus muscle of a cattle population, found HADHA gene associated with several terms related to the fatty acid metabolism, like oxidation of lipid and, accumulation of specific fatty acids. Furthermore, HADH enzymes are crucial in the mitochondrial beta-oxidation, participating in the two final steps of this process [51]. According to Zhang et al. [45] and Xu et al. [51], an up-regulation of the mitochondrial beta-oxidation process leads to a lower body fat content due to diminishing TG content. Maybe this could explain why some fatty acid metabolism terms and pathways appeared associated with REA, in our skeletal muscle study.

Adjacent to the pathways and terms associated with energy metabolism, we also found some associated with the muscle lipid metabolism, like the PPAR signaling pathway (bta03320). The PPAR pathway is responsible for the regulation of adipocyte tissue development, adipogenic differentiation and, lipogenesis. The PPARs (peroxisome proliferator-activated receptors) are nuclear receptors that take part in a number of biological processes, like skeletal muscle lipid oxidation, inflammation, mitochondrial respiration, energy homeostasis and, thermogenesis [52,53,54]. Several studies identified PPAR genes associated with fat traits in cattle [18, 47, 53,54,55].

Huang et al. [54] analyzing the transcriptomic profile of the subcutaneous adipose tissue from Wagyu and Holstein cattle confirmed the importance of the PPAR signaling pathway as a key regulator of lipid metabolism in bovine. These authors also found acetyl-CoA acyltransferase 1 (ACAT1) and acyl-CoA synthetase long-chain family member 1 (ACSL1) genes up-regulated in the BFT of Wagyu cattle. Here, the ACAT1 and ACSL1 genes were found in a gene module negatively associated with REA. These genes not only were verified in the PPAR signaling pathway; ACAT1 was identified in four GO terms and nine pathways (Additional file 4: Table S3), two of them only by WebGestalt (Table 3) – Butanoate metabolism (bta00650) and Lysine degradation (bta00310). The ACSL1 was found in the mitochondrion (GO:0005739) and, four KEGG pathways (Additional file 4: Table S3).

ACAT1 is a fatty acid deposition gene that catalyzes the conversion of cholesterol to cholesteryl esters [56,57,58]. Yue et al. [58] selected the ACAT1 as a candidate gene to study adipogenesis in bovine adipose-derived mesenchymal stem cells. Like ACAT1, ACSL1 is crucial for the lipid metabolism, contributing to fatty acid biosynthesis, transport, storage and degradation; and, taking part in the mitochondrial beta-oxidation process [54, 59]. Zhao et al. [60] evaluating different tissues of Qinchuan cattle verified that ACSL1 mRNA was highly expressed in the skeletal muscle (Longissimus thoracis) and subcutaneous fat, affirming that this gene may contribute to the determination of fatty acid composition in bovine skeletal muscle. Recently, Poleti et al. [61] found ACSL1 protein as differentially abundant when studying the intramuscular fat deposition in the LD muscle of this Nelore population.

Another noteworthy gene verified in the PPAR pathway is SLC27A6 (Solute Carrier Family 27 Member 6), part of the solute carrier superfamily (SLC). The solute carrier family 27A (SLC27A) is a group of molecules ubiquitously expressed, involved in the lipid metabolism by transporting fatty acid proteins [62, 63]. According to Melo et al. [63], despite the SLC27A family is essential for body lipid distribution, the SLC27A1 was found more expressed in the skeletal muscle than adipose tissue, suggesting that this gene has a critical role in the absorption and storage of fatty acids by the muscle.

Beside the SLC27A1, we also identified other SLC members in the Salmon module – SLC25A3, SLC25A4, SLC25A11, SLC25A12 and, SLC26A9 – associated to three cellular components and two oxidative stress diseases pathways (Additional file 4: Table S3). In previous work with this Nelore population, we have already found SLC genes more expressed in the group with the lowest GEBV for REA [64]. Junior et al. [1], considered SLC38A1 and 2 as candidate genes for muscle growth associated with REA, in a GWAS study with Nelore cattle. In another way, Costa et al. [50] identified the SLC37A4 involved in relevant lipid metabolism biological functions, such as the concentration of lipid, metabolism of triacylglycerol and homeostasis of cholesterol, when studying the bovine Longissimus muscle. Thus, the enrichment analysis of the Salmon module demonstrates the complex regulation of skeletal muscle energy and lipid metabolism, helping us to understand molecular insights occurring in the bovine LD muscle that can influence REA.

Conclusions

With the construction of co-expressed gene modules, we verified that inflammation and immune response pathways and biological processes could modulate the BFT in Nelore cattle. RSAD2, EIF2AK2, and PARP genes could be considered as putative candidate genes for BFT trait. For REA, we found that energy and lipid metabolism, mainly the fatty acid metabolism, were the major pathways regulating this trait in this cattle population. We highlight ACAT1 and ACSL1 as putative candidate genes, associated with the energy and lipid metabolism in the skeletal muscle. These results allow us a better comprehension of the molecular mechanisms that are behind these economically important traits, that lead to muscle and fat deposition in bovine.

Methods

Animals and phenotypes

A total of 385 Nelore steers descending 34 unrelated bulls, representing the main pedigree lineages of Brazilian Nelore cattle, raised between 2009 to 2011, were used in this study. All animals were raised in the same nutritional and handling conditions and, finished in feedlot. More details were provided in [19].

The animals were slaughtered at an average age of 25 months in a commercial abattoir located in Bariri (São Paulo, Brazil) under Federal Inspection Service (SIF) supervision and, Brazilian Ministry of Agriculture, Livestock and Food Supply (MAPA) regularization. As mentioned in our previous research [64], 5 g of the Longissimus dorsi (LD) muscle (12th-13th ribs) was collected from the right side of the carcasses at the time of slaughter and stored in liquid nitrogen until RNA-Sequencing analyses. A sample of the LD muscle (10th-13th ribs) was excised at 24 h after slaughter from the left side of the carcasses and transported to the Embrapa Pecuária Sudeste Laboratory (São Carlos, São Paulo, Brazil) to measure REA and BFT. The REA was dimensioned with a grid of points (values presented in cm2) and, the BFT was measured with a graduated ruler (values shown in mm). All the experimental procedures were approved by the Institutional Animal Care and Use Committee Guidelines from Embrapa (approval code CEUA 01/2013).

To the co-expression network analysis, we selected 48 animals with contrasting genomic estimated breeding values (GEBV, Additional file 1: Table S1) from the total population of 385 animals, where 12 represents the higher values and 12 the lowest values of REA (mean and SD values = 3.14 ± 0.56 and − 2.94 ± 0.59) and BFT (mean and SD values = 1.24 ± 0.28 and − 0.91 ± 0.12), respectively. After that, the lists of the animals were combined, and the repeated ones were removed, totaling 43 animals. The GEBV were calculated by GenSel program [65], based in the SNP marker information, obtained by the BovineHD 770 k BeadChip (Infinium BeadChip, Illumina, CA, USA), as described in [66]. The a priori genetic and residual variance values were obtained from the Bayes C analysis where the genetic and a priori residual variance was equal to 1 [67]. A new Bayes C analysis was performed with the previous values for genetic and residual variances to estimate the GEBV values for each animal. Correlation analysis between the phenotypic values of REA and BFT x the GEBV for these traits and, the GEBV for REA x GEBV for BFT were performed in the R program.

RNA-sequencing

  • RNA extraction: Total RNA was extracted from approximately 100 mg of muscle tissue (n = 43) using Trizol reagent (Life Technologies, Carlsbad, CA, USA), according to the manufacturer’s instructions. The RNA integrity was verified using Bioanalyzer 2100 (Agilent, Santa Clara, CA, USA). All the samples presented an RNA Integrity Number (RIN) greater than 7.

  • Library preparation: A total of 2 μg of total RNA of each sample was used for library preparation following the TruSeq RNA Sample Preparation kit v2 guide (Illumina, San Diego, CA, USA) protocol. Libraries were quantified utilizing quantitative PCR with the KAPA Library Quantification kit (KAPA Biosystems, Foster City, CA, USA), and libraries mean size was estimated by the Bioanalyzer 2100.

  • Sequencing, quality control and alignment: Samples were diluted for the same concentration and grouped in pools. The sequencing flowcell lanes were clustered with the TruSeq PE Cluster kit v3-cBot-HS (Illumina, San Diego, CA, USA) and then, sequenced using HiSeq2500 ultra-high-throughput sequencing system (Illumina, San Diego, CA, USA) with the TruSeq SBS kit v3-HS. A more detailed description can be found in [47]. SeqClean software (https://sourceforge.net/projects/seqclean/files/) was employed to remove the adapters sequences used in the library preparation step, and low-complexity reads. For the quality control, FastQC version 0.10.1 software (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/) was applied. TopHat version 2.1.0 software [68] was used to map the read alignment against the reference genome Bos taurus UMD3.1 (http://www.ensembl.org/Bos_taurus/Info/Index/). Lastly, read counts for all annotated genes were calculated adopting HTSeq software version 0.6.1 (https://htseq.readthedocs.io/en/release_0.10.0/) [69]. Only read sequences that uniquely align to know chromosomes were used in this study.

Co-expression network analysis

To perform the co-expression network analysis, we used the WGCNA (Weighted Correlation Network Analysis) package from R [8], with RNA-Sequencing data (n = 43) with their counts normalized by transcript per million (TPM). After the data input, a cleaning and preprocessing step were done to remove outlier samples and genes with excessive numbers of missing entries. So, from a list of 15,631 genes inputted, 14,529 ones were used to construct gene networks with the “blockwiseModules” function.

First, a matrix of similarity was constructed by calculating Pearson correlations, to measure the similarity between the gene expression profiles of all the samples. Then, the similarity matrix was transformed into an adjacency matrix (A) raised to a β exponent (soft threshold) based on the free-scale topology criterion. In this study, the β parameter was equal to 6, and the free-scale topology was R2 = 0.80. The topological overlap matrix (TOM) was used to define modules based on dissimilarity (1-TOM). The minimum and maximum module size (genes per module) were five and 12,000, respectively. Modules were merged based on the dissimilarity between their eigengenes, which is the first principal component of each module and, represents the gene expression profile within the module [8]. For modules grouping, was used a threshold of 0.25 corresponding to a correlation of 0.75. Finally, for each gene module was assigned a color, genes not assembled to any modules were grouped in the Grey module.

Module-trait associations were estimated using the correlation between the module eigengene (ME) and the GEBV values of REA and BFT, allowing the identification of modules highly correlated with the interest traits. Genes of modules with significant module-trait associations (p-value< 0.1), for at least one trait, were assigned for functional enrichment analysis.

Functional enrichment analysis

An Overrepresentation Enrichment Analysis (ORA) was performed by DAVID (Database for Annotation, Visualization and Integrated Discovery) version 6.8 [70]; and, WebGestalt 2017 (WEB-based GEne SeT AnaLysis Toolkit) [71]. DAVID software revealed all Gene Ontology (GO) terms (BP, CC, and MF) and KEGG pathways of the co-expressed genes with a False Discovery Rate (FDR) of 5%. WebGestalt was used to find additional relevant KEGG pathways (FDR 5%) that may not appear in DAVID analysis, to better understand the biological mechanisms involved in each trait studied. The Bos taurus genome was used as background for both analyses.

Abbreviations

BFT:

Backfat thickness

BP:

Biological process

CC:

Cellular Component

DAVID:

Database for Annotation, Visualization and Integrated Discovery

ER:

Endoplasmic reticulum

FDR:

False discovery rate

GEBV:

Genomic estimated breeding value

GO:

Gene Ontology

IFN:

Interferons

KEGG:

Kyoto Encyclopedia of Genes and Genomes

LD:

Longissimus dorsi

ME:

Module eigengene

MF:

Molecular function

REA:

Ribeye area

TPM:

Transcript per million

WebGestalt:

WEB-based GEne SeT AnaLysis Toolkit

WGCNA:

Weighted Correlation Network Analysis

References

  1. Junior GA, Costa RB, de Camargo GM, Carvalheiro R, Rosa GJ, Baldi F, Garcia DA, Gordo DG, Espigolan R, Takada L, et al. Genome scan for postmortem carcass traits in Nellore cattle. J Anim Sci. 2016;94(10):4087–95.

    Article  Google Scholar 

  2. Prado CS, Pádua JT, Corrêa MPC, Ferraz JBS, Miyagi ES, Resende LS. Comparação de diferentes métodos de avaliação da área de olho de lombo e cobertura de gordura em bovinos de corte. Ciência Animal Brasileira. 2004;5(3):141–9.

    Google Scholar 

  3. Lopes LS, Ladeira MM, Machado Neto OR, Paulino PVR, Chizzotti ML, Ramos EM, Oliveira DM. Characteristics of carcass and commercial meat cuts from Red Norte and Nellore young bulls finished in feedlot. Rev Bras Zootec. 2012;41(4):970–7.

    Article  Google Scholar 

  4. Bonin MN, Ferraz JB, Eler JP, Rezende FM, Cucco DC, Carvalho ME, Silva RC, Gomes RC, Oliveira EC. Sire effects on carcass and meat quality traits of young Nellore bulls. Genet Mol Res. 2014;13(2):3250–64.

    Article  CAS  Google Scholar 

  5. Fan B, Du Z-Q, Gorbach DM, Rothschild MF. Development and application of high-density SNP arrays in genomic studies of domestic animals. Asian-Aust Anima Sci. 2010;23(7):833–47.

    Article  CAS  Google Scholar 

  6. Royer AM, Shivers C, Riley DG, Elzo MA, Garcia MD. Single nucleotide polymorphisms associated with carcass traits in a population of Brahman and Brahman-influenced steers. Genet Mol Res. 2016;15(2). https://doi.org/10.4238/gmr.15028280.

  7. Wickramasinghe S, Cánovas A, Rincón G, Medrano JF. RNA-sequencing: a tool to explore new frontiers in animal genetics. Livest Sci. 2014;166:206–16.

    Article  Google Scholar 

  8. Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics. 2008;9:559.

    Article  Google Scholar 

  9. DiLeo MV, Strahan GD, den Bakker M, Hoekenga OA. Weighted correlation network analysis (WGCNA) applied to the tomato fruit metabolome. PLoS One. 2011;6(10):e26683.

    Article  CAS  Google Scholar 

  10. Sabino M, Carmelo VAO, Mazzoni G, Cappelli K, Capomaccio S, Ajmone-Marsan P, Verini-Supplizi A, Trabalza-Marinucci M, Kadarmideen HN. Gene co-expression networks in liver and muscle transcriptome reveal sex-specific gene expression in lambs fed with a mix of essential oils. BMC Genomics. 2018;19(1):236.

    Article  Google Scholar 

  11. Fuller TF, Ghazalpour A, Aten JE, Drake TA, Lusis AJ, Horvath S. Weighted gene coexpression network analysis strategies applied to mouse weight. Mamm Genome. 2007;18(6-7):463–72.

    Article  Google Scholar 

  12. Heiland DH, Demerath T, Kellner E, Kiselev VG, Pfeifer D, Schnell O, Staszewski O, Urbach H, Weyerbrock A, Mader I. Molecular differences between cerebral blood volume and vessel size in glioblastoma multiforme. Oncotarget. 2017;8(7):11083–93.

    Article  Google Scholar 

  13. Ponsuksili S, Du Y, Hadlich F, Siengdee P, Murani E, Schwerin M, Wimmers K. Correlated mRNAs and miRNAs from co-expression and regulatory networks affect porcine muscle and finally meat properties. BMC Genomics. 2013;14:533.

    Article  CAS  Google Scholar 

  14. Ponsuksili S, Siengdee P, Du Y, Trakooljul N, Murani E, Schwerin M, Wimmers K. Identification of common regulators of genes in co-expression networks affecting muscle and meat properties. PLoS One. 2015;10(4):e0123678.

    Article  Google Scholar 

  15. Alexandre PA, Kogelman LJ, Santana MH, Passarelli D, Pulz LH, Fantinato-Neto P, Silva PL, Leme PR, Strefezzi RF, Coutinho LL, et al. Liver transcriptomic networks reveal main biological processes associated with feed efficiency in beef cattle. BMC Genomics. 2015;16:1073.

    Article  Google Scholar 

  16. Kong RS, Liang G, Chen Y, Stothard P, Guan le L. Transcriptome profiling of the rumen epithelium of beef cattle differing in residual feed intake. BMC Genomics. 2016;17:592.

    Article  Google Scholar 

  17. Chen Y, Liu Y, Du M, Zhang W, Xu L, Gao X, Zhang L, Gao H, Li J, Zhao M. Constructing a comprehensive gene co-expression based interactome in Bos taurus. PeerJ. 2017;5:e4107.

    Article  Google Scholar 

  18. Oliveira GB, Regitano LCA, Cesar ASM, Reecy JM, Degaki KY, Poleti MD, Felicio AM, Koltes JE, Coutinho LL. Integrative analysis of microRNAs and mRNAs revealed regulation of composition and metabolism in Nelore cattle. BMC Genomics. 2018;19(1):126.

    Article  Google Scholar 

  19. Tizioto PC, Decker JE, Taylor JF, Schnabel RD, Mudadu MA, Silva FL, Mourão GB, Coutinho LL, Tholon P, Sonstegard TS, et al. Genome scan for meat quality traits in Nelore beef cattle. Physiol Genomics. 2013;45(21):1012–20.

    Article  CAS  Google Scholar 

  20. Yokoo MJ, Lobo RB, Araujo FR, Bezerra LA, Sainz RD, Albuquerque LG. Genetic associations between carcass traits measured by real-time ultrasound and scrotal circumference and growth traits in Nelore cattle. J Anim Sci. 2010;88(1):52–8.

    Article  CAS  Google Scholar 

  21. Bonin MN, Ferraz JB, Pedrosa VB, Silva SL, Gomes RC, Cucco DC, Santana MH, Campos JH, Barbosa VN, Castro FS, et al. Visual body-scores selection and its influence on body size and ultrasound carcass traits in Nellore cattle. J Anim Sci. 2015;93(12):5597–606.

    Article  CAS  Google Scholar 

  22. Diniz FB, Villela SD, Mourthe MH, Paulino PV, Boari CA, Ribeiro JS, Barroso JA, Pires AV, Martins PG. Evaluation of carcass traits and meat characteristics of Guzerat-crossbred bulls. Meat Sci. 2016;112:58–62.

    Article  Google Scholar 

  23. Tao X, Liang Y, Yang X, Pang J, Zhong Z, Chen X, Yang Y, Zeng K, Kang R, Lei Y, et al. Transcriptomic profiling in muscle and adipose tissue identifies genes related to growth and lipid deposition. PLoS One. 2017;12(9):e0184120.

    Article  Google Scholar 

  24. Clímaco SM, Ribeiro ELA, Mizubuti IY, Silva LDF, Barbosa MAAF, Bridi AM. Desempenho e características de carcaça de bovino de corte de quatro grupos genéticos terminados em confinamento. Rev Bras Zootec. 2011;40(7):1562–7.

    Article  Google Scholar 

  25. Maia MO, Susin I, Pires AV, Gentil RS, Ferreira EM, Mendes CQ, Alencar SM. Growth, carcass characteristics, chemical composition and fatty acid profile of the longissimus dorsi muscle in goat kids fed diets with castor oil. R Bras Zootec. 2012;41(11):2343–9.

    Article  Google Scholar 

  26. Schroder UJ, Staufenbiel R. Invited review: methods to determine body fat reserves in the dairy cow with special regard to ultrasonographic measurement of backfat thickness. J Dairy Sci. 2006;89(1):1–14.

    Article  CAS  Google Scholar 

  27. Yin Z, Deng T, Peterson LE, Yu R, Lin J, Hamilton DJ, Reardon PR, Sherman V, Winnier GE, Zhan M, et al. Transcriptome analysis of human adipocytes implicates the NOD-like receptor pathway in obesity-induced adipose inflammation. Mol Cell Endocrinol. 2014;394(1-2):80–7.

    Article  CAS  Google Scholar 

  28. Nakamura T, Furuhashi M, Li P, Cao H, Tuncman G, Sonenberg N, Gorgun CZ, Hotamisligil GS. Double-stranded RNA-dependent protein kinase links pathogen sensing with stress and metabolic homeostasis. Cell. 2010;140(3):338–48.

    Article  CAS  Google Scholar 

  29. Gregor MF, Hotamisligil GS. Inflammatory mechanisms in obesity. Annu Rev Immunol. 2011;29:415–45.

    Article  CAS  Google Scholar 

  30. Wensveen FM, Valentic S, Sestan M, Turk Wensveen T, Polic B. The “big bang” in obese fat: events initiating obesity-induced adipose tissue inflammation. Eur J Immunol. 2015;45(9):2446–56.

    Article  CAS  Google Scholar 

  31. Warfel JD, Bermudez EM, Mendoza TM, Ghosh S, Zhang J, Elks CM, Mynatt R, Vandanmagsar B. Mitochondrial fat oxidation is essential for lipid-induced inflammation in skeletal muscle in mice. Sci Rep. 2016;6:37941.

    Article  CAS  Google Scholar 

  32. Sanyal A, Naumann J, Hoffmann LS, Chabowska-Kita A, Ehrlund A, Schlitzer A, Arner P, Bluher M, Pfeifer A. Interplay between obesity-induced inflammation and cGMP signaling in white adipose tissue. Cell Rep. 2017;18(1):225–36.

    Article  CAS  Google Scholar 

  33. Severa M, Coccia EM, Fitzgerald KA. Toll-like receptor-dependent and -independent viperin gene expression and counter-regulation by PRDI-binding factor-1/BLIMP1. J Biol Chem. 2006;281(36):26188–95.

    Article  CAS  Google Scholar 

  34. Martire S, Navone ND, Montarolo F, Perga S, Bertolotto A. A gene expression study denies the ability of 25 candidate biomarkers to predict the interferon-beta treatment response in multiple sclerosis patients. J Neuroimmunol. 2016;292:34–9.

    Article  CAS  Google Scholar 

  35. Hinson ER, Cresswell P. The antiviral protein, viperin, localizes to lipid droplets via its N-terminal amphipathic alpha-helix. Proc Natl Acad Sci U S A. 2009;106(48):20452–7.

    Article  CAS  Google Scholar 

  36. Dogan A, Lasch P, Neuschl C, Millrose MK, Alberts R, Schughart K, Naumann D, Brockmann GA. ATR-FTIR spectroscopy reveals genomic loci regulating the tissue response in high fat diet fed BXD recombinant inbred mouse strains. BMC Genomics. 2013;14:386.

    Article  CAS  Google Scholar 

  37. Ramos-Lopez O, Riezu-Boj JI, Milagro FI, Martinez JA. DNA methylation signatures at endoplasmic reticulum stress genes are associated with adiposity and insulin resistance. Mol Genet Metab. 2018;123(1):50–8.

    Article  CAS  Google Scholar 

  38. Udumula MP, Babu MS, Bhat A, Dhar I, Sriram D, Dhar A. High glucose impairs insulin signaling via activation of PKR pathway in L6 muscle cells. Biochem Biophys Res Commun. 2017;486(3):645–51.

    Article  CAS  Google Scholar 

  39. Jokinen R, Pirnes-Karhu S, Pietilainen KH, Pirinen E. Adipose tissue NAD(+)-homeostasis, sirtuins and poly(ADP-ribose) polymerases -important players in mitochondrial metabolism and metabolic health. Redox Biol. 2017;12:246–63.

    Article  CAS  Google Scholar 

  40. Vida A, Marton J, Miko E, Bai P. Metabolic roles of poly(ADP-ribose) polymerases. Semin Cell Dev Biol. 2017;63:135–43.

    Article  CAS  Google Scholar 

  41. Mohamed JS, Hajira A, Pardo PS, Boriek AM. MicroRNA-149 inhibits PARP-2 and promotes mitochondrial biogenesis via SIRT-1/PGC-1alpha network in skeletal muscle. Diabetes. 2014;63(5):1546–59.

    Article  CAS  Google Scholar 

  42. Jin C, Flavell RA. Innate sensors of pathogen and stress: linking inflammation to obesity. J Allergy Clin Immunol. 2013;132(2):287–94.

    Article  CAS  Google Scholar 

  43. Costa D, Abreu JBR, Mourão RC, Silva JCG, Rodrigues VC, Sousa JCD, Marques RAFS. Características de Carcaça de Novilhos Inteiros Nelore e F1 Nelore x Holandês. Ciência Animal Brasileira. 2007;8(4):685–94.

    Google Scholar 

  44. Deshmukh AS, Murgia M, Nagaraj N, Treebak JT, Cox J, Mann M. Deep proteomics of mouse skeletal muscle enables quantitation of protein isoforms, metabolic pathways, and transcription factors. Mol Cell Proteomics. 2015;14(4):841–53.

    Article  CAS  Google Scholar 

  45. Zhang L, Keung W, Samokhvalov V, Wang W, Lopaschuk GD. Role of fatty acid uptake and fatty acid beta-oxidation in mediating insulin resistance in heart and skeletal muscle. Biochim Biophys Acta. 2010;1801(1):1–22.

    Article  CAS  Google Scholar 

  46. van Hall G. The physiological regulation of skeletal muscle fatty acid supply and oxidation during moderate-intensity exercise. Sports Med. 2015;45(Suppl 1):S23–32.

    Article  Google Scholar 

  47. Cesar AS, Regitano LC, Poleti MD, Andrade SC, Tizioto PC, Oliveira PS, Felicio AM, do Nascimento ML, Chaves AS, Lanna DP, et al. Differences in the skeletal muscle transcriptome profile associated with extreme values of fatty acids content. BMC Genomics. 2016;17(1):961.

    Article  Google Scholar 

  48. Marrades MP, Gonzalez-Muniesa P, Arteta D, Martinez JA, Moreno-Aliaga MJ. Orchestrated downregulation of genes involved in oxidative metabolic pathways in obese vs. lean high-fat young male consumers. J Physiol Biochem. 2011;67(1):15–26.

    Article  CAS  Google Scholar 

  49. Jeong J, Bong J, Kim GD, Joo ST, Lee HJ, Baik M. Transcriptome changes favoring intramuscular fat deposition in the longissimus muscle following castration of bulls. J Anim Sci. 2013;91(10):4692–704.

    Article  CAS  Google Scholar 

  50. Costa ASH, Costa P, Alves SP, Alfaia CM, Prates JAM, Vleck V, Cassar-Malek I, Hocquette JF, Bessa RJB. Does growth path influence beef lipid deposition and fatty acid composition? PLoS One. 2018;13(4):e0193875.

    Article  Google Scholar 

  51. Xu WD, Yang XY, Li DH, Zheng KD, Qiu PC, Zhang W, Li CY, Lei KF, Yan GQ, Jin SW, et al. Up-regulation of fatty acid oxidation in the ligament as a contributing factor of ankylosing spondylitis: a comparative proteomic study. J Proteome. 2015;113:57–72.

    Article  CAS  Google Scholar 

  52. Ehrenborg E, Krook A. Regulation of skeletal muscle physiology and metabolism by peroxisome proliferator-activated receptor delta. Pharmacol Rev. 2009;61(3):373–93.

    Article  CAS  Google Scholar 

  53. Doran AG, Berry DP, Creevey CJ. Whole genome association study identifies regions of the bovine genome and biological pathways involved in carcass trait performance in Holstein-Friesian cattle. BMC Genomics. 2014;15:837.

    Article  Google Scholar 

  54. Huang W, Guo Y, Du W, Zhang X, Li A, Miao X. Global transcriptome analysis identifies differentially expressed genes related to lipid metabolism in wagyu and Holstein cattle. Sci Rep. 2017;7(1):5278.

    Article  Google Scholar 

  55. Lim D, Chai HH, Lee SH, Cho YM, Choi JW, Kim NK. Gene expression patterns associated with peroxisome proliferator-activated receptor (PPAR) signaling in the longissimus dorsi of Hanwoo (Korean cattle). Asian-Australas J Anim Sci. 2015;28(8):1075–83.

    Article  CAS  Google Scholar 

  56. Chang TY, Li BL, Chang CCY, Urano Y. Acyl-coenzyme a:cholesterol acyltransferases. Am J Physiol Endocrinol Metab. 2009;297(1):E1–9.

    Article  CAS  Google Scholar 

  57. Sakashita N, Lei X, Kamikawa M, Nishitsuji K. Role of ACAT1-positive late endosomes in macrophages: cholesterol metabolism and therapeutic applications for Niemann-pick disease type C. J Med Investig. 2014;61(3-4):270–7.

    Article  Google Scholar 

  58. Yue Y, Zhang L, Zhang X, Li X, Yu H. De novo lipogenesis and desaturation of fatty acids during adipogenesis in bovine adipose-derived mesenchymal stem cells. In Vitro Cell Dev Biol Anim. 2018;54(1):23–31.

    Article  Google Scholar 

  59. Coleman RA, Lewin TM, Muoio DM. Physiological and nutritional regulation of enzymes of triacylglycerol synthesis. Annu Rev Nutr. 2000;20:77–103.

    Article  CAS  Google Scholar 

  60. Zhao ZD, Zan LS, Li AN, Cheng G, Li SJ, Zhang YR, Wang XY, Zhang YY. Characterization of the promoter region of the bovine long-chain acyl-CoA synthetase 1 gene: roles of E2F1, Sp1, KLF15, and E2F4. Sci Rep. 2016;6:19661.

    Article  CAS  Google Scholar 

  61. Poleti MD, Regitano LCA, Souza GHMF, Cesar ASM, Simas RC, Silva-Vignato B, Oliveira GB, Andrade SCS, Cameron LC, Coutinho LL. Longissimus dorsi muscle label-free quantitative proteomic reveals biological mechanisms associated with intramuscular fat deposition. J Proteome. 2018;179:30–41.

    Article  CAS  Google Scholar 

  62. Gallardo D, Amills M, Quintanilla R, Pena RN. Mapping and tissue mRNA expression analysis of the pig solute carrier 27A (SLC27A) multigene family. Gene. 2013;515(1):220–3.

    Article  CAS  Google Scholar 

  63. Melo C, Gallardo D, Quintanilla R, Zidi A, Castelló A, Díaz I, Amills M, Pena RN. An association analysis between polymorphisms of the pig solute carrier family 27A (SLC27A), member 1 and 4 genes and serum and muscle lipid traits. Livest Sci. 2013;152(2):143–6.

    Article  Google Scholar 

  64. Silva-Vignato B, Coutinho LL, Cesar ASM, Poleti MD, Regitano LCA, Balieiro JCC. Comparative muscle transcriptome associated with carcass traits of Nellore cattle. BMC Genomics. 2017;18(1):506.

    Article  Google Scholar 

  65. Fernando RL, Garrick DJ. GenSel – user manual for a portfolio of genomic selection related analyses. 3rd ed. Iowa: Iowa State University; 2009. p. 1–24.

    Google Scholar 

  66. Cesar AS, Regitano LC, Mourao GB, Tullio RR, Lanna DP, Nassu RT, Mudado MA, Oliveira PS, do Nascimento ML, Chaves AS, et al. Genome-wide association study for intramuscular fat deposition and composition in Nellore cattle. BMC Genet. 2014;15:39.

    Article  Google Scholar 

  67. Kizilkaya K, Garrick DJ, Fernando RL, Mestav B, Yildiz MA. Use of linear mixed models for genetic evaluation of gestation length and birth weight allowing for heavy-tailed residual effects. Genet Sel Evol. 2010;42:26.

    Article  Google Scholar 

  68. Kim D, Pertea G, Trapnell C, Pimentel H, Kelley R, Salzberg SL. TopHat2: accurate alignment of transcriptomes in the presence of insertions, deletions and gene fusions. Genome Biol. 2013;14(4):R36.

    Article  Google Scholar 

  69. Anders S, Pyl PT, Huber W. HTSeq—a Python framework to work with high-throughput sequencing data. Bioinformatics. 2015;31(2):166–9.

    Article  CAS  Google Scholar 

  70. Huang DW, Sherman BT, Lempicki RA. Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat Protoc. 2009;4(1):44–57.

    Article  CAS  Google Scholar 

  71. Wang J, Vasaikar S, Shi Z, Greer M, Zhang B. WebGestalt 2017: a more comprehensive, powerful, flexible and interactive gene set enrichment analysis toolkit. Nucleic Acids Res. 2017;45(W1):W130–7.

    Article  CAS  Google Scholar 

Download references

Acknowledgements

We thank Embrapa and University of São Paulo for the collaborative efforts, CAPES (Coordination for the Improvement of Higher Level Education) for the scholarship to BSV, and FAPESP (São Paulo Research Foundation) for providing financial support (process number 2012/23638-8).

Funding

The funding institutions were Embrapa (Macroprograma 1, 01/2005) and FAPESP (São Paulo Research Foundation, process number 2012/23638-8). LCAR and LLC were granted CNPq fellowships.

Availability of data and materials

The dataset supporting the conclusions of this article is available in the European Nucleotide Archive (ENA) repository (EMBL-EBI), under accession PRJEB13188 [https://www.ebi.ac.uk/ena/submit/sra/].

Author information

Authors and Affiliations

  1. College of Agriculture “Luiz de Queiroz”, University of São Paulo, Piracicaba, SP, 13418-900, Brazil

    Bárbara Silva-Vignato, Luiz L. Coutinho & Aline S. M. Cesar

  2. College of Animal Science and Food Engineering, University of São Paulo, Pirassununga, SP, 13635-900, Brazil

    Mirele D. Poleti

  3. Federal University of Lavras, Lavras, MG, 37200-000, Brazil

    Cristina T. Moncau

  4. Embrapa Pecuária Sudeste, São Carlos, SP, 13560-970, Brazil

    Luciana C. A. Regitano

  5. College of Veterinary Medicine and Animal Science, University of São Paulo, Pirassununga, SP, 13635-900, Brazil

    Júlio C. C. Balieiro

Authors
  1. Bárbara Silva-Vignato
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Luiz L. Coutinho
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Mirele D. Poleti
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Aline S. M. Cesar
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Cristina T. Moncau
    View author publications

    You can also search for this author in PubMed Google Scholar

  6. Luciana C. A. Regitano
    View author publications

    You can also search for this author in PubMed Google Scholar

  7. Júlio C. C. Balieiro
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

BSV, MDP and ASMC participated in the experimental design. BSV, MDP and ASMC performed data analysis. BSV drafted the manuscript. BSV, MDP, ASMC, CTM, JCCB, LCAR and LLC collaborated with interpretation and discussion of the results. LCAR, LLC and JCCB provided the experimental environment and coordination. All authors have read and approved the final manuscript.

Corresponding author

Correspondence to Bárbara Silva-Vignato.

Ethics declarations

Ethics approval

The animals (n = 385) used in this study came from herds of the Brazilian Agricultural Research Corporation (Embrapa). All the experimental procedures involving animals were approved by the Institutional Animal Care and Use Committee Guidelines from Embrapa (approval code CEUA 01/2013).

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Publisher’s Note

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

Additional files

Additional file 1:

Table S1. Test of means (t-test) of backfat thickness (BFT) and ribeye area (REA) between groups with High (H) and Low (L) genomic estimated breeding values (GEBV) for REA and BFT in the Longissimus dorsi muscle of Nelore steers. (XLS 35 kb)

Additional file 2:

Figure S1. Heatmap plot of the gene network using a subset of 400 genes. The heatmap plot depicts the Topological Overlap Matrix (TOM) among a subset of 400 genes from the analysis. Each row and column represents a single gene. The light colors represent the low overlap between modules, progressively darker red color represents higher overlap. Darker color blocks along the diagonal represent gene modules. The gene dendrogram and module assignment are shown above and along the left side of the graph. (PNG 107 kb)

Additional file 3:

Table S2. Functional enrichment analysis from the gene list of the Green Yellow module, performed by DAVID v6.8 (FDR < 0.05). The table contains the Gene Ontology category, identification and description, p-value adjusted for a false discovery rate of 5%, nominal p-value, number of genes and gene names for each category. (XLS 30 kb)

Additional file 4:

Table S3. Functional enrichment analysis from the gene list of the Salmon module, performed by DAVID v6.8 (FDR < 0.05). The table contains the Gene Ontology category, identification and description, p-value adjusted for a false discovery rate of 5%, nominal p-value, number of genes and gene names for each category. (XLS 38 kb)

Additional file 5:

Table S4. KEGG Pathways (FDR < 0.05) identified by WebGestalt 2017 from the gene list of the Salmon module. The table contains the KEGG Pathway identification, description, number of genes and gene names for each pathway. (XLS 31 kb)

Additional file 6:

Table S5. Gene list from the Ivory module. The table contains the Ensembl gene identification, gene symbol and gene description of the entire list of genes from the Ivory module. (XLS 28 kb)

Additional file 7:

Table S6. Gene list from the Light Yellow module. The table contains the Ensembl gene identification, gene symbol and gene description of the entire list of genes from the Light Yellow module. (XLS 38 kb)

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), 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 (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Silva-Vignato, B., Coutinho, L.L., Poleti, M.D. et al. Gene co-expression networks associated with carcass traits reveal new pathways for muscle and fat deposition in Nelore cattle. BMC Genomics 20, 32 (2019). https://doi.org/10.1186/s12864-018-5345-y

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12864-018-5345-y

Keywords

BMC Genomics

ISSN: 1471-2164

Contact us