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

Integrative analysis of microRNAs and mRNAs revealed regulation of composition and metabolism in Nelore cattle

BMC Genomics volume 19, Article number: 126 (2018) Cite this article

Abstract

Background

The amount of intramuscular fat can influence the sensory characteristics and nutritional value of beef, thus the selection of animals with adequate fat deposition is important to the consumer. There is growing knowledge about the genes and pathways that control the biological processes involved in fat deposition in muscle. MicroRNAs (miRNAs) belong to a well-conserved class of non-coding small RNAs that modulate gene expression across a range of biological functions in animal development and physiology. The aim of this study was to identify differentially expressed (DE) miRNAs, regulatory candidate genes and co-expression networks related to intramuscular fat (IMF) deposition. To achieve this, we used mRNA and miRNA expression data from the Longissimus dorsi muscle of 30 Nelore steers with high (H) and low (L) genomic estimated breeding values (GEBV) for IMF deposition.

Results

Differential miRNA expression analysis between animals with extreme GEBV values for IMF identified six DE miRNAs (FDR 10%). Functional annotation of the target genes for these microRNAs indicated that the PPARs signaling pathway is involved with IMF deposition. Candidate regulatory genes such as SDHAF4, FBXO17, ALDOA and PKM were identified by partial correlation with information theory (PCIT), phenotypic impact factor (PIF) and regulatory impact factor (RIF) co-expression approaches from integrated miRNA-mRNA expression data. Two DE miRNAs (FDR 10%), bta-miR-143 and bta-miR-146b, which were upregulated in the Low IMF group, were correlated with regulatory candidate genes, which were functionally enriched for fatty acid oxidation GO terms. Co-expression patterns obtained by weighted correlation network analysis (WGCNA), which showed possible interaction and regulation between mRNAs and miRNAs, identified several modules related to immune system function, protein metabolism, energy metabolism and glucose catabolism according to in silico analysis performed herein.

Conclusion

In this study, several genes and miRNAs were identified as candidate regulators of IMF by analyzing DE miRNAs using two different miRNA-mRNA co-expression network methods. This study contributes to the understanding of potential regulatory mechanisms of gene signaling networks involved in fat deposition processes measured in muscle. Glucose metabolism and inflammation processes were the main pathways found in silico to influence intramuscular fat deposition in beef cattle in the integrative mRNA-miRNA co-expression analysis.

Background

The amount of intramuscular fat (IMF) is an important characteristic associated with juiciness and taste of beef [1]. IMF deposition is associated with size and number of adipocytes, the balance between lipogenesis and lipolysis rate in muscle and changes in catabolic activities in different species [2,3,4,5]. In humans, skeletal muscle insulin resistance can be associated with IMF deposition [6], or more recently with lipid intermediates [7]. Overall meat quality can be impacted by many factors such as nutritional program, environment, age, sex and genetics. However, little attention has been paid to the role of microRNAs in the regulation of IMF deposition in cattle.

MicroRNAs (miRNAs) are endogenous non-coding (ncRNA) ribonucleic acids (RNAs) that are approximately twenty-two nucleotides in length [8]. These molecules modulate the expression of genes at the post-transcriptional level by blocking the translation of target mRNAs [8, 9]. MiRNAs play an important role in post-transcriptional gene regulation in many tissues and are associated with the control of several important biological processes related to lipid metabolism [10]. Understanding the regulatory functions of miRNA and other small RNAs on the expression of target genes that impact lipogenesis and adipogenesis is important to identify target molecules that influence fat deposition. Several studies have been published, which demonstrate the importance of miRNAs as potential biomarkers for variations in subcutaneous adipose tissue [11,12,13]. However, limited information about the importance of miRNA is available for IMF [14, 15]. Once identified, biomarkers could be used in animal breeding programs to improve meat quality and animal productivity [11, 13] and potentially contribute to our understanding of insulin resistance associated with human diseases such as obesity and diabetes [16, 17].

Although RNA-seq analyses can be helpful for genomic studies and can generate lists of expressed genes in specific tissues to ultimately detect differentially expressed (DE) genes, the biological interpretation of this data is still a challenge. Network approaches that integrate data have proven useful in the identification of complex transcriptional regulation. For example, hub genes, which are highly correlated with a large number of genes, have been shown to have key regulatory roles in gene expression networks [18,19,20]. Thus, co-expression analysis may be more sensitive at detecting biologically interesting pathways than analysis of DE genes expression [21]. Several network approaches are available for this purpose, such as the Weighted Gene Co-expression Network Analysis (WGCNA) and the Partial Correlation with Information Theory (PCIT) methods. The WGCNA method identifies gene correlation networks, i.e. gene clusters of biological significance, from expression profiling data [22]. The PCIT method identifies differences in pairs of correlated gene expression levels to measure a gene’s differential connectivity across levels of a phenotype [23]. Utilization of both PCIT and WGCNA have enabled a better understanding of the co-regulation of mRNAs and miRNAs for different phenotypes [21, 24,25,26,27] to better comprehend the biological mechanisms and regulatory processes in lipid metabolism.

In this study, skeletal muscle microRNA and mRNA expression data from animals with different IMF deposition were integrated with two well documented systems biology methodologies [22, 23]. These analyses indicate that miRNAs play a role in IMF deposition by modulating carbohydrate, lipid and immune response metabolic pathways in skeletal muscle.

Results

Phenotypic and sequencing data

The genetic variance, residual variance and heritability for intramuscular fat (IMF) obtained from this population were 0.196, 0.490 and 0.29 ± 0.16, respectively, as previously published [28]. The animals were ranked using genomic estimated breeding values (GEBV) for IMF values and fifteen animals with high IMF GEBV (H) and fifteen with low IMF GEBV (L) were selected for miRNA-Seq analysis (Table 1). This strategy, to select animals with extreme GEBV, was performed because the correlation (r) between the raw IMF values (% IMF) and GEBV was high (r = 0.76) [28] and the statistical T-test showed that the GEBV averages for groups were statistically different (p-value = 2.2e-16).

Table 1 Phenotypic data for intramuscular fat percentage (IMF), genomic estimated breeding values (GEBV) and number of normalized mapped miRNA reads for all animals
Full size table

A total of 32 million (M) sequence reads were obtained from an Illumina MiSeq. The average number of total reads per sample was one million. The total number of mapped reads was 24 M with an average of 84% reads mapped (Table 1).

Differentially expressed microRNAs and target genes identification

Twenty-six novel and 463 known miRNAs were identified using miRDeep2. Among all microRNAs identified, six of them were differentially expressed (DE) with a False Discovery Rate (FDR) of 10% (Table 2). Negative values of fold change indicate lower miRNA expression in animals with low IMF deposition and positive values indicate higher miRNA expression for this group. These six microRNAs targeted 2250 genes expressed in skeletal muscle based on IPA analysis (Additional file 1: Table S1). Of note, because bta-let-7f and bta-let-7a-5p belong to the same family of miRNAs and have the same seed sequence, they most likely target the same genes (Table 2).

Table 2 List of differentially expressed miRNAs between Low and High groups, based on genomic estimated breeding values (GEBV) for intramuscular fat, identified by miRDeep2 and the number of target genes obtained by IPA® for each miRNA
Full size table

Enrichment analysis of target genes from the DE microRNAs

Functional enrichment analyses of target genes by IPA showed networks and canonical pathways related to fatty acid metabolism (Table 3). Gene networks and the principal canonical pathways are described in Additional file 2: Table S2, Additional file 3: Table S3, Additional file 4: Figures S1, S2 and S3, and Additional file 5: Figures S4, S5 and S6. The most relevant gene network was “lipid metabolism, small molecule biochemistry, vitamin and mineral metabolism” that involved genes such as PPARGC1A, MYCN, ESR2 and ARL4D, that are targets of downregulated miRNAs and MED1, SMAD4, NEDD4 and MBOAT2, that are targets of upregulated miRNAs in the L group (Fig. 1).

Table 3 List of the top gene networks and signaling pathways related with lipid metabolism identified by IPA®
Full size table
Fig. 1
figure 1

Gene network “lipid Metabolism, Small Molecule Biochemistry, Vitamin and Mineral Metabolism” identified from the DE miRNA’s target genes list generated by IPA. Grey shapes represent target genes and the white shapes are other genes of the network that are not target genes. Solid lines mean direct interaction and dashed lines an indirect interaction between genes

