Co-Expression Networks Reveal Potential Regulatory Roles of miRNAs in Fatty Acid Composition of Nelore Cattle

Fatty acid (FA) content affects the sensorial and nutritional value of meat and plays a significant role in biological processes such as adipogenesis and immune response. It is well known that, in beef, the main FAs associated with these biological processes are oleic acid (C18:1 cis9, OA) and conjugated linoleic acid (CLA-c9t11), which may have beneficial effects on metabolic diseases such as type 2 diabetes and obesity. Here, we performed differential expression and co-expression analyses, weighted gene co-expression network analysis (WGCNA) and partial correlation with information theory (PCIT), to uncover the complex interactions between miRNAs and mRNAs expressed in skeletal muscle associated with FA content. miRNA and mRNA expression data were obtained from skeletal muscle of Nelore cattle that had extreme genomic breeding values for OA and CLA. Insulin and MAPK signaling pathways were identified by WGCNA as central pathways associated with both of these fatty acids. Co-expression network analysis identified bta-miR-33a/b, bta-miR-100, bta-miR-204, bta-miR-365-5p, bta-miR-660, bta-miR-411a, bta-miR-136, bta-miR-30-5p, bta-miR-146b, bta-let-7a-5p, bta-let-7f, bta-let-7, bta-miR 339, bta-miR-10b, bta-miR 486, and the genes ACTA1 and ALDOA as potential regulators of fatty acid synthesis. This study provides evidence and insights into the molecular mechanisms and potential target genes involved in fatty acid content differences in Nelore beef cattle, revealing new candidate pathways of phenotype modulation that could positively benefit beef production and human consumption.

Introduction

Fatty acid (FA) content is an important trait that can influence the sensorial and nutritional value of beef and plays a significant role in molecular and physiological processes. Despite that the consumption of beef fatty acids is being associated with metabolic diseases such as type 2 diabetes and obesity effects, such as altered blood lipid and lipoprotein content (Wood et al., 2008), beef has beneficial effects on human health due to its high nutritional value, being also an important source of oleic acid (OA) (Laaksonen et al., 2005a). Likewise, conjugated linoleic acids (CLAs) could have a range of nutritional benefits in the diet (Valsta et al., 2005).

Fatty acid biosynthesis biological processes are complex and dependent on several regulatory mechanisms, such as post-transcriptional regulation of gene expression (Nakamura and Nara, 2003). In this sense, miRNAs have been shown to block the translation of target mRNAs and thereby post-transcriptionally regulate adipogenesis and several other biological processes involved in fatty acid metabolism in bovine (Guo et al., 2017).

Transcriptomic studies have shown that the expression of many miRNAs is species specific and tissue specific, indicating that miRNAs may have potential roles in organ and tissue development, metabolism, immune response (Lawless et al., 2014), milk production traits, and fertility (Fatima and Morris, 2013). In beef cattle, differences in the expression pattern of miRNAs have been identified in animals with different amounts of subcutaneous fat, which could indicate a potential regulatory role of these molecules in the development of adipose tissue (Jin et al., 2010) and fat metabolism (Romao et al., 2014). These studies have identified numerous miRNAs expressed in cattle, but the miRNA regulatory mechanisms that underlie these phenotypes are unclear.

A recent integrative analysis of miRNA–mRNA co-expression in this Nelore population revealed several genes and miRNAs as candidate regulators of intramuscular fat deposition. Glucose metabolism and inflammation processes were the main pathways found to influence intramuscular fat deposition in Nelore beef cattle (Oliveira et al., 2018). Furthermore, previous RNAseq studies including this population identified differences in the skeletal muscle transcriptome profile associated with extreme values of fatty acid content. Oleic acid and CLA-c9t11 content had significant effects on the expression level of genes related to oxidative phosphorylation, cell growth, survival, and migration (Cesar et al., 2016). However, there are no studies about miRNA–mRNA co-expression related to fatty acid composition in bovines.

The integration of previous transcriptomic studies (Cesar et al., 2016) with miRNA expression data information may help provide a better understanding of the molecular mechanisms involved in the variation of FA content and deposition. Therefore, the goal of this study was to perform an integrative miRNA and mRNA expression analysis in skeletal muscle of beef cattle to unravel novel regulatory networks and signaling pathways involved in fatty acid biosynthesis and composition.

Material and Methods

Ethics Statement

Experimental procedures were carried out in accordance with the relevant guidelines provided by the Institutional Animal Care and Use Committee Guidelines of the Embrapa Pecuária Sudeste—Protocol CEUA 01/2013. The Ethical Committee of the Embrapa Pecuária Sudeste (São Carlos, São Paulo, Brazil) approved all experimental protocols (approval code CEUA 01/2013) prior to the conduction of the study.

Phenotypic Data