Full size image

Target genes enriched for PPAR-RXR signaling pathways (i.e. lipogenesis promoting) were negatively regulated by miRNAs which were upregulated in L group. Target genes associated with fatty acid oxidation were targets of downregulated miRNAs. Some important genes for lipid metabolism present in this pathway included: PPARa, PKA and ADIPOR2. These genes are targets of the downregulated miRNAs bta-let-7 and bta-miR-423 (i.e. downregulated in the L group). On the other hand, STAT5b and GPDH are targets of upregulated miRNAs (bta-miR-100 and bta-miR-143) in L group (Fig. 2).

Fig. 2
figure 2

The PPARα-RXR signaling pathway is over-represented in miRNA target genes identified by IPA. The shapes highlighted in purple represent the miRNA target genes and the white shapes represent the other genes of the pathway that are not target genes

Full size image

Co-expression analysis: PCIT - differential Hubbing

After data filtering by expression in high (H) IMF and low (L) IMF groups, a list with 383 miRNAs and 14,650 genes expressed in bovine skeletal muscle were used for Partial Correlation with Information Theory (PCIT) analysis, which allowed for the identification of ten positive and negative differentially hubbed (DH) genes and miRNAs (Additional file 6: Table S4; Table 4).

Table 4 List of the top ten positive and negative Differentially Hubbed (DH) genes and miRNAs, when comparing High and Low groups of GEBV for IMF
Full size table

The genes with a significant correlation with DH genes were used to construct co-expression networks and identify enriched GO terms (Additional file 7: Tables S5 and S6, Additional file 8: Tables S7 and S8). The GO terms enriched among all genes correlated to the top ten negative DH genes were most related to glucose metabolism (GO ID: 6006, GO ID: 6007, GO ID: 6096) (Fig. 3) and for the top ten positive DH genes the GO terms were related to protein and mRNA metabolism (GO ID: 6364, GO ID: 6350, GO ID: 30,163, GO ID: 30,162, GO ID: 51,603) (Fig. 4). The DE miRNAs bta-miR-143 and bta-miR-146b were upregulated in the low (L) IMF group (Table 2), were positively correlated with negative DH genes, and were associated with glucose and fatty acids catabolism (Fig. 3). The most important DH genes potentially involved in the regulation of lipid metabolism and protein metabolism are shown in Table 5. The co-expression networks of top DH genes were visualized by BioLayout (Fig. 5).

Fig. 3
figure 3

GO terms enriched from genes significantly correlated with negatively differentially hubbed (DH) genes. Bubble color indicates the user-provided p-value and the lower the p-value of the processes grouped in the category, the more intense is the coloring of bubbles; bubble size indicates the frequency of the GO term in the underlying GOA database. Highly similar GO terms are linked by edges in the graph, where the line width indicates the degree of similarity

Full size image
Fig. 4
figure 4

GO terms enriched from genes significantly correlated with positively differentially hubbed (DH) genes. Bubble color indicates the user-provided p-value and the lower the p-value of the processes grouped in the category, the more intense is the coloring of bubbles; bubble size indicates the frequency of the GO term in the underlying GOA database. Highly similar GO terms are linked by edges in the graph, where the line width indicates the degree of similarity

Full size image
Table 5 List of the top two differentially hubbed (DH) genes and the GO terms associated with them. The negative DH genes have higher number of connections in Low GEBV group and positive DH genes in High GEBV group
Full size table
Fig. 5
figure 5

Negative (a) and positive (b) differentially hubbed (DH) genes associated with lipid metabolism between the High GEBV and Low GEBV IMF groups. The blue edges represent negative correlations between hub gene and the correlated genes while the red edges represent positive correlations

Full size image

Co-expression analysis: Phenotype impact factor (PIF) and regulatory impact factor (RIF)

The PIF and RIF analyses were used to identify putative regulatory genes that may explain differences in phenotype between groups of animals, based on differences in gene expression (Additional file 9: Tables S9, S10 and S11). The most relevant genes for fatty acid metabolism found in the top RIF and PIF analyses and the GO terms associated with them are shown in Table 6. The genes with positive values for RIF 1 and 2 were PYGM, ENO3, ATP2A1, GAPDH and ALDOA, which were principally related to glucose metabolism and energy metabolism. Those with negative values of RIF were mostly miRNAs. Specifically miRNA bta-miR-143, which was identified as DE in this study (Table 1), was a miRNA with one of the most negative RIF2 value. The PIF analysis also identified ALDOA as a putative regulatory gene for the difference in fat deposition between H and L groups. The enrichment of GO terms was similar with those found in DH analysis (Additional file 8: Tables S7 and S8; Additional file 10: Table S12).

Table 6 List of the genes with the most extreme Phenotypic Impact Factor (PIF) and Regulatory Impact Factor (RIF) 1 and 2 values and the GO terms associated with them
Full size table

Co-expression analysis: WGCNA - miRNAs correlated with mRNA modules

The weighted correlation network analysis (WGCNA) methodology was applied in two different manners, first to integrate mRNAs and miRNAs by analyzing those modules (i.e. co-expression networks) that had a negative correlation between them and second by identifying modules that are important to phenotypic variation by correlating all modules identified in WGCNA with intramuscular fat (IMF) deposition. A total of 27 mRNA modules in high (H) IMF and 44 in the low (L) IMF group were identified. Furthermore, there were 14 miRNA modules in H and 22 in L. The grey module contained all genes not included in a correlated module (Additional file 11: Tables S13 and S14; Additional file 12: Tables S15 and S16).

After correlating all miRNA and mRNA modules with one another, those modules that were negatively correlated with one another were investigated further. Among all correlated modules, three miRNA modules were negatively correlated with five mRNA modules in the H group, while six miRNA modules were negatively correlated with seven mRNA modules in the L group (p-value > 0.05; Additional file 13: Table S17). The genes that composed each mRNA module were significantly over enriched for GO terms related to lipid metabolism (adj. p-value < 0.1; Additional file 14: Tables S18 and S19). These lipid metabolism GO terms that were enriched from mRNA modules were then used to construct mRNA-miRNA co-expressed networks for both groups H and L (Figs. 6 and 7). The miRNAs modules were enriched for GO terms based on the hub miRNA target genes of each module (Additional file 14: Tables S20 and S21).

Fig. 6
figure 6

Co-expression networks showing the negative correlation among miRNAs and biological processes enriched in mRNA modules in High IMF GEBV group. Colored circles represent hub miRNAs, with higher connectivity, inside each module and squares represent the GO terms associated with each mRNA module, represented by different letter color

Full size image
Fig. 7
figure 7

Co-expression networks show the negative correlation among miRNAs and biological processes enriched in mRNA modules in Low IMF GEBV group. Colored circles represent hub miRNAs, with higher connectivity, inside each module and squares represent the GO terms associated with each mRNA module, represented by different letter color

Full size image

The WGCNA results indicated that the cyan and lightgreen mRNA modules in the high (H) IMF group were enriched for inflammatory response and adaptive immune response (Additional file 14: Table S18). They were also negatively correlated with the black and pink miRNA modules (Fig. 6), whose target genes were associated with carboxylic acid transport and positive regulation of leukocyte migration (Additional file 14: Table S20). Unlike the high (H) group, the low (L) IMF group contained multiple miRNA modules that were negatively correlated with the same mRNA module (Additional file 13: Table S17, Fig. 7). The target genes of the most connected miRNA module (midnightblue) in the L group were associated with metal ion homeostasis and regulation of interleukin-1 beta (IL-1B) production (Additional file 14: Table S21). This module was negatively correlated with four mRNA modules (green, orange, yellowgreen and black) (Fig. 7). These four mRNA modules were enriched for protein kinase B (PKB) signaling cascade, positive regulation of fatty acid oxidation and negative regulation of protein ubiquitination (Additional file 14: Table S19).

The hub miRNAs in the green module from the high IMF group (e.g. bta-miR-106a, bta-miR-671) (Fig. 6) and midnightblue module from the low IMF group (e.g. bta-miR-122, bta-miR-1291) (Fig.7) identified in the WGCNA, were also found correlated with the DH genes in the PCIT analysis (Additional file 7: Tables S5 and S6). Additionally, the miR-1291 was both a hub in the WGCNA midnightblue module from low (L) IMF group (Fig. 7) and a negatively DH gene (Table 4). These results indicate that WGCNA and PCIT not only identify similar biological processes associated with fat deposition, but they can also identify the same miRNAs in co-expression networks.

Co-expression analysis: WGCNA – Modules correlated with phenotype

Besides the integration of mRNA and miRNA data, the correlation of modules with the phenotype (high or low IMF deposition) was also performed (Additional file 15: Figure S7; Additional file 16: Figure S8). Three mRNA modules in the high (H) IMF group and two in the low (L) IMF group, and three miRNA modules both the H and L groups were correlated with IMF (Table 7). In the H group, the black miRNA module and the cyan mRNA module were negatively correlated with each other (Additional file 13: Table S17).

Table 7 GO term enrichment of modules significantly correlated with the IMF trait, for mRNAs and miRNAs in the High (H) and Low (L) groups
Full size table

The enrichment analysis results of modules correlated with IMF are presented in Table 7. The table shows the biological processes downregulated by miRNAs modules correlated with phenotype and the biological processes associated with mRNAs modules correlated with phenotype in both groups. The GO terms associated with each mRNA and miRNA module that was significantly correlated with IMF are presented in Additional file 17: Tables S22, S23, S24 and S25.

The hub miRNAs in modules that were correlated with intramuscular fat (IMF) in both H and L groups downregulate processes related with lipid metabolism. The mRNA modules that were correlated with IMF in both groups were enriched for inflammatory response and immune system GO terms (Table 7).

Discussion

The regulation of lipid oxidation and biosynthesis is under strong feedback control in order to maintain homeostasis [29]. Intramuscular fat deposition in cattle is of economic importance and has been studied in several breeds, including Nelore under different nutritional conditions [30]. However, there is still limited knowledge on the molecular processes, in particular regarding microRNA involvement in the control of this trait [31]. It is important to remember that samples collected in this study are from the Longissimus dorsi (LD) muscle and thus contains myofibers and other cell types such as adipocytes, fibroblast, blood and nerve cells. Thus, the gene expression profile obtained can be influenced by the proportion of the different cell types in the sample and/or by genes being differentially expressed in a specific cell type.

Networks enriched for lipid and carbohydrate metabolism

Although research studying the influence of miRNAs on metabolism has been recently published [32], there has been a rapid growth in the number of identified miRNAs that are involved in the regulation of genes and signaling molecules responsible for maintaining lipid homeostasis [26, 33].

The comparison of miRNA expression between animals with different genetic potential for intramuscular fat (IMF) deposition resulted in the identification of six differentially expressed (DE) miRNAs involved in fatty acid metabolism and lipid deposition. The miR-423 and let-7 family, upregulated in the high (H) IMF group, have been previously reported to be associated with obesity in humans and disorders in glucose metabolism in mice [16, 34] and have been implicated as possible biomarkers for risk of obesity [16, 34]. The miRNAs upregulated in the low (L) IMF group, miR-100, miR-146 and miR-143, have been reported to control aspects of adipogenesis in humans [35,36,37,38]. By in vitro analysis Chen et al. [35] suggested that overexpression of mir-143 in adipose-derived stem cells (ADSCs) in rats could promote or inhibit adipogenesis by regulation of MAPK signaling pathway depending on the stage of development. Interestingly, they found that upregulation of miR-143 expression in early stages of adipogenesis blocks adipocyte differentiation, but when it happens later induces clonal expansion of adipose tissue.

The enrichment analysis of DE miRNAs’ target genes revealed that the PPAR pathway was overrepresented in IPA analysis (Additional file 5: Figures S4, S5 and S6). Furthermore, many important target genes related to lipid metabolism were present in the gene networks identified by IPA (Table 3 and Additional file 4: Figures S1, S2 and S3). Among them, PPARs were found, which are a class of ligand-activated transcription factors that have a well-known influence on lipid metabolism and glucose homeostasis [39,40,41,42,43]. The target genes of miRNAs downregulated in the low (L) IMF group were enriched for the GO term fatty acid oxidation involved in PPARa pathway. On the other hand, based on IPA enrichment results, miRNAs upregulated in this group would downregulate genes of lipogenesis and adipogenesis.

In a previous comparison of gene expression profile between Angus (higher IMF deposition capacity) and Nelore, Martins et al. (2015) [30] did not observed a difference in PPARγ mRNA expression; however, a higher amount of PPARγ protein was detected by western blot. This result highlights the limitation of individual mRNA gene expression analysis and shows the importance of a systems biology approach, in which co-expression analysis can indicate a modulation in a pathway even without detectable difference in individual gene expression.

To gain additional insights into the pathways impacted by miRNAs, we performed co-expression analysis by integration of the miRNA and mRNA data. The PCIT analysis revealed that the top negatively differentially hubbed (DH) genes, which had more connections in the low (L) IMF group, were correlated with genes associated mostly with carbohydrate metabolism (Table 5). The DH genes that may be the most relevant for IMF were SDHAF4 and bta-miR-24. SDHAF4 is important for the assembly of succinate dehydrogenase and plays a role in ATP synthesis by the electron transport chain. Both SDHAF4 and bta-miR-24 are associated with carbohydrate metabolism and could indicate that a change in myofiber type is associated with IMF. However, in our previous study in which we measured DE genes associated with IMF in these same animals [28], there was no significant difference in myosin heavy chain (MyHC) isoform expression. Interestingly, miR-24 has been shown to negatively regulate adipocyte differentiation and hepatic lipid accumulation in mice [44, 45].

The candidate regulatory genes identified by PIF and RIF that negatively regulate intramuscular fat (IMF) deposition were PKM, bta-miR-143 and bta-miR-26b. PKM is associated with glucose metabolism, while bta-miR-26b was related to control cholesterol efflux and lipogenesis in mice [46, 47] (Table 6). The target genes of bta-miR-143 were enriched for glutamate catabolism. Glutamate is a key component in cellular metabolism, and it is related to biosynthesis of lipids, because it is utilized in the citric acid cycle to produce ATP through α-ketoglutarate [48]. MiR-143 downregulates this process by blocking excess ATP production that could induce storage of lipids instead of undergoing lipid degradation. This co-expression analysis reaffirms the importance of the bta-miR-143 in control of fat deposition.

WGCNA revealed that the mRNA module in the high (H) IMF group that was positively correlated with IMF deposition (darkorange) was enriched for vitamin B6 metabolic process (Table 7), which is indirectly related to lipogenesis. Several enzyme reactions involved in fatty acid metabolism require vitamin B-6 as a coenzyme, such as the biosynthesis of sphingolipids [49], which are a class of lipids that are components of cell membranes. Moreover, the black and yellow miRNA modules in the high (H) IMF group were positively correlated with IMF deposition. The target genes for the miRNA in these modules are associated with carboxylic acid transport and regulation of translation (Table 7). Fatty acids are carboxylic acids and their transport into the mitochondria leads to activation of β-oxidation to produce energy. This indicates that miRNAs associated with high fat deposition are downregulating biological processes such as transport and catabolism of fatty acids, while miRNAs negatively associated with higher fat deposition downregulate glycolipid degradation. Interestingly in the high (H) IMF group, the black miRNA module and cyan mRNA module were negatively correlated with each other (Additional file 13: Table S17) and both were differently correlated with IMF, positively and negatively, respectively (Table 7). These miRNAs were associated with downregulation of fatty acid transport, while the mRNAs were associated with inflammation.

In the low (L) IMF group, the greenyellow and lightgreen miRNA modules were positively correlated with IMF. The miRNA in these modules downregulate genes enriched for lipid metabolic process and calcium ion homeostasis, respectively (Table 7). Calcium (Ca) participates in many signaling networks that contribute to modulation of enzyme function, including Ca-sensitive enzymes involved in lipolysis and lipogenesis [50]. The purple miRNA module was negatively correlated with IMF. Most of the miRNAs in this module were expressed at lower levels in lean animals. The target genes of these miRNA were associated with lipid storage. Overall, the miRNAs in co-expression networks associated with low IMF were related to lipid metabolism, lipolysis, lipogenesis and lipid storage.