Description of phenotypic data and genomic heritability for oleic acid (OA (C18:1 cis9) and conjugated linoleic acid (CLA-c9t11) content from Nelore steers were previously reported (Cesar et al., 2014), with genomic heritability mean values of 0.16 ± 0.11 and 0.04 ± 0.09, respectively.

Animals used in the RNAseq and miRNAseq analyses were selected for extreme values of OA and CLA based on the rank of their genomic estimated breeding values (GEBVs) in a larger population of 386 Nellore steers. The GEBVs were calculated by GenSel software using the Bayes B approach (Cesar et al., 2014). A total of 30 animals with extreme GEBVs for OA and CLA content were selected separately: the top 13 high (H-OA) and 15 low (L-OA) animals for OA content, and the top 15 high (H-CLA) and 15 low (L-CLA) animals for CLA content. The Nelore steers used in this study are exactly the same group of animals used in Cesar et al. (2016).

mRNA Expression Data

The processing and analysis of mRNA expression data from skeletal muscle from the same population of animals used in this study were previously described in Cesar et al. (2016), with the sample accession of mRNA expression data of PRJEB13188. In total, 16,710 genes were identified as expressed in skeletal muscle for OA, and 16,530 genes were expressed in skeletal muscle for CLA. These genes were used for co-expression analysis.

miRNA Expression Data

The processing and analysis of miRNA expression data from skeletal muscle and miRNA target predictions used in this study followed the same procedures previously described in De Oliveira et al. (2018). Therefore, we will give a brief description of the analyses carried out.

In brief, sequencing of miRNA cDNA libraries was conducted on a MiSeq (Illumina, San Diego, CA) with MiSeq Reagent Kit 50 cycles in the Laboratory Multiuser ESALQ in Piracicaba/SP/Brazil, according to the protocol described by Illumina. The FastQC tools (http://www.bioinformatics.babraham.ac.uk/projects/fastqc) and FASTX (http://hannonlab.cshl.edu/fastx-toolkit) were used to check the quality of reads, and reads were subjected to alignment to bovine genome reference UMD version 3.1 (Ensembl 84: Mar 2016) through the software miRDeep2 version 2.0.0.7 (Friedländer et al., 2008). The reads were then mapped to regions of the genome using the Bowtie tool (Langmead et al., 2009), built into miRDeep2 software.

Differentially Expressed miRNAs

Differentially expressed (DE) miRNAs were identified for OA and CLA content phenotypes from a total of 30 small RNA libraries derived from skeletal muscle (N = 13 H-OA, N = 15 L-OA and N = 15 H-CLA, N = 15 L-CLA) using DESeq2 software (Love et al., 2014). The Benjamini–Hochberg method (Benjamini and Hochberg, 1995) was used to control for the rate of false positive (FDR; 10%) due to the number of genes and miRNAs tested. We set an FDR threshold of 0.1 (i.e., 10% of false positives are expected) to correct for false positives, avoiding to lose information, as these are exploratory analyses that should indicate biological responses to be in the future verified.

miRNA Target Predictions and Functional Enrichment Analysis

The target genes of DE miRNAs from skeletal muscle were predicted with TargetScan (Agarwal et al., 2015) and miRanda (Betel et al., 2010) software. In order to predict the potential regulatory target transcripts, the target genes were filtered by skeletal muscle (Cesar et al., 2016) mRNA expression data previously analyzed on the same set of samples. Functional enrichment analysis of target genes was performed by WebGestalt (Wang et al., 2017) using Bos taurus and the overrepresentation enrichment analysis (ORA) as organism and method of interest.

miRNA and mRNA Co-Expression Approaches: Weighted Gene Co-Expression Network Analysis (WGCNA) and Partial Correlation With Information Theory (PCIT)

WGCNA

Co-expression networks were constructed by WGCNA (Langfelder and Horvath, 2008) v1.36 package in RStudio environment using miRNA (N = 13 H-OA, N = 15 L-OA and N = 15 H-CLA, N = 15 L-CLA) and mRNA (N = 13 H-OA, N = 15 L-OA and N = 15 H-CLA, N = 15 L-CLA) skeletal muscle expression data.

miRNA and mRNA networks were constructed separately for high (H) and low (L) fatty acid groups. miRNA network construction and module detection used the step-by-step signed network construction with a soft threshold of β = 6 (R ² > 0.90) and a minimum module size of 5. The same approach was adopted for mRNA signed network construction with a soft threshold β = 6 (R ² > 0.91) and a minimum module size of 30. 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 (Langfelder and Horvath, 2008; Betel et al., 2010). 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.

An integrative analysis was performed, in which miRNA module eigengenes (MEs) and mRNA MEs for high and low fatty acid groups were correlated with one another by calculating the Pearson correlations. miRNA and mRNA modules with a negative correlation and a p-value < 0.10 were selected for functional enrichment analysis. miRNA modules that were significantly correlated were then further explored to identify hub miRNAs. Hub miRNAs were selected based on the top five greatest module membership (MM) values.

In order to better understand the biological significance of the modules identified, the functional enrichment analysis of genes and miRNA target genes were performed by WebGestalt (Wang et al., 2017) web tool. The functional enrichment analysis used the list of target genes from hub miRNAs selected from miRNA modules that were negatively correlated with mRNA modules. Co-expression networks among hub miRNAs and the GO terms of the target genes were constructed in Cytoscape v.3.3.0 (Cline et al., 2007).

PCIT: Differential Hubbing (DH), Regulatory Impact Factor (RIF), and Phenotypic Impact Factor (PIF) Metrics

The gene list used for PCIT analyses included all miRNAs and mRNAs detected in our study, but only those with a direct and partial correlation greater than 0.90 were used for the differential hubbing (DH) analyses. The DH was computed by the difference of significant connections of mRNA and miRNA between the high and low OA and CLA groups. Functional enrichment analysis was performed on the top five differential hubbing genes, or for miRNAs based on the list of predicted target genes, to determine if these specific genes/miRNAs have biological relevance to fatty acid composition.

The regulatory impact factor (RIF) and phenotypic impact factor (PIF) scores were calculated as described in Reverter et al. (2010) to predict which transcripts were potential regulators of gene expression differences between the high and low OA and CLA groups. Regulatory impact factor 1 (RIF1) ranks genes as potential regulators of networks based largely on changes in correlations between two different states, while regulatory impact factor 2 (RIF2) ranks genes with more emphasis on how expression level changes between two different states (Reverter et al., 2010). Phenotypic impact factor (PIF) values were used to identify and rank genes based on the magnitude of gene expression and the difference in the expression of that gene between two treatments (Reverter et al., 2010). The RIF calculations presented here were modified from the original method, and the complete list of expressed mRNA or miRNA was tested as potential regulators, and only mRNAs or miRNAs with a significant partial correlation of 0.90 from PCIT were included in the RIF and PIF score estimates; as described in Cesar et al. (2016). PIF score estimates were ranked to select the top 10 PIF regulators of all genes in the dataset.

Results

Phenotypic and miRNA Expression Data

Genomic estimated breeding values (GEBVs) and the number of normalized mapped miRNA reads for Nelore steers genetically divergent for two FA content, oleic acid (C18:1 cis9) and conjugated linoleic acid (CLA-c9t11), are shown in Table 1. Sample sizes of extreme Nelore steers used for OA and CLA analyses were slightly different due to data availability. A Student’s t-test was previously performed by Cesar et al. (2016) to evaluate the mean differences between the high and low FA groups, and significant differences (p-value < 0.05) were observed for both OA and CLA content.

TABLE 1

Table 1 Genomic estimated breeding values (GEBVs) and number of normalized mapped miRNA reads for Nelore steers with divergent fatty acid content groups.

On average, 85% of miRNA reads for OA and CLA were mapped to the B. taurus UMD 3.1 genome assembly (Ensembl 84: Mar 2016). In total, 404 OA and 386 CLA mature miRNAs were detected by miRDeep2 software (Supplementary Table 1).

Differentially Expressed miRNAs and Target Gene Identification

We identified 137 and 131 unique mature miRNA sequences with non-zero expression levels according to DESeq criteria for samples in OA (Supplementary Table 2) and CLA (Supplementary Table 3) groups, respectively, which were used in differential expression and co-expression analyses. The miRNAs bta-miR-126-5p and bta-miR-2419-5p were upregulated in the L-OA and H-CLA groups, respectively (Table 2).

TABLE 2

Table 2 Differentially expressed miRNAs identified by miRDeep2 for Nelore steers with divergent fatty acid content groups and the number of predicted target genes.

Among the 3,619 bta-miR-126-5p target genes (Supplementary Table 4) identified in the bovine genome, 162 were previously identified as DE in this same population (Cesar et al., 2016). SCD (stearoyl-CoA desaturase), CDS2 (CDP-diacylglycerol synthase 2), FAR2 (fatty acyl CoA reductase 2), and NAB1 (NGFI-A binding protein 1) were included in this list of DE genes, well known to be related to biological processes associated with fatty acid composition. Regarding the 365 bta-miR-2419-5p target genes (Supplementary Table 4), 16 were also previously identified as DE (Cesar et al., 2016). Among these genes are included CAV3 (caveolin 3), JMJD1C (jumonji domain containing 1C), FOXO6 (forkhead box O6), and PRKAG2 (protein kinase AMP-activated non-catalytic subunit gamma 2). These bta-miR-126-5p DE target genes were also enriched for genes involved in biological processes associated with fatty acid composition.

miRNA and mRNA Co-Expression

Weighted Gene Co-Expression Network Analysis (WGCNA)

WGCNA was used to identify potential regulatory networks related to OA and CLA content in skeletal muscle. To this end, miRNA and mRNA expression data from animals with extreme GEBVs of OA and CLA were evaluated separately. After quality control, a total of 137 miRNAs and 16,710 mRNAs were analyzed for OA network construction, while a total of 131 miRNAs and 16,530 mRNAs were used for CLA network construction.

A total of 52 mRNA modules were identified in the H-OA group (Figure S1A), while in the L-OA group, 95 mRNA modules were identified (Figure S1B). Regarding miRNA, 12 and 9 modules were identified, respectively, in the H-OA (Figure S2A) and L-OA (Figure S2B) groups. For CLA, 52 and 28 mRNA modules in the H-CLA (Figure S3A) and L-CLA (Figure S3B) groups were identified, respectively. CLA miRNA network analysis identified six and five miRNA modules in the H-CLA (Figure S4A) and L-CLA (Figure S4B), respectively.

In order to investigate miRNA–mRNA interactions, mRNA module eigengenes (MEs), which represent the sum of gene expression profiles of each module, were correlated with miRNA MEs. The mRNA and miRNA modules with negative correlations and a nominal p-value < 0.10 were selected for further investigation. The focus on negative correlation between mRNA and miRNA modules was based on the fact that the predominant and canonical effect of miRNAs on gene expression is through mRNA downregulation, which would equate to a negative miRNA–mRNA expression correlation (Guo et al., 2010). In the H-OA group, four negative correlations between mRNA and miRNA MEs were identified, while in the L-OA group, eight negative correlations were observed (Table 3). In the H-CLA group, seven negative correlations between miRNA and mRNA MEs were observed, while in the L-CLA group, three negative correlations were observed (Table 3).

TABLE 3

Table 3 Signaling pathways of miRNA module eigengenes (MEs) negatively correlated with mRNA MEs for Nelore steers with divergent fatty acid content groups.

Hub miRNAs are the miRNAs with higher connectivity inside the module and are probably more informative (Filteau et al., 2013). In order to find hub miRNAs involved in the co-expression networks, we selected the top five miRNAs representing the greatest module membership (MM) values from miRNA modules negatively correlated with mRNA modules. Table 3 shows signaling pathways obtained from WebGestalt software based on enrichment of the genes from miRNA modules and the target genes for hub miRNAs.

In the H-OA group, genes from mRNA modules were not significantly enriched for any Gene Ontology terms, while miRNA target genes from magenta and turquoise modules were significantly enriched for insulin resistance (Supplementary Table 5). miRNA target genes from black module were enriched for the MAPK signaling pathway (Supplementary Table 5). In the L-OA group, genes from the dark olive green module were significantly enriched for fatty acid degradation processes, while miRNA target genes from turquoise, pink, and blue modules were enriched for AMPK, insulin signaling pathway, and proteoglycans in cancer, respectively (Supplementary Table 5).

In the H-CLA group, genes from the mRNA modules brown and plum were significantly enriched for insulin resistance and steroid biosynthesis, respectively, while miRNA target genes from red, blue, and black modules were enriched for insulin resistance and insulin signaling pathway (Supplementary Table 6). In the L-CLA group, miRNA target genes from blue and red modules were enriched for insulin and MAPK signaling pathway, respectively (Supplementary Table 6).

Figures 1 and 2 show co-expression networks for hub miRNAs enriched for signaling pathways related to FA composition in H-OA (Figure 1) and L-OA groups (Figure 2) and in H-CLA (Figure 3) and L-CLA (Figure 4) groups in Nelore cattle.

FIGURE 1

Figure 1 Co-expression networks of H-OA group from Nelore cattle. Colored diamonds represent the top five hub miRNAs within each module, and colored rectangles represent the signaling pathways associated (p-value ≤ 0.10) with hub miRNA target genes.

FIGURE 2

Figure 2 Co-expression networks of L-OA group from Nelore cattle. Colored diamonds represent the top five hub miRNAs within each module, and colored rectangles represent the signaling pathways associated (p-value ≤ 0.10) with hub miRNA target genes.

FIGURE 3

Figure 3 Co-expression networks of H-CLA group from Nelore cattle. Colored diamonds represent the top five hub miRNAs within each module, and colored rectangles represent the signaling pathways associated (p-value ≤ 0.10) with hub miRNA target genes.

FIGURE 4

Figure 4 Co-expression networks of L-CLA group from Nelore cattle. Colored diamonds represent the top five hub miRNAs within each module, and colored rectangles represent the signaling pathways associated (p-value ≤ 0.10) with hub miRNA target genes.

Partial Correlation With Information Theory (PCIT)

PCIT, regulatory impact factor (RIF), and phenotypic impact factor (PIF) were used to score genes as potential regulators of signaling pathways between high and low FA groups. Differential hubbing (DH) represents the number of significant partial correlations that a gene has between two phenotypic states (Hudson et al., 2009). DH values of all genes and miRNAs used in PCIT analysis are in Supplementary Table 7. Table 4 shows the top five negative and positive extreme hubbing genes when comparing high and low OA and CLA groups. Functional enrichment analysis was performed on the top five DH genes to determine if they have biological relevance to fatty acid composition. The top 10 negatively and positively hubbed genes for OA and CLA content and their associated GO terms are shown in Supplemental Table 8.

TABLE 4

Table 4 Top negative and positive extreme hubbing genes groups for oleic acid (OA) and conjugated linoleic acid (CLA) content.

The bta-mir-339a and FRAT1 are among the top negatively hubbed genes for OA, associated with glycogen synthase kinase-3 binding protein (IPR008014) and Wnt signaling pathway (bta04310). GAL3ST3 and ATP6V0E1 are among the top positively hubbed genes, associated with the glycolipid biosynthetic process (GO: 0009247) and generation of metabolites and energy, respectively. KAT5 is among the top negatively hubbed genes for CLA, associated with proteasome-mediated ubiquitin-dependent protein catabolic process (GO: 0043161), while TMEM115 is associated with protein glycosylation (GO: 0006486). PSMG1 is among the top positively hubbed genes and associated with proteasome assembly (GO: 0043248).

For OA, a total of 14,900 transcripts had negative RIF1 values, whereas 1,948 transcripts had positive values. For RIF2, 16,684 transcripts had negative values, while 127 transcripts had positive values (Supplementary Table 9). For CLA, a total of 14,171 transcripts had negative RIF1 values, while 1,891 transcripts had positive values. For RIF2, 108 transcripts had negative values, while 16,512 transcripts had positive values (Supplementary Table 10).

Table 5 shows the top five negative and positive genes identified by RIF1 and RIF2 score by contrasting H and L groups for OA content. The RIF1 analysis identified putative regulators for OA content GO terms associated with muscling, e.g., muscle contraction and actin filament organization (TPM1, TPM2, and MYL1), and with fatness, e.g., glycolytic process and ATP biosynthetic process (ALDOA). The RIF2 analysis identified bta-mir-10b as a negative putative regulator, as well as the genes ACTN2 and TNNT1, associated with skeletal muscle contraction GO terms. The same group of genes identified as positive RIF1 regulators is among the top positive RIF2 regulators. Functional enrichment analysis was performed using DAVID software on the top five genes to determine their biological relevance on fatty acid composition (Supplementary Table 9).

TABLE 5

Table 5 Top negative and positive genes identified by regulatory impact factor 1 (RIF1) and regulatory impact factor 2 (RIF2) score for oleic acid (OA) content.

Table 6 shows the top negative and positive RIF1 and RIF2 genes when high and low groups for CLA content were contrasted. RIF1 analysis identified putative negative regulators for CLA content genes related to regulation of gene expression (GNRH1) and muscling, e.g., actin filament organization (TRPV4). The genes ALDOA, CKM, and TPM were identified as positive RIF1 regulators and were also identified as negative RIF2 regulators. RIF2 analysis identified positive putative regulator genes enriched for GO terms associated with muscling (MYH1, DES, and PDLIM3). Functional enrichment analysis on the top five genes is presented in Supplementary Table 10.

TABLE 6

Table 6 Top negative and positive genes identified by regulatory impact factor 1 (RIF1) and regulatory impact factor 1 (RIF2) score for conjugated linoleic acid (CLA) content.

For OA and CLA, PIF analysis identified 13,480 and 15,531 transcripts, respectively (adjusted p-value < 0.05), which could significantly impact fatty acid composition. Table 7 lists the top 10 regulators ranked by PIF analysis for OA and CLA content. Functional enrichment analysis to determine the biological relevance on fatty acid composition of the top 10 genes is presented in Supplementary Table 11.

TABLE 7

Table 7 Top 10 regulators identified by phenotypic impact factor (PIF) analysis with FDR adjusted p-value for oleic acid (OA) and conjugated linoleic acid (CLA) content.

Interestingly, among the top 10 regulators identified by PIF analyses for OA and CLA content is the same group of genes (ACTA1, TTN, MYH1, ALDOA, CKM, and NEB) with Gene Ontology (GO) terms associated with muscling (ACTA1, MYH, and MYLPF) and fatness and fiber type (ALDOA).

Discussion

The purpose of this study was to investigate the complex interactions between miRNAs and mRNAs in bovine skeletal muscle associated with variation in oleic acid (OA) and conjugated linoleic acid (CLA-c9t11) content. This was accomplished by performing differential miRNA expression analysis and two different gene co-expression approaches, weighted gene co-expression network analysis (WGCNA) and partial correlation with information theory (PCIT). Furthermore, we sought to identify the key drivers of gene expression networks. Analyses were focused on these two fatty acids due to their importance in many biological processes and beneficial effects on metabolic diseases and human health (Laaksonen et al., 2005b). Also, a significant number of differentially expressed genes in response to OA (1134) and CLA (872) content were identified in this same population previously (Cesar et al., 2016). OA is a monounsaturated fatty acid present in membrane phospholipids, triglycerides, and cholesterol and is associated with protection against heart disease (Laaksonen et al., 2005b), while CLA antidiabetic effects are mediated via anti-inflammatory processes in white adipose tissue (Moloney et al., 2007).

Integration of miRNA and mRNA co-expression from next-generation sequencing data from the same individuals revealed the possible regulatory roles of miRNAs on fatty acid composition. WGCNA allowed the identification of miRNA and mRNA modules that were negatively correlated with each other, which may indicate that FA composition is modulated by specific miRNA–mRNA interactions. However, only one mRNA module in the L-OA group and two mRNA modules in the H-CLA group presented functional enrichment for fatty acid degradation, insulin resistance, and steroid biosynthesis. Therefore, the subsequent analyses were focused on the identification of key miRNAs that may be involved in co-expression networks and thereby in the regulation of fatty acid composition.

From differential expression analysis, the SCD gene, identified as a putative target gene of bta-miR-126-5p, was found to be upregulated in the H-OA content group in a previous skeletal muscle RNAseq study (Cesar et al., 2016). SCD is as a key enzyme in de novo lipogenesis (Scaglia and Igal, 2005), and its upregulation has been associated with deposition of unsaturated FAs. Thus, the upregulated bta-miR-126-5 in the L-OA group observed herein may explain the reduced SCD mRNA levels observed before (Cesar et al., 2016). These results suggest that SCD expression level is related to OA content and may suggest a new role for bta-miR-126-5p in the regulation of SCD. It is relevant to note that only negative miRNA–mRNA correlations were considered herein. However, it is important to emphasize that this was an exploratory study and should be complemented by other in vitro and in vivo analyses to better discriminate the mechanisms of miRNA regulation of gene expression.

Functional enrichment analysis indicated a relationship between insulin, insulin resistance, adipocytokine signaling pathway, and non-alcoholic fatty liver disease for PRKAG2 gene, a target gene of bta-miR-2419-5p. The major effects of insulin on muscle and adipose tissue are related to carbohydrate, lipid, and protein metabolism (Dimitriadis et al., 2011). Insulin can decrease the rate of lipolysis, stimulate fatty acid synthesis, increase the uptake of triglycerides from blood into muscle and adipose tissue, and decrease the rate of fatty acid oxidation in muscle and liver (Dimitriadis et al., 2011). Metabolic diseases such as obesity and coronary disorders can be the consequence of insulin resistance, i.e., the inability of insulin to drive glucose into muscle and other tissues, which can be caused by excessive body fat deposition (Dimitriadis et al., 2011). Assuming that bta-miR-2419-5p expression level regulates CLA content and that this miRNA can regulate PRKAG2 expression, it is possible that bta-miR-2419-5p can regulate insulin expression to ultimately respond to CLA content. The transcription factors (TFs) FOXO1, FOXO3, and FOXO4 are also targets of bta-miR-2419-5p. Therefore, they may also be regulated by CLA content. FOXO1 expression has been previously associated with CLA content in this population (Cesar et al., 2016).

Co-expression network analysis integrating miRNA and mRNA expression data revealed two miRNA modules (magenta and turquoise) in the H-OA group whose potential target genes were significantly enriched for the GO term insulin resistance. As discussed previously, it is known that insulin resistance can be caused by excessive fat deposition (Ortega et al., 2013). The miRNA magenta module grouped both bta-miR-181a and bta-miR-33a/b. The bta-miR-181a is from the same family as bta-miR-181b, which has been reported to regulate the biosynthesis of bovine milk by targeting ACSL1, an important enzyme of milk lipid synthesis (Lian et al., 2016). Furthermore, bta-miR-33a/b has been reported to contribute to the regulation of fatty acid metabolism and the insulin signaling pathway (Dávalos et al., 2011), which indicates that the insulin signaling pathway may be impacted in the H-OA group.

The miRNA turquoise module contained both bta-miR-146b and bta-miR-26b. In a recent integrative analysis of miRNA–mRNA expression related to intramuscular fat (IMF) deposition in animals from this population, the bta-miR-146b was found to be downregulated in individuals with high IMF content, while the bta-miR-26b was identified by PCIT as a candidate regulatory gene that negatively regulates IMF deposition (Oliveira et al., 2018). IMF represents the amount of fat accumulated between muscle fibers or within muscle cells, and it is a determinant factor that affects meat quality (Wood et al., 2008). Despite the fact that in the present study IMF deposition between high and low groups was not statistically different (Cesar et al., 2016), we could not ignore the fact that OA, CLA, and other FAs comprise IMF. It has been reported that as intramuscular lipid content accumulates, there is a concomitant elevation in the concentration of oleic acid (Smith et al., 2006). Other hub miRNAs such as bta-miR-196a and bta-miR-30f (from the same family of hubs bta-miR-30d and bta-miR-30a-5p) have been reported to have a higher expression level in cattle that have higher amounts of IMF (Guo et al., 2017). Therefore, taking together the results from these and previous studies, we can highlight the role of bta-miR-146b, bta-miR-26b, bta-miR-30d, and bta-miR-196a in regulating muscle fatty acid composition.

Lipid accumulation can activate the immune system and inflammatory pathways due to the secretion of proinflammatory molecules by adipocytes (Lau, 2005). The mitogen-activated protein kinase (MAPK) cascade, enriched from the miRNA black module, is highly conserved and involved in various cellular functions, including responses to proinflammatory stimuli (Soares-Silva et al., 2016). Networks and DE genes related to immune system and inflammatory response were previously associated with high amounts of IMF in this population (Cesar et al., 2015; Oliveira et al., 2018). Among the hub miRNAs from the black module was the bta-miR100, which was previously found to be downregulated in animals with high IMF (Oliveira et al., 2018).

In the L-OA group, the target genes of two miRNA modules (pink and green) were enriched for GO terms associated with the insulin signaling pathway. Insulin is a hormone with a direct effect on lipid metabolism (Dimitriadis et al., 2011), which could explain why modules present in the high and low fatty acid groups are associated with insulin-related terms. However, we can observe an overrepresentation of the insulin signaling pathway in the low fatty acid group, while in the high fatty acid group, we can observe an overrepresentation of the insulin resistance pathway. As discussed before, insulin resistance has been linked to excessive body fat deposition and obesity (Dimitriadis et al., 2011), supporting our findings.

Hub miRNAs from the pink module (bta-miR-204 and bta-miR-365-5p) and from the green module (bta-miR-660) have been previously associated with adipose tissue in cattle (Gu et al., 2007). Other hub miRNAs from the green module (bta-miR-411a and bta-miR-136) were expressed at a higher level in Wagyu compared with Holstein cattle (Guo et al., 2017). Wagyu cattle accumulate large amounts of marbling and specifically monounsaturated fatty acids, of which oleic acid is primarily responsible for the soft fat (Smith et al., 2009). Taken together, these results may indicate a possible role of these miRNAs in post-transcriptional regulation of OA content.

Target genes of the miRNA turquoise module were enriched for GO terms associated with the AMPK signaling pathway. AMPK is a basic regulator of cellular and body energy metabolism and may enhance activity of mitochondrial proteins involved in oxidative metabolism (Thomson et al., 2008). Cesar et al. (2016) reported that several canonical pathways of oxidative phosphorylation were upregulated in animals with H-OA content. Therefore, the enrichment of target genes associated with the AMPK signaling pathway in L-OA might indicate a post-transcriptional regulation of this pathway resulting in downregulated oxidative metabolism in these animals, complementing the findings by Cesar et al. (2016). Still on the turquoise module, the hub bta-miR-146b, which is in the same miRNA family as bta-miR-146a, has been correlated with target genes that are functionally enriched for GO terms associated with fatty acid oxidation (Oliveira et al., 2018).

Target genes of the miRNA blue module were enriched for GO terms associated with proteoglycans in cancer pathways. Proteoglycans have been shown to be key macromolecules that contribute to biology of various types of cancer (Iozzo and Sanderson, 2011). Previous studies have reported an important contribution of OA intake to human health, with protective effects against cancer development (Schwartz et al., 2008). In this sense, genes related to cancer were found to be upregulated in animals with low CLA content (Cesar et al., 2016). These findings may be evidence of miRNA modulation of a carcinogenic pathway. In the miRNA blue module as well, a hub miRNA, bta-miR-21-5p, has been reported to be an important regulator of bovine mammary lipogenesis and metabolism (Li et al., 2015). As in meat, fatty acid composition can influence the nutritional quality of milk and milk fat (Soyeurt et al., 2008). From miRNA yellow module, target genes were enriched for Wnt signaling pathway. Wnt is a member of signal transduction pathways, which regulates crucial aspects of cell fate determination (Komiya and Habas, 2008). It has been reported that the knockdown of a key enzyme in fatty acid synthesis (FASN) could attenuate the Wnt signaling pathway via downregulation of specific genes (Wang et al., 2016).

Two miRNA modules (red and blue) identified in the H-CLA group were enriched for GO terms associated with insulin signaling pathway and insulin resistance. The hub bta-miR-30b-5p from red module has been reported to regulate muscle cell differentiation (Zhang et al., 2016), while hub bta-miR-10b from blue module and bta-miR-146b from grey module have been reported to have a higher expression in mammary tissue (Wicik et al., 2016). Moreover, bta-miR-146b was associated with the GO terms of Wnt and inflammatory pathways. The canonical Wnt signaling pathway was also enriched in the L-OA group.

In the L-CLA group, target genes in the miRNA blue module were enriched for the insulin signaling pathway and target genes of the miRNA red module for MAPK signaling pathway GO terms. Besides being involved with inflammatory processes, the MAPK signaling pathway is involved in the activation of PPARα (peroxisome proliferator-activated receptor) by adiponectin, stimulating fatty acid oxidation in muscle cells (Myeong et al., 2006). Adiponectin is an adipocytokine secreted by adipocytes with its beneficial effects on insulin resistance and metabolic disorders (Myeong et al., 2006). Interestingly, hub miRNAs from the red module such as bta-let-7a-5p, bta-let-7f, and bta-let-7e have been associated with insulin-like growth factor receptor signaling pathway (Wicik et al., 2016).

To better understand the biological processes that influence muscle FA composition, PCIT analysis was conducted. In this analysis, we identified the potential regulators that could be involved in the gene expression changes in skeletal muscle due to OA and CLA content. Negative and positive regulators are defined based on the number of significant partial correlations that a gene has between two states (Reverter and Chan, 2008).

Differential hubbing (DH) analysis identified bta-mir-339a as one of the top five negative regulators, with more connections in the L-OA group. The target genes of bta-mir-339a were enriched for 10 significant pathways, including the MAPK signaling pathway, which was also identified by WGCNA. Furthermore, bta-mir 339 has been reported to be expressed at a higher level in bovine adipose tissue than in other tissues (Gu et al., 2007). DH analysis also identified the FRAT1 gene as a potential negative regulator, which is associated with the Wnt signaling pathway. As discussed previously, the Wnt signaling pathway would be affected by the knockdown of FASN, resulting in lower fatty acid synthesis, which complements our findings from WGCNA in the L-OA group. Among the top 10 positive differentially hubbed genes is the ATP6V0E gene, from the same family of ATP6V1D gene, which has been reported as a factor mediating hepatic steatosis (Nakadera et al., 2016), a metabolic syndrome frequently associated with obesity and diabetes. Furthermore, the GAL3ST3 gene has been associated with lipid biosynthetic processes (Suzuki et al., 2001).

For CLA content, DH analysis identified KAT5 as extreme negative and PSMG1 as extreme positive hubbed genes associated with GO terms for proteasome pathways. The proteasome is a large protein complex responsible for degradation of intracellular proteins (Tanaka, 2009), and proteasome dysfunction has been associated with oxidative stress and insulin sensitivity in human obesity (Díaz-Ruiz et al., 2015). Gonçalves et al. (2018) otherwise concluded that the proteasome pathway may be a potential regulator of beef tenderness in this population. The total lipid content of muscle has a recognized role in beef tenderness, and the concentration of fatty acids is positively correlated with the palatability of beef (Wood et al., 2008). The proteasome pathway was also previously associated with OA content in this population (Cesar et al., 2016), supporting our findings of genes related to proteasome pathways as potential regulators.

Bta-miR-10b was identified as a CLA candidate regulator by both RIF2 and PIF analyses, being also pointed as a hub miRNA by WGCNA. Although its target genes were not enriched for any pathways specifically related to fatty acid metabolism, this miRNA has been previously correlated with backfat thickness and adipose tissue in cattle (Gu et al., 2007; Jin et al., 2010). Bta-mir-486, which was identified by PIF analysis as a top regulator for CLA content, has been associated with skeletal muscle growth (Jing et al., 2015) and was recently found to be downregulated in feed efficient animals of this same cattle population (De Oliveira et al., 2018).

ACTA1 is the gene with the highest PIF rank, which indicates that it may be the most important gene related to both OA and CLA variations in this population. ACTA1 expression is specific to muscle fibers, with an essential role in muscle contraction and cell morphology (Pollard and Cooper, 1986). ALDOA gene was the fourth and second major regulator identified by RIF and PIF analyses for OA and CLA content, respectively. This gene is involved in adipogenic differentiation, which is critical for intramuscular fat deposition and meat quality (Li et al., 2016). The same group of genes (TPM2, CKM, TPM1, and MYL1), including ALDOA, was identified as positive RIF1 and RIF2 regulators for OA content, indicating the relevance of these genes in fatty acid metabolism. Moreover, Oliveira et al. (2018) also identified the ALDOA gene as a putative regulatory for the differences in IMF deposition, which, taken together, may indicate that ALDOA is an important gene regulator of fatty acid deposition and composition in Nelore cattle.

In this integrative analysis, insulin resistance, insulin, and MAPK signaling pathways were overrepresented in high and low fatty acid groups. These signaling pathways have been linked to adipocyte differentiation and lipogenesis in cattle (Guo et al., 2017). Based on the literature and our results, insulin and inflammatory processes are influencing OA and CLA composition in Nelore cattle. This study also indicates that hub miRNAs like bta-miR-33a/b, bta-miR-100, bta-miR-204, bta-miR-365-5p, bta-miR-660, bta-miR-411a, bta-miR-136, bta-miR-30-5p, bta-miR-146b, bta-let-7a-5p bta-let-7f, and bta-let-7e are involved with these biological processes. Among the results pointed out by both RIF and PIF analyses, the bta-mir 339, bta-mir-10b, bta-miR 486, and genes ACTA1 and ALDOA are the most relevant regulators for muscle fatty acid composition in Nelore cattle.

Fat and fatty acids, whether in adipose tissue or muscle, contribute to various aspects of meat quality and are central to the nutritional value of meat (Wood et al., 2008). Furthermore, they can have beneficial effects on human health. OA consumption is associated with low levels of low-density lipoprotein (LDL) or “bad cholesterol,” which in turn may reduce atherosclerosis risk and diabetes occurrence. Further, OA consumption could increase levels of high-density lipoprotein (HDL) in blood (Smith et al., 2009), whereas CLA consumption may contribute to reduced body fat, cardiovascular diseases, and cancer and can modulate inflammatory responses (Dilzer and Park, 2012).

Conclusions

In the present study, signaling pathways, miRNAs, and gene regulators related to fatty acid composition in Nelore cattle were identified by miRNA expression and gene co-expression network approaches. Although some of these potential regulators have been previously linked to fatty acid composition, the complex miRNA–mRNA regulatory network has never been reported so far. This study improves our understanding of the molecular mechanisms controlling intramuscular muscle fat composition in bovines, revealing new candidate networks regulating OA and CLA phenotypes, which could positively benefit beef production.

Ethics Statement

Experimental procedures were carried out in accordance with the relevant guidelines provided by the Institutional Animal Care and Use Committee Guidelines of the Embrapa Pecuria Sudeste—Protocol CEUA 01/2013. The Ethical Committee of the Embrapa Pecuria Sudeste (Sao Carlos, Sao Paulo, Brazil) approved all experimental protocols (approval code CEUA 01/2013) prior to the conduction of the study.

Author Contributions

PO, LC, and LR conceived and designed the experiment; PO, LC, AC, GM, AZ, and LC performed the experiments; PO, AC, WD, and MS performed analysis; PO, AC, WD, BA, JK, JR, and LR interpreted the results; PO, AC, WD, MS, BA, JK, JR, and LR drafted and revised the manuscript. All authors read and approved the final manuscript.

Funding

This study was conducted with funding from EMBRAPA, São Paulo Research Foundation (FAPESP, grant number: 2012/23638-8), and scholarship to PO (grant numbers: 2014/22235-1 and 2016/03291-4), the National Council for Scientific and Technological Development (CNPq, grant numbers: 449172/2014-7, 303754/2016-8, 444374/2014-0, and 309004/2016-0), and fellowships to LR, LC, and GM.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We thank the São Paulo Research Foundation (FAPESP) and the National Council for Scientific and Technological Development (CNPq) for providing financial support. We also thank the Iowa State University for accepting the first author as a visiting scholar. We thank Dr. Dante P. D. Lanna, Dr. Michele L. Nascimento, and Dr. Amália S. Chaves for monitoring the feedlots.

Abbreviations

CLA, conjugated linoleic acid; DH, differential hubbing; FA, fatty acid; GEBV, genomic estimated breeding value: IMF, intramuscular fat content: ME, module eigengene: MM, module membership: OA, oleic acid: PCIT, partial correlation with information theory: PIF, phenotypic impact factor: RIF1, regulatory impact factor 1: RIF2, regulatory impact factor 2: WGCNA, weighted gene co-expression network analysis:

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2019.00651/full#supplementary-material

References

Agarwal, V., Bell, G. W., Nam, J. W., Bartel, D. P. (2015). Predicting effective microRNA target sites in mammalian mRNAs. Elife. 4, e05005. doi: 10.7554/eLife.05005

CrossRef Full Text | Google Scholar

Benjamini, Y., Hochberg, Y. (1995). Controlling the false discovery rate: a practical and powerful approach to multiple testing. J. R. Stat. Soc. Series B Stat. 57 (1), 289–300. doi: 10.1111/j.2517-6161.1995.tb02031.x

CrossRef Full Text | Google Scholar

Betel, D., Koppal, A., Agius, P., Sander, C., Leslie, C. (2010). Comprehensive modeling of microRNA targets predicts functional non-conserved and non-canonical sites. Genome Biol. 11 (8), R90. doi: 10.1186/gb-2010-11-8-r90

PubMed Abstract | CrossRef Full Text | Google Scholar

Cesar, A. S. M., Regitano, L. C. A., Koltes, J. E., Fritz-Waters, E. R., Lanna, D. P. D., Gasparin, G., et al. (2015). Putative regulatory factors associated with intramuscular fat content. PLoS One. 10 (6), e0128350. doi: 10.1371/journal.pone.0128350

PubMed Abstract | CrossRef Full Text | Google Scholar

Cesar, A. S. M., Regitano, L. C. A., Poleti, M. D., Andrade, S. C. S., Tizioto, P. C., Oliveira, P. S. N., et al. (2016). Differences in the skeletal muscle transcriptome profile associated with extreme values of fatty acids content. BMC Genomics. 17 (1), 961. doi: 10.1186/s12864-016-3306-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Cesar, A. S., Regitano, L. C., Tullio, R. R., Lanna, D. P., Nassu, R. T., Mudado, M., et al. (2014). Genome-wide association study for intramuscular fat deposition and composition in Nellore cattle. BMC Genet. 15, 39. doi: 10.1186/1471-2156-15-39

PubMed Abstract | CrossRef Full Text | Google Scholar

Cline, M. S., Smoot, M., Cerami, E., Kuchinsky, A., Landys, N., Workman, C., et al. (2007). Integration of biological networks and gene expression data using Cytoscape. Nat. Protoc. 2 (10), 2366–2382. doi: 10.1038/nprot.2007.324

PubMed Abstract | CrossRef Full Text | Google Scholar

Dávalos, A., Goedeke, L., Smibert, P., Ramírez, C. M., Warrier, N. P., Andreo, U., et al. (2011). miR-33a/b contribute to the regulation of fatty acid metabolism and insulin signaling. Proc. Natl. Acad. Sci. U. S. A. 108, 9232–9237. doi: 10.1073/pnas.1102281108

PubMed Abstract | CrossRef Full Text | Google Scholar

De Oliveira, P. S. N., Coutinho, L. L., Tizioto, P. C., Cesar, A. S. M., de Oliveira, G. B., Diniz, W. J., et al., et al. (2018). An integrative transcriptome analysis indicates regulatory mRNA–miRNA networks for residual feed intake in Nelore cattle. Sci. Rep. 8 (1), 17072. doi: 10.1038/s41598-018-35315-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Díaz-Ruiz, A., Guzmán-Ruiz, R., Moreno, N. R., García-Rios, A., Delgado-Casado, N., Membrives, A., et al. (2015). Proteasome dysfunction associated to oxidative stress and proteotoxicity in adipocytes compromises insulin sensitivity in human obesity. Antioxid. Redox Signal. 23 (7), 597–612. doi: 10.1089/ars.2014.5939

PubMed Abstract | CrossRef Full Text | Google Scholar

Dilzer, A., Park, Y. (2012). Implication of conjugated linoleic acid (CLA) in human health. Crit. Rev. Food Sci. Nutr. 52 (6), 488-513. doi: 10.1080/10408398.2010.501409

PubMed Abstract | CrossRef Full Text | Google Scholar

Dimitriadis, G., Mitron, P., Lambadiari, V., Maratou, E., Raptis, S. A. (2011). Insulin effects in muscle and adipose tissue. Diabetes Res. Clin. Pract. Suppl 1, S52–9. doi: 10.1016/S0168-8227(11)70014-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Fatima, A., Morris, D. G. (2013). MicroRNAs in domestic livestock. Physiol. Genomics. 45 (16), 685–696. doi: 10.1152/physiolgenomics.00009.2013

PubMed Abstract | CrossRef Full Text | Google Scholar

Filteau, M., Pavey, S. A., St-Cyr, J., Bernatchez, L. (2013). Gene coexpression networks reveal key drivers of phenotypic divergence in lake whitefish. Mol. Biol. Evol. 30 (6), 1384–1396. doi: 10.1093/molbev/mst053

PubMed Abstract | CrossRef Full Text | Google Scholar

Friedländer, M. R., Chen, W., Adamidi, C., Maaskola, J., Einspanier, R., Knespel, S., et al. (2008). Discovering microRNAs from deep sequencing data using miRDeep. Nat. Biotechnol. 26 (4), 407–415. doi: 10.1038/nbt1394

PubMed Abstract | CrossRef Full Text | Google Scholar

Gonçalves, T. M., de Almeida Regitano, L. C., Koltes, J. E., Cesar, A. S. M., da Silva Andrade, S. C., Mourão, G. B., et al. (2018). Gene co-expression analysis indicates potential pathways and regulators of beef tenderness in nellore cattle. Front. Genet. 9, 441. doi: 10.3389/fgene.2018.00441

PubMed Abstract | CrossRef Full Text | Google Scholar

Gu, Z., Eleswarapu, S., Jiang, H. (2007). Identification and characterization of microRNAs from the bovine adipose tissue and mammary gland. FEBS Lett. 581 (5), 981–988. doi: 10.1016/j.febslet.2007.01.081

PubMed Abstract | CrossRef Full Text | Google Scholar

Guo, H., Ingolia, N. T., Weissman, J. S., Bartel, D. P. (2010). Mammalian microRNAs predominantly act to decrease target mRNA levels. Nature. 466 (7308), 835–840. doi: 10.1038/nature09267

PubMed Abstract | CrossRef Full Text | Google Scholar

Guo, Y., Zhang, X., Huang, W., Miao, X. (2017). Identification and characterization of differentially expressed miRNAs in subcutaneous adipose between Wagyu and Holstein cattle. Sci. Rep. 7, 44026. doi: 10.1038/srep44026

PubMed Abstract | CrossRef Full Text | Google Scholar

Hudson, N. J., Reverter, A., Dalrymple, B. P. (2009). A differential wiring analysis of expression data correctly identifies the gene containing the causal mutation. PLoS Comput. Biol. 5 (5), e1000382. doi: 10.1371/journal.pcbi.1000382

PubMed Abstract | CrossRef Full Text | Google Scholar

Iozzo, R. V., Sanderson, R. D. (2011). Proteoglycans in cancer biology, tumour microenvironment and angiogenesis. J. Cell. Mol. Med. 15 (5), 1013–1031. doi: 10.1111/j.1582-4934.2010.01236.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Jin, W., Dodson, M. V., Moore, S. S., Basarab, J. A., Guan, L. L. (2010). Characterization of microRNA expression in bovine adipose tissues: a potential regulatory mechanism of subcutaneous adipose tissue development. BMC Mol. Biol. 11, 29. doi: 10.1186/1471-2199-11-29

PubMed Abstract | CrossRef Full Text | Google Scholar

Jing, L., Hou, Y., Wu, H., Miao, Y., Li, X., Cao, J., et al. (2015). Transcriptome analysis of mRNA and miRNA in skeletal muscle indicates an important network for differential residual feed intake in pigs. Sci. Rep. 5, 11953. doi: 10.1038/srep11953

PubMed Abstract | CrossRef Full Text | Google Scholar

Komiya, Y., Habas, R. (2008). Wnt signal transduction pathways. Organogenesis. 4 (2), 68–75. doi: 10.4161/org.4.2.5851

PubMed Abstract | CrossRef Full Text | Google Scholar

Laaksonen, D. E., Lindström, J., Lakka, T. A., Eriksson, J. G., Niskanen, L., Wikström, K., et al. (2005a). Physical activity in the prevention of type 2 diabetes: the Finnish Diabetes prevention study. Diabetes. 54 (1), 158–165. doi: 10.2337/diabetes.54.1.158

PubMed Abstract | CrossRef Full Text | Google Scholar

Laaksonen, D. E., Nyyssönen, K., Niskanen, L., Rissanen, T. H., Salonen, J. T. (2005b). Prediction of cardiovascular mortality in middle-aged men by dietary and serum linoleic and polyunsaturated fatty acids. Arch. Intern. Med. 165 (2), 193–199. doi: 10.1001/archinte.165.2.193

PubMed Abstract | CrossRef Full Text | Google Scholar

Langfelder, P., Horvath, S. (2008). WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics. 9, 559. doi: 10.1186/1471-2105-9-559

PubMed Abstract | CrossRef Full Text | Google Scholar

Langmead, B., Trapnell, C., Pop, M., Salzberg, S. (2009). Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. Genome Biol. 10 (3), R25. doi: 10.1186/gb-2009-10-3-r25

PubMed Abstract | CrossRef Full Text | Google Scholar

Lau, D. C. W. (2005). Adipokines: molecular links between obesity and atheroslcerosis. Am. J. Physiol. Heart Circ. Physiol. 288 (5), H2031–H2041. doi: 10.1152/ajpheart.01058.2004

PubMed Abstract | CrossRef Full Text | Google Scholar

Lawless, N., Vegh, P., O’Farrelly, C., Lynn, D. J. (2014). The role of microRNAs in bovine infection and immunity. Front. Immunol. 5, 611. doi: 10.3389/fimmu.2014.00611

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, R., Beaudoin, F., Ammah, A. A., Bissonnette, N., Benchaar, C., Zhao, X., et al. (2015). Deep sequencing shows microRNA involvement in bovine mammary gland adaptation to diets supplemented with linseed oil or safflower oil. BMC Genomics. 16, 884. doi: 10.1186/s12864-015-1965-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, X., Huang, K., Chen, F., Li, W., Sun, S., Shi, X. E., et al. (2016). Verification of suitable and reliable reference genes for quantitative real-time PCR during adipogenic differentiation in porcine intramuscular stromal-vascular cells. Animal. 10 (6), 947–952. doi: 10.1017/S1751731115002748

PubMed Abstract | CrossRef Full Text | Google Scholar

Lian, S., Guo, J. R., Nan, X. M., Ma, L., Loor, J. J., Bu, D. P. (2016). MicroRNA Bta-miR-181a regulates the biosynthesis of bovine milk fat by targeting ACSL1. J. Dairy Sci. 99 (5), 3916–3924 doi: 10.3168/jds.2015-10484

PubMed Abstract | CrossRef Full Text | Google Scholar

Love, M. I., Huber, W., Anders, S. (2014). Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 15 (12), 550. doi: 10.1186/s13059-014-0550-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Moloney, F., Toomey, S., Noone, E., Nugent, A., Allan, B., Loscher, C. E., et al. (2007). Antidiabetic effects of cis-9, trans-11-conjugated linoleic acid may be mediated via anti-inflammatory effects in white adipose tissue. Diabetes. 56 (3), 574–582. doi: 10.2337/db06-0384

PubMed Abstract | CrossRef Full Text | Google Scholar

Myeong, J. Y., Gha, Y. L., Chung, J. J., Young, H. A., Seung, H. H., Jae, B. K. (2006). Adiponectin increases fatty acid oxidation in skeletal muscle cells by sequential activation of AMP-activated protein kinase, p38 mitogen-activated protein kinase, and peroxisome proliferator-activated receptor{alpha}. Diabetes. 55 (9), 2562–2570. doi: 10.2337/db05-1322

PubMed Abstract | CrossRef Full Text | Google Scholar

Nakadera, E., Yamashina, S., Izumi, K., Inami, Y., Sato, T., Fukushima, H., et al. (2016). Inhibition of mTOR improves the impairment of acidification in autophagic vesicles caused by hepatic steatosis. Biochem. Biophys. Res. Commun. 469 (4), 1104–1110. doi: 10.1016/j.bbrc.2015.12.010

PubMed Abstract | CrossRef Full Text | Google Scholar

Nakamura, M. T., Nara, T. Y. (2003). Essential fatty acid synthesis and its regulation in mammals. Prostaglandins Leukot. Essent. Fatty Acids. 68 (2), 145–150. doi: 10.1016/S0952-3278(02)00264-8

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Ortega, F. J., Mercader, J. M., Catalán, V., Moreno-Navarrete, J. M., Pueyo, N., Sabater, M., et al. (2013). Targeting the circulating microRNA signature of obesity. Clin. Chem. 59 (5), 781–792. doi: 10.1373/clinchem.2012.195776

PubMed Abstract | CrossRef Full Text | Google Scholar

Pollard, T. D., Cooper, J. A. (1986). Actin and actin-binding proteins. Annu. Rev. Biochem. 55, 987–1035. doi: 10.1146/annurev.bi.55.070186.005011

PubMed Abstract | CrossRef Full Text | Google Scholar

Reverter, A., Chan, E. K. F. (2008). Combining partial correlation and an information theory approach to the reversed engineering of gene co-expression networks. Bioinformatics. 24 (21), 2491–2497. doi: 10.1093/bioinformatics/btn482

PubMed Abstract | CrossRef Full Text | Google Scholar

Reverter, A., Hudson, N. J., Nagaraj, S. H., Pérez-Enciso, M., Dalrymple, B. P. (2010). Regulatory impact factors: unraveling the transcriptional regulation of complex traits from expression data. Bioinformatics. 26 (7), 896–904. doi: 10.1093/bioinformatics/btq051

PubMed Abstract | CrossRef Full Text | Google Scholar

Romao, J. M., Jin, W., He, M., McAllister, T., Guan, L. L. (2014). MicroRNAs in bovine adipogenesis: genomic context, expression and function. BMC Genomics. 15, 137. doi: 10.1186/1471-2164-15-137

PubMed Abstract | CrossRef Full Text | Google Scholar

Scaglia, N., Igal, R. A. (2005). Stearoyl-CoA desaturase is involved in the control of proliferation, anchorage-independent growth, and survival in human transformed cells. J. Biol. Chem. 280 (27), 25339–25349. doi: 10.1074/jbc.M501159200

PubMed Abstract | CrossRef Full Text | Google Scholar

Schwartz, G. J., Fu, J., Astarita, G., Li, X., Gaetani, S., Campolongo, P., et al. (2008). The lipid messenger OEA links dietary fat intake to satiety. Cell Metab. 8 (4), 281–288. doi: 10.1016/j.cmet.2008.08.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Smith, S. B., Lunt, D. K., Chung, K. Y., Choi, C. B., Tume, R. K., Zembayashi, M. (2006). Adiposity, fatty acid composition, and delta-9 desaturase activity during growth in beef cattle. Anim. Sci. J. 77 (5), 478–486. doi: 10.1111/j.1740-0929.2006.00375.x

CrossRef Full Text | Google Scholar

Smith, S., Gill, C., Lunt, D., Brooks, M. (2009). Regulation of fat and fatty acid composition in beef cattle. Asian-australas. J. Anim. Sci. 22 (9), 1225–1233. doi: 10.5713/ajas.2009.r.10

CrossRef Full Text | Google Scholar

Soares-Silva, M., Diniz, F. F., Gomes, G. N., Bahia, D. (2016). The mitogen-activated protein kinase (MAPK) pathway: role in immune evasion by trypanosomatids. Front. Microbiol. 7, 183. doi: 10.3389/fmicb.2016.00183

PubMed Abstract | CrossRef Full Text | Google Scholar

Soyeurt, H., Dardenne, P., Dehareng, F., Bastin, C., Gengler, N. (2008). Genetic parameters of saturated and monounsaturated fatty acid content and the ratio of saturated to unsaturated fatty acids in bovine milk. J. Dairy Sci. 91 (9), 3611–3626. doi: 10.3168/jds.2007-0971

PubMed Abstract | CrossRef Full Text | Google Scholar

Suzuki, A., Hiraoka, N., Suzuki, M., Angata, K., Misra, A. K., McAuliffe, J., et al. (2001). Molecular cloning and expression of a novel human β-Gal-3-O-sulfotransferase that acts preferentially on N-acetyllactosamine in N- and O-glycans. J. Biol. Chem. 276 (26), 24388–24395. doi: 10.1074/jbc.M103135200

PubMed Abstract | CrossRef Full Text | Google Scholar

Tanaka, K. (2009). The proteasome: overview of structure and functions. Proc. Jpn. Acad., Ser. B. 85 (1), 12–36. doi: 10.2183/pjab.85.12

CrossRef Full Text | Google Scholar

Thomson, D. M., Herway, S. T., Fillmore, N., Kim, H., Brown, J. D., Barrow, J. R., et al. (2008). AMP-activated protein kinase phosphorylates transcription factors of the CREB family. J. Appl. Physiol. 104 (2), 429–438. doi: 10.1152/japplphysiol.00900.2007

PubMed Abstract | CrossRef Full Text | Google Scholar

Valsta, L. M., Tapanainen, H., Männistö, S. (2005). Meat fats in nutrition. Meat Sci. 70 (3), 525–530. doi: 10.1016/j.meatsci.2004.12.016

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, H., Xi, Q., Wu, G. (2016). Fatty acid synthase regulates invasion and metastasis of colorectal cancer via Wnt signaling pathway. Cancer Med. 5 (7), 1599–1606. doi: 10.1002/cam4.711

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Wicik, Z., Gajewska, M., Majewska, A., Walkiewicz, D., Osińska, E., Motyl, T. (2016). Characterization of microRNA profile in mammary tissue of dairy and beef breed heifers. J. Anim. Breed. Genet. 133 (1), 31–42. doi: 0.1111/jbg.12172

Google Scholar

Wood, J. D., Enser, M., Fisher, A. V., Nute, G. R., Sheard, P. R., Richardson, R. I., et al. (2008). Fat deposition, fatty acid composition and meat quality: a review. Meat Sci. 78 (4), 343–358 doi: 10.1016/j.meatsci.2007.07.019

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, W. W., Sun, X. F., Tong, H. L., Wang, Y. H., Li, S. F., Yan, Y. Q., et al. (2016). Effect of differentiation on microRNA expression in bovine skeletal muscle satellite cells by deep sequencing. Cell. Mol. Biol. Lett. 21, 8. doi: 10.1186/s11658-016-0009-x

PubMed Abstract | CrossRef Full Text | Google Scholar