Networks related to immune system and inflammatory response

It is known that lipid accumulation in obesity activates the immune system which leads to an inflammatory state due to secretion of proinflammatory molecules by adipocytes [51]. Genes associated with inflammatory response were identified as target genes of DE miRNAs that were upregulated in the low (L) IMF group, which was enriched for the PPAR-RXR signaling pathway (Fig. 2). These genes mediate signal transduction from members of the interleukin-1 (IL-1) family. IL-1, which is regulated by PPARa, can induce and regulate a network of proinflammatory cytokines that initiate inflammatory responses [52]. Using the same population of animals as utilized here, Cesar et al. [28] previously reported that DE genes were associated with inflammatory response.

The WGCNA results of co-expressed mRNA-miRNA modules in the high (H) IMF group (Fig. 6) were enriched for GO terms associated with inflammatory response and adaptive immune response. These same modules were indirectly correlated with carboxylic acid transport and leukocyte migration (Additional file 13: Table S17). Adaptive immune cells have been reported to be increased in obese mice and humans, which can trigger a sequence of proinflammatory reactions and could be associated with impaired glucose tolerance and insulin sensitivity [53,54,55,56]. This indicates that mRNA and miRNA co-expression networks constructed for animals with differences in lipid accumulation (high IMF group) may be involved in pathways that regulate the immune system and inflammation and are correlated with lipid metabolism.

Conclusion

A combination of DE and co-expression based analyses indicate that lipid metabolism, glucose metabolism and inflammatory response are the main biological processes associated with IMF deposition in Nelore cattle. The miRNAs identified in this study were not only associated with extreme intramuscular fat deposition levels, but also participate in co-expression networks that may affect mRNA expression and metabolic pathways modulating fat deposition. Furthermore, both co-expression approaches could construct similar miRNAs networks that were correlated with genes and pathways important for IMF accumulation. This study allowed us to better understand the potential roles of miRNA regulation and interaction in fat deposition and revealed new candidate regulatory genes and miRNAs associated with IMF.

Methods

Animals and phenotypic data

Genotypic and phenotypic data were collected on 310 Nelore steers sired by 34 unrelated sires that represent the main breeding lineages in Brazilian Nelore from an experimental breeding herd from EMBRAPA between 2009 and 2011 [57]. The animals were raised in feedlots under identical nutrition and handling conditions until slaughter at an average age of 25 months. Samples from Longissimus dorsi (LD) muscle located between the 12th and 13th ribs were collected at two time points: at slaughter for RNA sequencing analysis, and 24 h after slaughter for the intramuscular fat (IMF) deposition measurement as described below [28].

Approximately 100 g of muscle were lyophilized and ground to measure IMF deposition using an Ankom XT20 extractor and the AOCS procedure (official Procedure Am 5–04) as described Cesar et al. [57]. Animals with extreme values for intramuscular fat (IMF) deposition were selected based on their genomic estimated breeding values (GEBV) [57]. GEBV was predicted by Genomic Best Linear Unbiased Prediction (GBLUP) methodology, which was conducted using ASREML software [58]. A group of 30 animals were selected (fifteen with high IMF GEBV values and fifteen with low IMF GEBV values) for mRNA and miRNA analyses.

RNA extraction and RNA-sequencing

Total RNA was isolated from 100 mg of LD muscle samples from 30 steers with extreme GEBV values. The extraction of total RNA was performed using the Trizol reagent (Invitrogen) according to the protocol described by Chomczynski and Sacchi [59]. After extraction, total RNA was quantified by spectrophotometer (NanoDrop 200 - Thermo Scientific. Wilmington. Delaware, USA). The integrity of the RNA was verified by size separation on a 1% agarose gel and analysis on a Bioanalyzer 2100 (Agilent Technologies - Santa Clara, CA, USA) with the RNA 6000 Nano kit. All samples had an RNA integrity number (RIN) greater than or equal to 8. Then samples were diluted to a final concentration of 200 ng/μL. Sequencing libraries were generated with the TruSeq® smallRNA Sample Preparation kit (Illumina - San Diego, USA). The concentration of the cDNA libraries was determined with the KAPA Library Quantification Kit (KAPA Biosystems) and then samples were sequenced on a Miseq machine (Illumina), using MiSeq Reagent Kit v3 (150 cycles), generating around 1 million reads/sample.

Reads filtering and miRNAs identification

After sequencing, data quality was evaluated with FastQC [60] and filtered by Phred score quality using FASTX-Toolkit [61] software, where the minimum quality Phred score was 28. Then, the miRDeep2 [62] program was used to identify and quantify miRNAs, using the default parameters. The sequences were mapped against the bovine reference genome Bos taurus UMD 3.1 and compared with miRBase database (v. 21) [63].

Differentially expressed miRNAs

In order to identify differentially expressed (DE) miRNAs between the L and H groups, the total count data of each miRNA was analyzed with the DESeq2 package [64], using a statistical model that fitted contemporary group (animal origin and year that the animal enter the experiment) as a categorical fixed effect and age at slaughter of an animal as a covariate. To remove variation due to the preparation of sequencing libraries, the expression data were normalized by library size, as described in the manual of the DESeq2 package [64]. The Benjamini-Hochberg (BH) [65] methodology was used to control the False Discovery Rate (FDR) of DE at 10%.

Identification of miRNA target genes and enrichment analysis of DE miRNAs

The miRNA target genes were obtained from the MicroRNA Target Filter tool of QIAGEN’s Ingenuity Pathway Analysis (IPA®, Redwood City-CA) that uses TargetScan, miRecords and TarBase as the miRNA target genes databases. The IPA® uses information from predicted targets of mammalian microRNAs, based on the fact that target sites are usually conserved because miRNAs are highly evolutionarily conserved [66,67,68,69]. After this first approach, to obtain the target genes by IPA, we selected only those target genes that were expressed in our animals, using the data generated by RNA-seq of skeletal muscle samples of the same animals from previous study [28]. The functional enrichment of target genes was also performed by IPA software to identify enriched metabolic pathways and gene networks associated with lipid metabolism.

PCIT and differential hubbing (DH) network analysis

To improve the functional annotation of miRNA and mRNA interactions in a systems biology context, the Partial Correlation with Information Theory (PCIT) analyses [70, 71] were conducted on the combined list of miRNAs (383) and mRNAs (14,650) after normalization of expression level by DESeq2. The miRNAs and mRNAs were filtered to select only those expressed in animals in both H and L IMF groups. The mRNA expression data utilized in this study was previously published by Cesar et al. [28].

PCIT was used to evaluate the specific behavior or co-expression between all miRNAs and genes and from this information, differential connectivity or hubbing (DH) [71] was calculated. Differential hubbing is the difference in the number of significant partial correlations (connections) a gene has between two different treatments, in this case compared between H and L groups and filtering those correlations higher than 0.9. BioLayout Express3D [72] software was used to visualize gene networks.

PIF and RIF analysis

To identify putative candidate regulators responsible for the differences observed in phenotypes, the Phenotypic Impact Factor (PIF) and Regulatory Impact Factor (RIF) approaches were performed [23, 70, 73]. PIF gives a ‘weight’ for the contribution and importance of genes to the differences involved between phenotypes, based exclusively on their numerical properties. RIF is based on differences in the regulator’s correlations and it represents the relative importance of genes/miRNAs on the phenotypically relevant part of the network. The RIF1 value is based largely on changes in correlation between two treatments levels (i.e. differential wiring). The RIF2 value allows genes to be ranked as potential regulators based on the expression changes of a regulator and how it can affect the expression of other genes in the network due to treatment differences [23].

WGCNA

The same list of genes and miRNAs used in the PCIT analysis was utilized to run the R package WGCNA [22]. This analysis constructs clusters of highly correlated genes and miRNAs in modules and allows the correlation of them to each other and to a trait (i.e. intramuscular fat (IMF) deposition). In contrast to the analysis performed with PCIT, the WGCNA was done separately for genes and then for miRNAs.

Modules of mRNA

To construct clusters of genes, pair-wise Pearson correlation coefficients were first calculated between all expressed transcripts to generate a signed similarity. To emphasize (weight) stronger correlations and punish weaker correlations, the signed similarity matrix was then raised to the lowest power β that approximated a scale-free network topology (R2 > 0.90) to generate an adjacency matrix [74]. The β’s used to construct the mRNA modules from the L and H IMF group’s expression data were 12 and 8, respectively (Additional file 18: Figure S9). The topological overlap distance calculated from the adjacency matrix is then clustered with the average linkage hierarchical clustering. The default minimum cluster merge height of 0.25 was retained. The clusters created by WGCNA are called modules, and the minimum number of genes in a module was set to 30. Each module represents a group of highly interconnected genes with similar expression profiles across the samples and the expression profile pattern is distinct from those of other modules [22]. Modules were named by a conventional color scheme and genes not classified in a correlated module were grouped in the grey module. After modules were defined, the module Eigengene (MEs) values were calculated. The Eigengene of a module is defined as the eigenvector associated with the first principal component of the expression matrix representing the expression profile of all genes within a given module [22, 75].

Modules of miRNA

The steps for constructing miRNA co-expression modules were as described above. After generating the signed similarity matrix, a power β value was chosen to generate the adjacency matrix, the β’s used to construct the miRNA modules were 9 for the L and 4 for the H IMF groups (Additional file 19: Figure S10). The topological overlap distance was calculated and a minimum module size of five miRNAs was chosen. Five was chosen as the minimum module size for the miRNAs due to the smaller size of the miRNA transcriptome relative to the mRNA transcriptome [22, 75].

Correlation between mRNA and miRNA modules

An integrative analysis was performed correlating the ME of miRNAs with the ME of mRNAs, for each group. Those modules with a negative correlation higher than − 0.4 with a p-value < 0.05 were used for enrichment analysis. The co-expression networks among hub miRNAs, representing the whole module, and the GO terms of mRNAs inside the correlated modules were constructed in Cytoscape v.3.3.0 0 [76].

Correlation of modules with phenotype

Using the ME, the Module-Trait relationships were estimated by calculating the Pearson’s correlations between the ME and the animals’ phenotypic information (i.e. % IMF) to select potential biologically interesting modules that could explain the phenotypic differences between groups. To avoid losing information and expand biological response, we set a p-value threshold of correlation of 0.1 to select modules correlated with trait.

MiRNA target gene identification and enrichment analysis of co-expression data

The general gene enrichment of GO terms for biological processes was made using BiNGO (Biological Networks Gene Ontology), tool for Cytoscape v.3.3.0 [76] and REVIGO [77] to visualize clusters of GO terms. The Benjamini-Hochberg (BH) [65] methodology was used as a multiple testing correction to control the False Discovery Rate (FDR) at 10%. For miRNAs, the combined results from miRanda and TargetScan approaches were used to identify the target genes and these genes were also filtered by skeletal muscle RNA-seq data of previous study [28] to do the enrichment. MiRanda is a method for target site identification from sequence information [78]. It compares the miRNAs complementarity to 3’UTR regions of genome. Using a perl script we generated a fasta file with all 3’UTR regions of bovine genome (UMD 3.1) from ENSEMBL to use as input in miRanda. The TargetScan was performed for mammals and customized by species (cow/Bos taurus) (http://www.targetscan.org/vert_70/).

The enrichment of miRNA modules identified by WGCNA was conducted using the target genes information of specific hub miRNAs in each miRNA module. In this case the hub miRNAs were those with the highest Modular Membership (MM) value for the module, which means that these miRNAs have higher connectivity inside the module and are probably more informative [18].

Abbreviations

BH:

Benjamin-Hochberg

DE:

Differentially expressed

DH:

Differentially hubbed

FC:

Fold Change

FDR:

False discovery rate

GBLUP:

Genomic best linear unbiased prediction

GEBV:

Genomic estimated breeding values

GO:

Gene ontology

GOA:

Gene ontology annotation

H:

High IMF group

IMF:

Intramuscular fat

IPA:

Ingenuity pathways analysis

L:

Low IMF group

LD:

Longissimus dorsi

ME:

Module eigengene

MM:

Modular membership

ncRNA:

non-coding RNA

PCIT:

Partial correlation with information theory

PIF:

Phenotypic impact factor

RIF:

Regulatory impact factor

RNA:

Ribonucleic acid

WGCNA:

Weighted correlation network analysis

References

  1. Hocquette JF, Gondret F, Baéza E, Médale F, Jurie C, Pethick DW. Intramuscular fat content in meat-producing animals: development, genetic and nutritional control, and identification of putative markers. Animal. 2010;4(2):303–19.

    Article  CAS  PubMed  Google Scholar 

  2. Cianzo D, Topel D, Whitehurst G, Beitz D, Self H. Adipose tissue growth and cellularity: changes in bovine adipocyte and number. J Anim Sci. 1985;60:970–6.

    Article  Google Scholar 

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

    Article  CAS  PubMed  Google Scholar 

  4. Hernández P, Navarro JL, Toldrá F. Effect of frozen storage on lipids and lipolytic activities in the longissimus dorsi muscle of the pig. Z Lebensmittel Zitschrift fur Lebensmittel Untersuchung und Forschung. 1999;208:110–5.

    Article  Google Scholar 

  5. Martínez-Álvaro M, Agha S, Blasco A, Hernández P. Muscle lipid metabolism in two rabbit lines divergently selected for intramuscular fat. J Anim Sci. 2017;95(6):2576–84.

    PubMed  Google Scholar 

  6. Nawrocki AR, Scherer PE. The delicate balance between fat and muscle: adipokines in metabolic disease and musculoskeletal inflammation. Curr Opin Pharmacol. 2004;4(3):281–9.

    Article  CAS  PubMed  Google Scholar 

  7. Kitessa MS, Abeywardena YM. Lipid-induced insulin resistance in skeletal muscle: the chase for the culprit Goes from Total intramuscular fat to lipid intermediates, and finally to species of lipid intermediates. Nutrients. 2016;8(8):466-79.

  8. Bartel DP. MicroRNAs: genomics, biogenesis, mechanism, and function. Cell. 2004;116(2):281–97.

    Article  CAS  PubMed  Google Scholar 

  9. Lee RC, Feinbaum RL, Ambros V. The C. Elegans heterochronic gene lin-4 encodes small RNAs with antisense complementarity to lin-14. Cell. 1993;75(5):843–54.

    Article  CAS  PubMed  Google Scholar 

  10. Coutinho LL, Matukumalli LK, Sonstegard TS, Van Tassell CP, Gasbarre LC, Capuco AV, Smith TPL. Discovery and profiling of bovine microRNAs from immune-related and embryonic tissues. Physiol Genomics. 2007;29(1):35.

    Article  CAS  PubMed  Google Scholar 

  11. Romao JM, Jin W, He M, McAllister T, Guan LL. MicroRNAs in bovine adipogenesis: genomic context, expression and function. BMC Genomics. 2014;15(1):137.

    Article  PubMed  PubMed Central  Google Scholar 

  12. Jin W, Dodson MV, Moore SS, Basarab JA, Guan LL. Characterization of microRNA expression in bovine adipose tissues: a potential regulatory mechanism of subcutaneous adipose tissue development. BMC Mol Biol. 2010;11(1):29.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  13. Sun J, Zhang B, Lan X, Zhang C, Lei C, Chen H. Comparative transcriptome analysis reveals significant differences in MicroRNA expression and their target genes between adipose and muscular tissues in cattle. PLoS One. 2014;9(7):e102142.

    Article  PubMed  PubMed Central  Google Scholar 

  14. Wang H, Zheng Y, Wang G, Li H. Identification of microRNA and bioinformatics target gene analysis in beef cattle intramuscular fat and subcutaneous fat. Mol BioSyst. 2013;9:2154–62.

    Article  CAS  PubMed  Google Scholar 

  15. Guan L, Hu X, Liu L, Xing Y, Zhou Z, Liang X, Yang Q, Jin S, Bao J, Gao H, et al. Bta-miR-23a involves in adipogenesis of progenitor cells derived from fetal bovine skeletal muscle. Sci Rep. 2017;7:43716.

    Article  PubMed  PubMed Central  Google Scholar 

  16. Ortega FJ, Mercader JM, Catalán V, Moreno-Navarrete JM, Pueyo N, Sabater M, Gómez-Ambrosi J, Anglada R, Fernández-Formoso JA, Ricart W, et al. Targeting the circulating MicroRNA signature of obesity. Clin Chem. 2013;59(5):781.

    Article  CAS  PubMed  Google Scholar 

  17. Li H, Zhang Z, Zhou X, Wang Z, Wang G, Han Z. Effects of MicroRNA-143 in the differentiation and proliferation of bovine intramuscular preadipocytes. Mol Biol Rep. 2011;38(7):4273–80.

    Article  CAS  PubMed  Google Scholar 

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

    Article  CAS  PubMed  Google Scholar 

  19. Lim D, Lee S-H, Kim N-K, Cho Y-M, Chai H-H, Seong H-H, Kim H. 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.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

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

    Google Scholar 

  21. Schroyen M, Steibel JP, Koltes JE, Choi I, Raney NE, Eisley C, Fritz-Waters E, Reecy JM, Dekkers JCM, Rowland RRR, et al. Whole blood microarray analysis of pigs showing extreme phenotypes after a porcine reproductive and respiratory syndrome virus infection. BMC Genomics. 2015;16(1):516.

    Article  PubMed  PubMed Central  Google Scholar 

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

    Article  PubMed  PubMed Central  Google Scholar 

  23. Reverter A, Hudson NJ, Nagaraj SH, Pérez-Enciso M, Dalrymple BP. Regulatory impact factors: unraveling the transcriptional regulation of complex traits from expression data. Bioinformatics. 2010;26(7):896–904.

    Article  CAS  PubMed  Google Scholar 

  24. 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(1):533.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  25. Nunez YO, Truitt JM, Gorini G, Ponomareva ON, Blednov YA, Harris RA, Mayfield RD. Positively correlated miRNA-mRNA regulatory networks in mouse frontal cortex during early stages of alcohol dependence. BMC Genomics. 2013;14(1):725.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  26. Yang X, Koltes JE, Park CA, Chen D, Reecy JM. Gene co-expression network analysis provides novel insights into Myostatin regulation at three different mouse developmental Timepoints. PLoS One. 2015;10(2):e0117607.

    Article  PubMed  PubMed Central  Google Scholar 

  27. Buchanan JW, Reecy JM, Garrick DJ, Duan Q, Beitz DC, Koltes JE, Saatchi M, Koesterke L, Mateescu RG. Deriving gene networks from SNP associated with triacylglycerol and phospholipid fatty acid fractions from Ribeyes of Angus cattle. Front Genet. 2016;7:116.

    Article  PubMed  PubMed Central  Google Scholar 

  28. Cesar ASM, Regitano LCA, Koltes JE, Fritz-Waters ER, Lanna DPD, Gasparin G, Mourão GB, Oliveira PSN, Reecy JM, Coutinho LL. Putative regulatory factors associated with intramuscular fat content. PLoS One. 2015;10(6):e0128350.

    Article  PubMed  PubMed Central  Google Scholar 

  29. Raghow R, Yellaturu C, Deng X, Park EA, Elam MB. SREBPs: the crossroads of physiological and pathological lipid homeostasis. Trends Endocrinol Metab. 2008;19(2):65–73.

    Article  CAS  PubMed  Google Scholar 

  30. Martins TS, Sanglard LMP, Silva W, Chizzotti ML, Rennó LN, Serão NVL, Silva FF, Guimarães SEF, Ladeira MM, Dodson MV, et al. Molecular factors underlying the deposition of intramuscular fat and collagen in skeletal muscle of Nellore and Angus cattle. PLoS One. 2015;10(10):e0139943.

    Article  PubMed  PubMed Central  Google Scholar 

  31. Baik M, Kang HJ, Park SJ, Na SW, Piao M, Kim SY, Fassah DM, Moon YS. TRIENNIAL GROWTH AND DEVELOPMENT SYMPOSIUM: molecular mechanisms related to bovine intramuscular fat deposition in the longissimus muscle. J Anim Sci. 2017;95(5):2284–303.

    CAS  PubMed  Google Scholar 

  32. Xu P, Vernooy SY, Guo M, Hay BA. The drosophila MicroRNA Mir-14 suppresses cell death and is required for normal fat metabolism. Curr Biol. 2003;13(9):790–5.

    Article  CAS  PubMed  Google Scholar 

  33. Rottiers V, Näär AM. MicroRNAs in metabolism and metabolic disorders. Nat Rev Mol Cell Biol. 2012;13(4):239–50.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  34. Zhu H, Shyh-Chang N, Segrè Ayellet V, Shinoda G, Shah Samar P, Einhorn William S, Takeuchi A, Engreitz Jesse M, Hagan John P, Kharas Michael G, et al. The Lin28/let-7 Axis regulates glucose metabolism. Cell. 2011;147(1):81–94.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  35. Chen L, Hou J, Ye L, Chen Y, Cui J, Tian W, Li C, Liu L. MicroRNA-143 regulates Adipogenesis by modulating the MAP2K5–ERK5 signaling. Sci Rep. 2014;4:3819.

    Article  PubMed  PubMed Central  Google Scholar 

  36. Chen L, Dai Y-M, Ji C-B, Yang L, Shi C-M, Xu G-F, Pang L-X, Huang F-Y, Zhang C-M, Guo X-R. MiR-146b is a regulator of human visceral preadipocyte proliferation and differentiation and its expression is altered in human obesity. Mol Cell Endocrinol. 2014;393(1–2):65–74.

    Article  CAS  PubMed  Google Scholar 

  37. Esau C, Kang X, Peralta E, Hanson E, Marcusson EG, Ravichandran LV, Sun Y, Koo S, Perera RJ, Jain R, et al. MicroRNA-143 regulates adipocyte differentiation. J Biol Chem. 2004;279(50):52361–5.

    Article  CAS  PubMed  Google Scholar 

  38. Ortega FJ, Moreno-Navarrete JM, Pardo G, Sabater M, Hummel M, Ferrer A, Rodriguez-Hermosa JI, Ruiz B, Ricart W, Peral B, et al. MiRNA expression profile of human subcutaneous adipose and during adipocyte differentiation. PLoS One. 2010;5(2):e9022.

    Article  PubMed  PubMed Central  Google Scholar 

  39. Evans RM, Barish GD, Wang Y-X. PPARs and the complex journey to obesity. Nat Med. 2004;10(4):355–61.

    Article  CAS  PubMed  Google Scholar 

  40. Lefebvre P, Chinetti G, Fruchart J-C, Staels B. Sorting out the roles of PPARα in energy metabolism and vascular homeostasis. J Clin Investig. 2006;116(3):571–80.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  41. Koh E, Kim M, Park J, Kim H, Youn J, Park H, Youn J, Lee K. Peroxisome proliferator-activated receptor (PPAR)-α activation prevents diabetes in OLETF rats. Diabetes. 2003;42:2331–7.

  42. Stienstra R, Duval C, Müller M, Kersten S. PPARs, Obesity, and Inflammation. PPAR Research. 2007;2007:10.

  43. Harrington W, Britt C, Wilson J, Milliken N, Binz J, Lobe D, Oliver W, Lewis M, Ignar D. The effect of PPARα, PPARδ, PPARγ, and PPARpan agonists on BodyWeight, BodyMass, and serum lipid profiles in diet-induced obese AKR/JMice. PPAR Res. 2007;2007:1-13.

  44. Ng R, Wu H, Xiao H, Chen X, Willenbring H, Steer CJ, Song G: Inhibition of MicroRNA-24 Expression in Liver Prevents Hepatic Lipid Accumulation and Hyperlipidemia. Hepatology (Baltimore, Md) 2014, 60(2):554–564.

  45. Kang M, Yan LM, Li YM, Zhang WY, Wang H, Tang AZ, Ou HS. inhibitory effect of microRNA-24 on fatty acid-binding protein expression on 3T3-L1 adipocyte differentiation. Genetics and molecular research : GMR. 2013;12:5267–77.

    Article  CAS  PubMed  Google Scholar 

  46. Sun D, Zhang J, Xie J, Wei W, Chen M, Zhao X. MiR-26 controls LXR-dependent cholesterol efflux by targeting ABCA1 and ARL7. FEBS Lett. 2012;586(10):1472–9.

    Article  CAS  PubMed  Google Scholar 

  47. Fu X, Dong B, Tian Y, Lefebvre P, Meng Z, Wang X, Pattou F, Han W, Wang X, Lou F, et al. MicroRNA-26a regulates insulin sensitivity and metabolism of glucose and lipids. J Clin Invest. 2015;125(6):2497–509.

    Article  PubMed  PubMed Central  Google Scholar 

  48. Newsholme P, Lima MMR, Procopio J, Pithon-Curi TC, Doi SQ, Bazotte RB, Curi R. Glutamine and glutamate as vital metabolites. Braz J Med Biol Res. 2003;36:153–63.

    Article  CAS  PubMed  Google Scholar 

  49. Combs G. Vitamin B6. In: The Vitamins: Fundamental Aspects in Nutrition and Health. San Diego: Elsevier; 2008. p. 313–29.

    Google Scholar 

  50. Clapham DE. Calcium Signaling. Cell. 2007;131(6):1047–58.

    Article  CAS  PubMed  Google Scholar 

  51. Lau DCW, Dhillon B, Yan H, Szmitko PE, Verma S. Adipokines: molecular links between obesity and atheroslcerosis. Am J Physiol Heart Circ Physiol. 2005;288(5):H2031.

    Article  CAS  PubMed  Google Scholar 

  52. Dinarello CA. Interleukin-1 in the pathogenesis and treatment of inflammatory diseases. Blood. 2011;117(14):3720.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  53. Winer S, Chan Y, Paltser G, Truong D, Tsui H, Bahrami J, Dorfman R, Wang Y, Zielenski J, Mastronardi F, et al. Normalization of obesity-associated insulin resistance through immunotherapy. Nat Med. 2009;15(8):921–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  54. Rausch ME, Weisberg S, Vardhana P, Tortoriello DV. Obesity in C57BL/6J mice is characterized by adipose tissue hypoxia and cytotoxic T-cell infiltration. Int J Obes. 2007;32(3):451–63.

    Article  Google Scholar 

  55. Nishimura S, Manabe I, Nagasaki M, Eto K, Yamashita H, Ohsugi M, Otsu M, Hara K, Ueki K, Sugiura S, et al. CD8+ effector T cells contribute to macrophage recruitment and adipose tissue inflammation in obesity. Nat Med. 2009;15(8):914–20.

    Article  CAS  PubMed  Google Scholar 

  56. Sun S, Ji Y, Kersten S, Qi L. Mechanisms of inflammatory responses in obese adipose tissue. Annu Rev Nutr. 2012;32:261–86.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  57. Cesar ASM, Regitano LCA, Mourão GB, Tullio RR, Lanna DPD, Nassu RT, Mudado MA, Oliveira PSN, do Nascimento ML, Chaves AS, et al. Genome-wide association study for intramuscular fat deposition and composition in Nellore cattle. BMC Genet. 2014;15(1):39.

    Article  PubMed  PubMed Central  Google Scholar 

  58. GILMOUR A, GOGEL B, CULLIS B, THOMPSON R. ASReml User Guide Release 3.0. In: Hemel Hempstead. HP1 1ES UK: VSN International Ltd; 2009.

    Google Scholar 

  59. Chomczynski P, Sacchi N. Single-step method of RNA isolation by acid guanidinium thiocyanate-phenol-chloroform extraction. Anal Biochem. 1987;162(1):156–9.

    Article  CAS  PubMed  Google Scholar 

  60. ANDREWS S. FastQC: A quality control tool for High Throughput Sequence Data. 2010. http://www.bioinformatics.babraham.ac.uk/projects/fastqc.

  61. FASTX-Toolkit [http://hannonlab.cshl.edu/fastx_toolkit/].

  62. Friedländer MR, Mackowiak SD, Li N, Chen W, Rajewsky N. miRDeep2 accurately identifies known and hundreds of novel microRNA genes in seven animal clades. Nucleic Acids Res. 2012;40(1):37–52.

    Article  PubMed  Google Scholar 

  63. Griffiths-Jones S, Grocock RJ, van Dongen S, Bateman A, Enrigh AJ. miRBase: microRNA sequences, targets and gene nomenclature. Nucleic Acids Res. 2006;34:D140–4.

    Article  CAS  PubMed  Google Scholar 

  64. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15(12):550.

    Article  PubMed  PubMed Central  Google Scholar 

  65. Benjamini Y, Hochberg Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J R Stat Soc. 1995;57:289-300.

  66. Friedman RC, KK-H F, Burge CB, Bartel DP. most mammalian mRNAs are conserved targets of microRNAs. Genome Res. 2009;19(1):92–105.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  67. Naeem A, Zhong K, Moisá SJ, Drackley JK, Moyes KM, Loor JJ. Bioinformatics analysis of microRNA and putative target genes in bovine mammary tissue infected with streptococcus uberis1. J Dairy Sci. 2012;95(11):6397–408.

    Article  CAS  PubMed  Google Scholar 

  68. Li R, Beaudoin F, Ammah AA, Bissonnette N, Benchaar C, Zhao X, Lei C, Ibeagha-Awemu EM. Deep sequencing shows microRNA involvement in bovine mammary gland adaptation to diets supplemented with linseed oil or safflower oil. BMC Genomics. 2015;16:884.

    Article  PubMed  PubMed Central  Google Scholar 

  69. Li R, Dudemaine P-L, Zhao X, Lei C, Ibeagha-Awemu EM. Comparative analysis of the miRNome of bovine milk fat, whey and cells. PLoS One. 2016;11(4):e0154129.

    Article  PubMed  PubMed Central  Google Scholar 

  70. Hudson NJ, Reverter A, Dalrymple BP. A differential wiring analysis of expression data correctly identifies the gene containing the causal mutation. PLoS Comput Biol. 2009;5(5):e1000382.

    Article  PubMed  PubMed Central  Google Scholar 

  71. Reverter A, Chan EKF. Combining partial correlation and an information theory approach to the reversed engineering of gene co-expression networks. Bioinformatics. 2008;24(21):2491–7.

    Article  CAS  PubMed  Google Scholar 

  72. Theocharidis A, van Dongen S, Enright AJ, Freeman TC. Network visualization and analysis of gene expression data using BioLayout Express3D. Nat Protocols. 2009;4(10):1535–50.

    Article  CAS  PubMed  Google Scholar 

  73. Zampieri M, Soranzo N, Altafini C. Discerning static and causal interactions in genome-wide reverse engineering problems. Bioinformatics. 2008;24(13):1510–5.

    Article  CAS  PubMed  Google Scholar 

  74. Zhang B, Horvath S. A general framework for weighted gene co-expression network analysis. Stat Appl Genet Mol Biol. 2005;4:1-37.

  75. Langfelder P, Horvath S. Eigengene networks for studying the relationships between co-expression modules. BMC Syst Biol. 2007;1(1):54.

    Article  PubMed  PubMed Central  Google Scholar 

  76. Maere S, Heymans K, Kuiper M. BiNGO: a Cytoscape plugin to assess overrepresentation of gene ontology categories in biological networks. Bioinformatics. 2005;21(16):3448–9.

    Article  CAS  PubMed  Google Scholar 

  77. Supek F, Bošnjak M, Škunca N, Šmuc T. REVIGO summarizes and visualizes long lists of gene ontology terms. PLoS One. 2011;6(7):e21800.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  78. Enright AJ, John B, Gaul U, Tuschl T, Sander C, Marks DS. MicroRNA targets in drosophila. Genome Biol. 2004;5(1):R1.

    Article  Google Scholar 

Download references

Acknowledgements

This study was conducted with funding from EMBRAPA (Macroprogram 1, 01/2005) and FAPESP (processes number 2015/00617-3 and 2015/24688-7). The authors would like to acknowledge the collaborative efforts among EMBRAPA, University of São Paulo, and Iowa State University.

Funding

This study was conducted with funding from EMBRAPA (Macroprogram 1, 01/2005) and FAPESP (processes number 2015/00617–3 and 2015/24688–7). LCAR and LLC were granted CNPq fellowships. The foment agencies had no role in the experimental design and samples collection, as well as data analysis and interpretation, and in writing of this manuscript.

Availability of data and materials

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

Author information

Authors and Affiliations

  1. Department of Animal Science, University of São Paulo, Piracicaba, SP, 13418-900, Brazil

    Gabriella B. Oliveira, Aline S. M. Cesar, Karina Y. Degaki, Mirele D. Poleti, Andrezza M. Felício & Luiz L. Coutinho

  2. Embrapa Southeast Livestock, São Carlos, SP, 13560-970, Brazil

    Luciana C. A. Regitano

  3. Department of Animal Science, Iowa State University, Ames, IA, 50011, USA

    James M. Reecy

  4. Department of Animal Science, University of Arkansas, Fayetteville, AR, 72701, USA

    James E. Koltes

Authors
  1. Gabriella B. Oliveira
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Luciana C. A. Regitano
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Aline S. M. Cesar
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. James M. Reecy
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Karina Y. Degaki
    View author publications

    You can also search for this author in PubMed Google Scholar

  6. Mirele D. Poleti
    View author publications

    You can also search for this author in PubMed Google Scholar

  7. Andrezza M. Felício
    View author publications

    You can also search for this author in PubMed Google Scholar

  8. James E. Koltes
    View author publications

    You can also search for this author in PubMed Google Scholar

  9. Luiz L. Coutinho
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

GBO, ASMC and MDP participated in the experimental design and RNA-sequencing. GBO, ASMC, AMF and JEK performed data analysis. GBO drafted the manuscript. GBO, ASMC, JMR, KD and LLC collaborated with interpretation and discussion of the results and LCAR, JMR and LLC provided the experimental environment and coordination. All authors have read and approved the final manuscript.

Corresponding author

Correspondence to Luiz L. Coutinho.

Ethics declarations

Ethics approval and consent to participate

All experimental procedures involving steers were approved by the Institutional Animal.

Care and Use Committee Guidelines from Brazilian Agricultural Research Corporation –.

EMBRAPA (Macroprograma 1, 01/2005) and sanctioned by the president Dr. Rui Machado.

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.

Full list of DE miRNAs target genes. The table contains the Fold change (FC) of each DE miRNA and the ENSEMBL and gene symbol of their target genes. (XLSX 101 kb)

Additional file 2: Table S2.

Full list of gene networks constructed by IPA using the miRNA’s target genes list. The table contains the gene symbol gene symbol of all genes inside network, the number of target genes inside each network and diseases and functions that network could participate. (XLSX 14 kb)

Additional file 3: Table S3.

List of most relevant canonical pathways enriched by IPA using the miRNA’s target genes list. The table contains the name of pathway and the target genes that are present in each pathway. (XLSX 19 kb)

Additional file 4: Figures S1, S2 and S3.

Gene networks most relevant for lipid metabolism that were constructed by IPA using the miRNA’s target genes list. Grey shapes represent target genes and the white shapes are other genes of the network that are not target genes. (PDF 915 kb)

Additional file 5: Figures S4, S5 and S6.

The most relevant pathways enriched by IPA using the miRNA’s target genes list. The shapes highlighted in purple represent the miRNAs target genes and the white shapes represent the other genes of the pathway that are not target genes. (PDF 3964 kb)

Additional file 6: Tables S4.

Full list of differential hubbing (DH) values of all genes and miRNAs used in PCIT analysis. The table contains the ENSEMBL and gene symbol of each gene. (XLSX 441 kb)

Additional file 7: Tables S5 and S6.

Full list of genes significantly correlated with differential hubbed (DH) genes: genes correlated with positively DH genes (Table S5) and negatively DH genes (Table S6). The table contains the ENSEMBL of genes and the correlation value of each gene. (XLSX 745 kb)

Additional file 8: Tables S7 and S8.

Gene ontology (GO) enrichment analysis by Cytoscape of genes significantly correlated with differentially hubbed (DH) genes: negative DH genes (Table S7) and positive DH genes (Table S8). The table contains gene symbol and ENSEMBL of each DH gene, gene ontology term ID, p-value and adjusted p-value of GO enrichments, the description of each GO term and the genes present in each enriched biological process. (XLSX 45 kb)

Additional file 9: Tables S9, S10 and S11.

Full list of PIF values and adjusted p-values of all genes in the dataset (Table S9), RIF1 score (Table S10) and RIF2 score (Table S11). (XLSX 1975 kb)

Additional file 10: Table S12.

Gene ontology (GO) enrichment analysis by Cytoscape of genes highly correlated with top PIF and RIF genes and miRNAs. The table contains gene symbol and ENSEMBL of each top gene, gene ontology term ID, p-value and adjusted p-value of GO enrichments, the description of each GO term and the genes present in each enriched biological process. (XLSX 76 kb)

Additional file 11: Tables S13 and S14.

Full list of all mRNA modules and module membership (MM) values for each gene in dataset: modules constructed in the High intramuscular fat (IMF) group (Table S13), modules constructed in the Low IMF group (Table S14). (XLSX 18066 kb)

Additional file 12: Tables S15 and S16.

Full list of all miRNA modules and module membership (MM) values for each miRNA in dataset: modules constructed in the High intramuscular fat (IMF) group (Table S15), modules constructed in the Low IMF group (Table S16). (XLSX 260 kb)

Additional file 13: Table S17.

List of miRNA modules negatively correlated with mRNA modules and the GO terms associated with them in High IMF (H) and Low IMF (L) groups. The table contains the module’s name, the correlation value, p-value of correlation, number of miRNAs and mRNAs in each module and the GO term enriched for each module. (XLSX 11 kb)

Additional file 14: Tables S18, S19, S20 and S21.

Gene ontology (GO) enrichment analysis by Cytoscape of genes in mRNA modules negatively correlated with miRNA modules: enrichment of genes from mRNA modules constructed in the High intramuscular fat (IMF) group (Table S18) and in the Low IMF group (Table S19). Gene ontology enrichment analysis by Cytoscape of target genes of hub miRNAs in modules negatively correlated with mRNA modules: enrichment of target genes of miRNA modules constructed in the High IMF group (Table S20) and in the Low IMF group (Table S21). The table contains the module’s name, gene ontology ID, p-value and adjusted p-value of GO, the description of them and the genes present in each biological process enriched. (XLSX 71 kb)

Additional file 15: Figure S7.

Correlation between the mRNA modules and IMF. A: modules identified in the High IMF GEBV (H) group. B: modules identified in the Low IMF GEBV (L) group. Modules with intense red color have a higher correlation (close to + 1) and those with intense green color have a more negative correlation (close to − 1). (PDF 173 kb)

Additional file 16: Figure S8.

Correlation between miRNA modules and intramuscular fat (IMF). A: modules identified in the High IMF GEBV (H) group. B: modules identified in the Low IMF GEBV (L) group. Modules with intense red color have a higher positive correlation (close to + 1) and those with an intense green color have a more negative correlation (close to − 1). (PDF 118 kb)

Additional file 17: Tables S22, S23, S24 and S25.

Gene ontology (GO) enrichment analysis by Cytoscape of genes in mRNA modules correlated with the intramuscular fat (IMF) trait: enrichment of genes from mRNA modules correlated with High IMF (Table S22) and Low IMF groups (Table S23). Gene ontology enrichment analysis by Cytoscape of target genes of hub miRNAs in modules correlated with the IMF trait: enrichment of target genes of miRNA modules correlated with High intramuscular fat (Table S24) and Low IMF groups (Table S25). The table contains module’s name, gene ontology ID, p-value and adjusted p-value of GO, the description of them and the genes present in each biological process enriched. (XLSX 98 kb)

Additional file 18: Figure S9.

Scale free topology model and mean connectivity of the mRNA networks based on the power value β. (A) H group. (B) L group. (PDF 92 kb)

Additional file 19: Figure S10.

Scale free topology model and mean connectivity of the miRNA networks based on the power value β. (A) H group. (B) L group. (PDF 85 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

Oliveira, G.B., Regitano, L.C.A., Cesar, A.S.M. et al. Integrative analysis of microRNAs and mRNAs revealed regulation of composition and metabolism in Nelore cattle. BMC Genomics 19, 126 (2018). https://doi.org/10.1186/s12864-018-4514-3

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12864-018-4514-3

Keywords

BMC Genomics

ISSN: 1471-2164

Contact us