Abstract
Meat quality is a complex trait that is influenced by genetic and environmental factors, which includes mineral concentration. However, the association between mineral concentration and meat quality, and the specific molecular pathways underlying this association, are not well explored. We therefore analyzed gene expression as measured with RNA-seq in Longissimus thoracis muscle of 194 Nelore steers for association with three meat quality traits (intramuscular fat, meat pH, and tenderness) and the concentration of 13 minerals (Ca, Cr, Co, Cu, Fe, K, Mg, Mn, Na, P, S, Se, and Zn). We identified seven sets of co-expressed genes (modules) associated with at least two traits, which indicates that common pathways influence these traits. From pathway analysis of module hub genes, we further found an over-representation for energy and protein metabolism (AMPK and mTOR signaling pathways) in addition to muscle growth, and protein turnover pathways. Among the identified hub genes FASN, ELOV5, and PDE3B are involved with lipid metabolism and were affected by previously identified eQTLs associated to fat deposition. The reported hub genes and over-represented pathways provide evidence of interplay among gene expression, mineral concentration, and meat quality traits. Future studies investigating the effect of different levels of mineral supplementation in the gene expression and meat quality traits could help us to elucidate the regulatory mechanism by which the genes/pathways are affected.
Introduction
Meat is an important source of nutrients in the human diet. Meat quality traits such as intramuscular fat content (IMF), mineral concentration, and fatty acid profile influence consumer purchase decision (; ) and human health (Pighin et al., 2016). Mineral deficiency, mainly iron and zinc (Ritchie and Roser, 2018), and protein deficiency (), have been reported as worldwide health hazards. In addition, IMF, meat pH, and muscle mineral concentration also affect meat tenderness, flavor, and juiciness, which are major sensory traits related to eat satisfaction (; ; Pannier et al., 2014).
Brazil is one of the largest exporters of meat and meat products, and the Brazilian cattle herd is mainly composed of Nelore and its crosses (). Despite being well adapted to tropical climate, Nelore cattle has typically less tender and marble meat when compared with European breeds due to several genetic and environmental factors (; Tizioto et al., 2015). Genome-wide association (GWAS) of SNPs (Tizioto et al., 2013, 2015; ) and copy-number variations (CNVs) (Silva et al., 2016) in conjunction with transcriptomic studies (; Silva-Vignato et al., 2017; ; ), have illustrated the genetic factors affecting complex traits in Nelore. However, growing evidence suggested interplay among gene expression, mineral concentration, and meat quality traits, which are still unclear.
Multi-omic data integration has been useful to reveal potential causal and regulatory mechanisms underlying complex animal production, reproduction and welfare traits [reviewed in Suravajhala et al. (2016)]. Integrating genomic, transcriptomic, and phenotype data has contributed to an improved understanding of complex traits by identifying regulatory candidate genes and biological functions (Ponsuksili et al., 2013; ; ; ). Based on that, carried out a GWAS combined with gene network analysis for association with the carcass, meat quality traits and mineral concentration. Among the identified pathways, the authors pointed out calcium-related processes, apoptosis, and TGF-beta signaling involved with these traits.
Genome-wide association and differential gene expression analyses have been fruitful in investigating the role of genes in complex phenotypes. However, biological systems are a result of complex interactions among genes and multiple regulatory mechanisms, which are not explored in the above-mentioned studies. To address the relationship between transcriptome and traits, co-expression networks have been successfully employed. This approach allows to identify and cluster highly connected genes and associate them to the phenotypes, shedding light on the common pathways underlying these traits as well as the main regulators (). To date, there is no information about this approach integrating meat quality traits and mineral concentration in beef cattle. In addition, we still have a lack of knowledge about the interplay among gene expression, mineral concentration, and meat quality traits. Thus, to explore regulatory pathways, putative gene regulators, and to study their relationship with muscle and mineral metabolism in Nelore skeletal muscle, we integrated gene expression, eQTL variation, mineral concentration (macro and micro minerals), and meat quality traits (intramuscular fat, shear force, and meat pH) based on a network approach.
Materials and Methods
Ethics Statement
The Institutional Animal Care and Use Committee (IACUC) from the Empresa Brasileira de Pesquisa Agropecuária (EMBRAPA – Pecuária Sudeste) approved all experimental procedures involving the animals used in this study.
Animals and Phenotyping
A total of two hundred Nelore steers (produced at Embrapa Pecuária Sudeste, São Carlos – Brazil) were used in this study. The experimental design, production system, and animal management were previously described (Tizioto et al., 2015; ). Briefly, animals were raised in the grazing system until 21 months of age when they were taken to three feedlots under similar nutritional and sanitary management. The Nelore steers with an average age of 25 months were harvested at commercial facilities after about 90 days of feeding and the Longissimus thoracis (LT) muscle samples were collected.
The steaks (2.5 cm) harvested as a cross-section of the LT muscle (11th and 13th ribs) collected at slaughter were used to measure the beef quality traits as described (Tizioto et al., 2013; ). The traits evaluated were tenderness (Warner-Bratzler shear force – WBSF7, kg) measured 7 days after slaughter, meat pH measured 24 h after slaughter along with intramuscular fat (IMF%) (Tizioto et al., 2013).
Tissue samples were used for total RNA extraction () and mineral measurement (Tizioto et al., 2014). The concentration of macro minerals [calcium (Ca), magnesium (Mg), phosphorus (P), potassium (K), sodium (Na), sulfur (S)] and micro minerals [chromium (Cr), cobalt (Co), copper (Cu), manganese (Mn), selenium (Se), iron (Fe), and zinc (Zn)] were measured using inductively coupled plasma-optical emission spectrometry (ICP OES; Vista Pro-CCD ICP OES1, radial view, Varian, Mulgrave, Australia) as described by Tizioto et al. (2014).
Genome Expression Profile, Sequencing, and Data Processing
The LT muscle samples were collected immediately after slaughter, snap frozen in liquid nitrogen and kept at -80°C until RNA extraction. To extract RNA, approximately 100 mg of frozen tissue was used, and total RNA was purified using Trizol® standard protocol (Life Technologies, Carlsbad, CA, United States). The mRNA concentration and quality were evaluated in the Bioanalyzer 2100® (Agilent, Santa Clara, CA, United States).
The Illumina TruSeq® RNA Sample Preparation Kit v2 Guide (San Diego, CA, United States) protocol was used to generate cDNA libraries for each sample using 2 μg of total RNA as input. Library preparation and sequencing were conducted by ESALQ Genomics Center (Piracicaba, São Paulo, Brazil). cDNA libraries were purified and validated using Agilent 2100 Bioanalyzer (Santa Clara, CA, United States). Paired-end (PE) sequencing was performed on Illumina Hiseq 2500® (San Diego, CA, United States) platform following the standard protocols. The samples were multiplexed and run on multiple lanes to obtain 2 × 100 bp reads.
The PE reads were filtered using the Seqyclean package version 1.4.13 (1Zhbannikov et al., 2017), which removed all reads with a mean quality under 24, length under 65 bp, as well as the adapter sequences. Quality control (QC) of raw RNA-Seq reads was carried out with FastQC version 0.11.2 (2) and MultiQC version1.4 (3).
Read mapping and gene counting were carried out by STAR aligner version 2.5.4b () using a reference genome (Bos taurus, ARS-UCD1.2) and gene annotation file (release 106) obtained from NCBI (NCBI, 2018). One sample with mapping rate lower than 70% was removed out for further analyses.
The data editing was done using the Bioconductor package edgeR version 3.20.9 (Robinson et al., 2010). Taking into account that low expressed genes are less reliable and indistinguishable from sampling noise (Tarazona et al., 2015), the read counts per gene were normalized to counts per million (cpm function). The genes with less than one cpm in more than 90% of the samples were filtered out. Gene counts were normalized applying the variance stabilizing transformation (VST) from DESeq2 version 1.18.1 ().
Potential biases due to technical variation in gene expression among samples were evaluated by applying a Principal Component Analysis (PCA) and hierarchical clustering on normalized data using NOISeq version 2.22.1 (Tarazona et al., 2015). A linear model was fitted in order to adjust the gene expression matrix for batch effect (flow cell). To this end, the removeBatchEffect function from Limma (version 3.34.9) R package (Ritchie et al., 2015) was used. Three samples were identified as outliers. Thus, 12 known housekeeping genes were selected based on the literature (ACTB, API5, EIF2B2, GAPDH, GUSB, HMBS, PGK1, PPIA, RPL13A, VAPB, YWHAZ) to evaluate their variability on the samples. The housekeeping genes expression confirmed these samples as outliers, and therefore, they were filtered out.
Network Gene Co-expression Analysis
A co-expression approach was applied using the WGCNA R package version 1.63 (). The method adopted for constructing the networks included two steps: First, a similarity co-expression network was calculated with Pearson’s correlation for all genes, followed by transformation to a signed adjacency matrix (AM) by using the soft thresholding power β, to which co-expression similarity is raised. Based on the criteria of approximating scale-free topology, we chose the power of β = 12 such that the resulting network satisfies the scale-free topology (linear regression model fitting index R2 = 0.80).
Outlier animals (n = 2) were identified based on hierarchical clustering and filtered out (as they had a lower number of counts compared to other samples) after WGCNA quality control, as suggested by the WGCNA authors. Accordingly, 194 animals and 11,996 genes were used to construct an undirected, signed network. Topological overlap measure (TOM) was computed from AM where TOM was converted to dissimilarity TOM. Based on TOM dissimilarity, we used the dynamic tree cut v.1.63.1 () to identify the modules as the branches of the resulting dendrogram. As parameters, the minimum size per module was set to 50 genes with a high sensitivity to split the clusters (deepSplit = 4). Genes with a similar expression pattern across samples were grouped into the same module and arbitrarily labeled by number.
Weighted Gene Co-expression Network Analysis was used for summarizing the obtained modules by a concept of eigengene. Eigengenes are the first principal component of the expression matrix for each module and represent the weighted average of expression profile for each module. Modules highly correlated were merged based on the ME dissimilarity threshold of 0.2 leading to the final set of modules for constructing the network.
Trait Association Analysis and Module Selection
After the phenotypic data were mean-centered and scaled, a linear model was fitted to analyze the association between the expression profiles of the MEs and the phenotypes (). The model included the place of birth, the season of production, and animal’s age, according to the equation:
Where:
yijkl: is the expression level of the eigengene in each module (n = 23);
μ : is the intercept of ME;
Ci : is the fixed effect for the place of birth (three levels = CPPSE, IMA, NOHO);
Gj: is the fixed effect for the season of production (three levels = 2009, 2010, 2011);
Ak : is the covariate for the animal’s age, in days;
Tl : is the trait observation for each animal;
εijkl: is the random residual effect associated with each observation.
Modules associated with at least two beef quality or mineral traits (p ≤ 0.05) were selected for further analyses.
Pathway Over-Representation Analysis
Pathway analysis was performed using ClueGO version 2.5.1 to identify gene KEGG pathways over-represented in the selected modules (). Redundant terms were grouped based on the kappa score = 0.4 (). The p-value was calculated and corrected with a Bonferroni step down. Only pathways with a p-value (pV) p ≤ 0.05 were selected. These analyses were carried out based on the B. taurus annotation, and the network visualization was performed on Cytoscape version 3.6.1 (Shannon et al., 2003).
Hub Gene Selection
Highly connected genes (hub genes) are supposed to be the main regulators in the network and have a pivotal biological role concerning the associated trait (, ). Hub genes in the associated modules were selected based on the module membership ≥ 0.8 (). Among them, hub genes partaking in over-represented biological pathways previously identified were retained. Moreover, over-representation pathway analysis including all hub genes was applied following the approach previously described.
Integration of eQTL and Co-expression Modules
A list of eQTLs from the same population and dataset () evaluated in this work was provided. The dataset included 1,268 cis- and 10,334 trans-eQTLs based on the association between 461,466 SNPs and the expression level of 11,808 genes from 192 animals. Since the eQTLs have a known effect on gene expression, the eQTLs that target the hub genes (MM ≥ 0.8) in the selected modules were evaluated. A Fisher’s exact test was applied to assess the module under/over-representation (FDR = 0.05).
Results
We applied a network-based approach to identify relevant genes and pathways associated with meat quality and mineral concentration in Nelore cattle (Figure 1). Based on the transcriptomic profiles of skeletal muscle samples of 194 steers, we constructed a signed weighted gene co-expression network with WGCNA (). From co-expressed modules and pathway analysis, we thereby identified several hub genes significantly associated with meat quality traits and mineral concentration.
FIGURE 1
Descriptive Statistics and Correlation Estimates
We analyzed gene expression levels as measured with RNA-seq for association with three meat quality traits (intramuscular fat, meat pH, and tenderness) and the concentration of 13 minerals (Ca, Cr, Co, Cu, Fe, K, Mg, Mn, Na, P, S, Se, and Zn) available for a varying number of samples (ranging from 57 to 194, Supplementary Table S1). The genetic variance and heritability for the traits evaluated here, obtained from this population, ranged from low to moderate as previously published (Tizioto et al., 2013, 2015). A summary of descriptive statistics for each trait is in Supplementary Table S1 and Figure 2.
FIGURE 2

Box plot of mineral concentration (macro and micro) and meat quality traits. IMF% – Percentage of intramuscular fat content; WBSF7 – Warner-Bratzler Shear Force after 7 days of meat aging; The data are in log10 scale.
We performed clustering analysis to identify similarities between traits (Figure 3 – top). We identified four clusters as follows: cluster 1 (WBSF7 and Cr), cluster 2 (Co, Cu, Mn, and IMF), cluster 3 (Fe, Ca, S, Zn, Na, P, Mg, and K), and cluster 4 (pH and Se). The pair-wise correlation within all traits is provided in Supplementary Figure S1. Significant and strong correlation ranged from 0.45 to 0.99 among minerals in the cluster 3 (p ≤ 0.05). We identified positive correlation among IMF with some minerals (Ca = 0.25, Cu = 0.23, Mn = 0.24, K = 0.17, Na = 0.3, S = 0.18, and Zn = 0.23) (p ≤ 0.05). Meat pH was positively correlated with Se (r = 0.29), whereas negatively associated with Fe (-0.17), Mg (-0.22), P (-0.25), K (-0.21), Na (-0.26), S (-0.17), and Zn (-0.22) (Supplementary Figure S1). No significant correlation was observed between tenderness (WBSF7), IMF, and meat pH.
FIGURE 3

Hierarchical clustering of phenotypic correlation between traits (top) and module-trait association analysis (bottom). Modules are labeled by number on the y-axis with the number of contained genes in parenthesis. Each column represents a trait as indicated on the corresponding dendrogram branch. For significantly associated modules (p ≤ 0.05), the coefficient from the linear model is given within the cell.
Data Processing and Co-expression Network Construction
On average, a total of 13 million of 100 bp paired-end reads per sample were generated. Around 96.71% of unique reads were mapped to the reference B. taurus genome (ARS-UCD1.2). Taking into account that low expressed genes are less reliable and indistinguishable from sampling noise (Tarazona et al., 2015), we filtered out the genes with less than one cpm in more than 90% of the samples. In addition, four samples were removed because they had a mapping rate lower than 70% or showed high variability on the housekeeping genes expression (see methods). Thus, we used 11,996 genes and 194 samples for the co-expression analysis.
Considering the WGCNA assumptions, the weighted network starts from the level of thousands of genes, identifies modules of co-expressed genes, summarizes the module expression profile as the first principal component (ME), and relates the MEs with the trait of interest (
From clustering 11,996 genes with WGCNA, we obtained 23 modules labeled by number (Figure 3). The module size ranged from 69 genes (M9) to 2,008 genes (M14) (Figure 3 – bottom). The proportion of variance explained by the eigengenes ranged from 0.18 (M20) to 0.53 (M5) (Supplementary Table S2).
Trait Association and Pathway Enrichment Analysis
We performed an association analysis to identify the relationship between network and traits. This analysis measures the strength of the effect and the direction of the association between the module (eigengenes) and the trait. Thus, if the association is positive, it means the trait increases with increasing “eigengene expression” or vice-versa. We selected seven modules (M1, M5, M6, M7, M8, M9, and M17), associated with at least two traits (p ≤ 0.05) (Figure 3 – bottom) once we also want to point out shared pathways among traits. We found the highest number of significant associated modules between M5 (ten associations; negative with IMF, and the concentration of Mn, Fe, Ca, S, Zn, Na, P, Mg, and K), followed by M8 (nine associations; positive with Cr, negative with IMF, and the minerals of cluster 3, except Zn). The average expression profile of M17 module showed association with three traits (positive with WBSF7, Co, and Mn) along with M7 (negative with Na and IMF, and positive with Cr concentration). For the modules M1, M6, and M9, we found an association with two traits. We identified a positive association among M6 and M9 with Cr concentration while a negative association was observed between M9 with IMF, and M6 with Fe concentration. M1 was positively associated with the concentration of Cr and Co. The modules with none or only one trait association were not included for further analysis.
The module membership values for all the genes for selected modules are given in Supplementary Table S3. We carried out a pathway over-representation analysis on ClueGo version 2.5.1 for the seven selected modules (Table 1 and Supplementary Table S4) to identify meaningful metabolic pathways involved with meat quality traits and mineral concentration. We detected several pathways (p-Value ≤ 0.05, group p-value corrected with Bonferroni step down) mainly related to energy and protein metabolism, such as AMPK and mTOR signaling pathways.
Table 1
| Modulea | eQTLsb | TGEc | Hub genesd | Enriched pathwayse |
|---|---|---|---|---|
| M1 (78) | 13 | 8 (4) | (10) RSAD2, LOC100139670, MX1, IFIH1, OAS1Y, DDX58, EIF2AK2, IRF7, LOC512486, IFI16 |
|
| M5 (88) | 45 | 15 (3) | (9) PTPRC, CTSS, DOCK2, IL10RA, FERMT3, ITGB2, LCP2, CORO1A, PTPRJ |
|
| M6 (190) | 89 | 20 (6) | (20) GNAI1, PLIN1, FABP4, PPP1R1B, PCK2, ADIPOQ, PDE3B, TKT, MGST1, LIPE, ELOVL6, MRAS, ACACA, GNG2, HACD2, FASN, FBP1, ELOVL5, PCK1, G6PD |
|
| M7 (704) | 276 | 76 (2) | (11) CAV1, COL4A2, GNAI2, SEPT9, TNFRSF1A, YWHAB, COL4A1, TMSB4X, SPTAN1, MSN, PARVA |
|
| M8 (1,200) | 414 | 129 (5) | (10) COL1A2, DCN, FN1, COL1A1, FGFR1, COL3A1, CD44, ITGB5, DSE, CTSK |
|
| M9 (69) | 16 | 9 (2) | (1) JUND |
|
| M17 (975) | 126 | 66 (2) | (21) ASH1L, BIRC6, MED13, HERC1, MED13L, KMT2A, HUWE1, KMT2C, ATM, HERC2, KMT2D, MED1, GNAQ, SMC1A, TPR, NSD1, UBR5, ARHGAP5, KAT6A, PRDM2, ITCH |
|
| Total – 3304 | 979f | 323 (24) | 82 | |
Module characterization.
The table shows hub genes and eQTL information for each module found to be significantly associated with two or more traits in Figure 3. aSelected modules with the number of contained genes in parenthesis. beQTLs – Number of eQTLs associated with genes in a module (Based on
Hub Gene Selection, Pathway Analysis, and Integration With eQTLs
Highly connected genes are likely to play an important role both in the network’s topology and biological pathways. In this way, we combined a pathway-based gene analysis for each selected module (Supplementary Table S4) and gene connectivity measure (MM ≥ 0.8) (Supplementary Table S3) selecting 82 hub genes (Table 1, see methods). Further, taking advantage of an eQTL study carried out in the same population (
To gain further insights into their functions as well as to integrate the pathways among the modules, we carried out a KEGG pathway analysis. Considering a kappa score = 0.4 and p-Value ≤0.05 (Figure 4 and Supplementary Table S6), we clustered the identified pathways into eight groups. The pathways related to energy metabolism were clustered together and included AMPK, peroxisome proliferator-activated receptors (PPAR), insulin, glucagon, and adipocytokine signaling pathways. We also identified ubiquitin-mediated proteolysis and biosynthesis of fatty acids pathways over-represented in this network.
FIGURE 4

Network clusters based on over-represented KEGG pathways of hub genes associated with mineral concentration and meat quality traits. Functionally related groups partially overlap and are arbitrarily colored. The node size represents the pathway enrichment significance.
Discussion
In this study, we analyzed genome-wide co-expression in skeletal muscle for association with mineral concentration and meat quality traits. Skeletal muscle metabolism is an integrated system dependent on the efficient coordination of gene expression, which are tightly regulated (Smith et al., 2013). We found several co-expression modules associated with two or more minerals, meat tenderness, and IMF, which indicates that common pathways influence these traits. From pathway analysis of module hub genes, we further found an over-representation for energy and protein metabolism (AMPK and mTOR). These pathways have been reported as the main drivers regulating energy balance in muscle (Smith et al., 2013). AMPK and mTOR are metabolically linked, nutritional and hormonal responsive, with an intricated relationship with insulin, thyroid hormone (TH), and TGF-beta signaling pathways (Xu et al., 2012), which were reported here as well. In addition, these pathways have been associated with muscle development, fat deposition, and beef quality traits (
Phenotype Correlation and Co-expression Network Analysis
In agreement with previous reports, we found several minerals positively correlated with IMF, but negatively correlated with meat pH. For instance, Cu-supplemented Angus were found with reduced back fat and reduced serum cholesterol level (
Co-expression analysis resulted in 23 modules from which we considered seven modules for further analysis based on their association with at least two traits. The genes in modules like M5, M7, M8, and M9 were associated with IMF and several minerals suggesting a certain extent of co-regulation. It’s well known that minerals are essential in a wide range of biological processes. Here, we provide evidences that mineral content and meat quality traits are interrelated, as well as interplay with specific genes and pathways (as discussed below).
Variation in eQTL loci can explain a substantial fraction of variation observed on the gene expression level (Wang and Michoel, 2016). It has been observed that variation in eQTL loci is associated with concerted expression changes of many genes in co-expression clusters, thereby also impact the phenotype. Screening the detected co-expression modules, we found 323 genes affected by at least one eQTL. Despite 132 eQTLs targeting more than one module, most of the eQTLs were module-specific. However, no significant over-/under-representation (Fisher’s exact test) was detectable in these modules, suggesting that other regulatory mechanisms are involved. Despite that, the expression level of six hub genes was found affected by trans-eQTLs. These genes are involved with lipid metabolism [fatty acid synthase (FASN), and ELOVL fatty acid elongase 5 (ELOVL5), phosphodiesterase 3B (PDE3B)], immune system [lymphocyte cytosolic protein 2 (LCP2), and interleukin 10 receptor subunit alpha (IL10RA)], and actin remodeling (dedicator of cytokinesis 2 – DOCK2).
Pathway Analysis
Over-representation pathway analysis in the selected modules (Table 1 and Supplementary Table S4) yielded glycosaminoglycan biosynthesis and degradation, lysosome, and steroid biosynthesis in the M8 module. Phagosome, cell adhesion molecular pathways, and NOD-like receptor signaling pathway were found enriched in M1 and M5. For the M17, enriched pathways included protein synthesis pathways such as mTOR, PI3K-Akt, TH, and AMPK signaling. We also found protein degradation pathways enriched in the M17 module such as ubiquitin-mediated proteolysis. TGF-beta signaling and osteoclast differentiation were enriched in the M9 module. Energy metabolism pathways were found enriched in M6, including glycolysis, fatty acid biosynthesis, AMPK, and insulin signaling. Ras, PI3K-Akt signaling pathways, and protein processing were found enriched in M7.
We also carried out cross-module enrichment analysis considering all hub genes, which indicated that the AMPK signaling pathway plays an important role for muscle mineral metabolism and meat quality traits. The genes of the AMPK pathway were also associated with IMF, Cr, and Fe. Furthermore, the AMPK pathway was also found enriched in genes of M17 (associated with WBSF7, Co, and Mg) and M6 (associated with Cr and Fe concentration).
Energy and Lipid Metabolism
AMP-activated protein kinase signaling is a major regulator of the cellular energy status, protein metabolism, and muscle metabolism (
The co-expression of these genes, as well as the negative association between Fe concentration and lipid metabolism, were reported in our previous RNA-Seq work where FASN, THRSP, and FABP4 were shown to be downregulated in animals with low Fe concentration in muscle (
In addition to factors that increase the intracellular cyclic AMP level (Omar et al., 2009), Cr increases AMPK activity and positively affects the insulin sensitivity in skeletal muscle cells (
Supplementing goats with Cr decreased the expression level of ACACA, FASN, and FABP4 (Sadeghi et al., 2015) as measured by RT-PCR. Furthermore, increased Longissimus muscle area and reduced fat thickness was associated with a downregulation of ACACA expression in Cr-supplemented goats (
Muscle Development, Structure, and Proteolysis
As part of the TGF-beta signaling pathway, we identified the transforming growth factor beta 3 (TGFB3), which is involved in muscle proliferation, differentiation, and growth (Nishimura, 2015). However, muscle hypertrophy results from a balance of protein turnover in which AMPK signaling negatively affects the protein synthesis (
We found ubiquitin-mediated proteolysis enriched across modules as well as for genes in the M17 (Supplementary Table S3), which was associated with WBSF7, Co, and Mn. Proteolytic enzymes are important for protein turnover and postmortem meat aging (
Conclusion
We demonstrated transcriptional relationships among mineral concentration and meat quality traits in the skeletal muscle of Brazilian Nelore cattle. We identified 82 hub genes across seven co-expression modules which seem to be critical for this interplay. The AMPK and mTOR signaling pathways were hereby found to link mineral and muscle metabolism in Nelore cattle. Future studies investigating different levels of mineral supplementation, the mineral interaction, and their effect in the gene expression and meat quality traits could help us to elucidate the regulatory mechanism by which the genes/pathways are affected.
Statements
Data availability statement
All relevant data are within the paper and its Supporting Information files. All sequencing data is available in the European Nucleotide Archive (ENA) repository (EMBL-EBI), under accession PRJEB13188, PRJEB10898, and PRJEB19421 (https://www.ebi.ac.uk/ena/submit/sra/). All additional datasets generated and analyzed during this study may be available upon request from the corresponding author on reasonable request.
Ethics statement
The Institutional Animal Care and Use Committee (IACUC) from the Empresa Brasileira de Pesquisa Agropecuária (EMBRAPA – Pecuária Sudeste) approved all experimental procedures involving the animals used in this study.
Author contributions
WD, PT, LR, LC, and HK conceived the idea of this research. WD, GM, LG, FB, and AC, JA, PdO, PT carried out the bioinformatics and data analysis. WD, GM, PB, FB, HK, JA, LC, LR collaborated with the interpretation of results, discussion and review the manuscript. WD, PB, and GM drafted the manuscript. All authors have read and approved the final manuscript and agreed to be accountable for the content of the work.
Funding
This study was conducted with funding from EMBRAPA (Macroprograma 1, 01/2005), FAPESP (grant# 2012/23638-8), and by the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior – Brasil (CAPES) – Finance Code 001. LR and LC were granted CNPq fellowships. WD was granted by São Paulo Research Foundation (FAPESP) grant# 2015/09158-1 and grant# 2017/20761-7 scholarships. Federal University of São Carlos (PROAP/PNPD) granted funding for publishing.
Acknowledgments
We are thankful to Bruno G. N. Andrade for the server management and support; Dr. Ana Rita Araújo Nogueira and Dr. Caio F. Gromboni for the mineral data; and the Technical University of Denmark (DTU Compute) for accepting the first author as a visiting scholar.
Conflict of interest
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.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2019.00210/full#supplementary-material
FIGURE S1Correlation matrix of mineral concentration and meat quality traits in Nelore cattle. Each cell displays the correlation value when significant (p ≤ 0.05). The matrix is color-coded by correlation according to the color legend.
TABLE S1Summary statistics of meat quality traits and mineral concentration in Nelore cattle.
TABLE S2The proportion of variance explained by the module eigengene (MEs).
TABLE S3Gene list and module membership for each selected module. Spreadsheet tabs are divided by module.
TABLE S4Summary of pathway analysis from ClueGo for genes clustered into the selected modules. Spreadsheet tabs are divided by module.
TABLE S5Genes targeted by eQTLs for each selected module. Genes with MM ≥ 0.8 are highlighted in bold. Spreadsheet tabs are divided by module.
TABLE S6Summary of pathway analysis from ClueGo for hub genes.
Abbreviations
- AMPK
AMP-activated protein kinase
- CPM
Counts per million
- ECM
Extra Cellular Matrix
- IMF
Intramuscular Fat Content
- ME
Module eigengene
- MM
Module Membership
- QC
Quality Control
- WBSF7
Warner-Bratzler Shear Force after 7 days of meat aging
- WGCNA
Weighted Gene Co-expression Network Analysis
Footnotes
1.^https://github.com/ibest/seqyclean
2.^https://www.bioinformatics.babraham.ac.uk/projects/fastqc/
References
1
ABIEC (2018). Perfil da pecuária no Brazil. São Paulo. Available at: http://abiec.siteoficial.ws/images/upload/sumario-pt-010217.pdf
2
AhlbergC. M.SchiermiesterL. N.HowardT. J.CalkinsC. R.SpanglerM. L. (2014). Genome wide association study of cholesterol and poly- and monounsaturated fatty acids, protein, and mineral content of beef from crossbred cattle.Meat Sci.98804–814. 10.1016/j.meatsci.2014.07.030
3
AndersS.HuberW. (2010). Differential expression analysis for sequence count data.Genome Biol.11:R106. 10.1186/gb-2010-11-10-r106
4
AndrewsS. (2010). FastQC: a Quality Control Tool for High Throughput Sequence Data. Available at: https://www.bioinformatics.babraham.ac.uk/projects/fastqc/
5
BindeaG.MlecnikB.HacklH.CharoentongP.TosoliniM.KirilovskyA.et al (2009). ClueGO: a Cytoscape plug-in to decipher functionally grouped gene ontology and pathway annotation networks.Bioinformatics251091–1093. 10.1093/bioinformatics/btp101
6
CalvoL.ToldráF.RodríguezA. I.López-BoteC.ReyA. I. (2017). Effect of dietary selenium source (organic vs. mineral) and muscle pH on meat quality characteristics of pigs.Food Sci. Nutr.594–102. 10.1002/fsn3.368
7
CammackR.WrigglesworthJ. M.BaumH. (1990). “Iron-dependent enzymes in mammalian systems,” inIron:Transport and Storage, edsPonkaP.SchulmanH. M.WoodworthR. C.RichterG. W. (Boca Raton, FL: CRC Press), 17–39.
8
CesarA. S. M.RegitanoL. C. A.KoltesJ. E.Fritz-WatersE. R.LannaD. P. D.GasparinG.et al (2015). Putative regulatory factors associated with intramuscular fat content.PLoS One10:e0128350. 10.1371/journal.pone.0128350
9
CesarA. S. M.RegitanoL. C. A.MourãoG. B.TullioR. R.LannaD. P. D.NassuR. T.et al (2014). Genome-wide association study for intramuscular fat deposition and composition in Nellore cattle.BMC Genet.15:39. 10.1186/1471-2156-15-39
10
CesarA. S. M.RegitanoL. C. A.ReecyJ. M.PoletiM. D.OliveiraP. S. N.De OliveiraG. B.et al (2018). Identification of putative regulatory regions and transcription factors associated with intramuscular fat content traits.BMC Genomics19:499. 10.1186/s12864-018-4871-y
11
ClugstonG. A.SmithT. E. (2002). Global nutrition problems and novel foods Global nutrition problems and novel foods S101.Asia Pac. J. Clin. Nutr.11100–111. 10.1046/j.1440-6047.11.s6.1.x
12
ConstantinoC.de Azambuja RibeiroE. L.BridiA. M.TarsitanoM. A.de CastroF. A. B.Fernandes JúniorF.et al (2014). Performance, carcass and meat quality of ewes supplemented with magnesium oxide.Rev. Bras. Zootec.4327–35. 10.1590/S1516-35982014000100005
13
CunninghamB. A.MoncurJ. T.HuntingtonJ. T.KinlawW. B. (1998). “Spot 14” protein: a metabolic integrator in normal and neoplastic cells.Thyroid8815–825. 10.1089/105072502761016494
14
DegermanE.AhmadF.ChungY. W.GuirguisE.OmarB.StensonL.et al (2011). From PDE3B to the regulation of energy homeostasis.Curr. Opin. Pharmacol.11676–682. 10.1016/j.coph.2011.09.015
15
DinizW. J.CoutinhoL. L.TiziotoP. C.CesarA. S. M.GromboniC. F.NogueiraA. R. A.et al (2016). Iron content affects lipogenic gene expression in the muscle of nelore beef cattle.PLoS One11:e0161160. 10.1371/journal.pone.0161160
16
DobinA.DavisC. A.SchlesingerF.DrenkowJ.ZaleskiC.JhaS.et al (2013). STAR: ultrafast universal RNA-seq aligner.Bioinformatics2915–21. 10.1093/bioinformatics/bts635
17
DuM.ShenQ. W.UnderwoodK. R.TongJ. F.ZhuM. (2009). “AMP-activated protein kinase in muscle growth, fat deposition, and meat quality,” inApplied Muscle Biology and Meat Science, edsDuM.McCormickR. J. (Boca Raton, FL: CRC Press), 360.
18
EngleT. E.SpearsJ. W.ArmstrongT. A.WrightC. L.OdleJ. (2000). Effects of dietary copper source and concentration on carcass characteristics and lipid and cholesterol metabolism in growing and finishing steers.J. Anim. Sci.78:1053. 10.2527/2000.7841053x
19
EwelsP.MagnussonM.LundinS.KällerM. (2016). MultiQC: summarize analysis results for multiple tools and samples in a single report.Bioinformatics323047–3048. 10.1093/bioinformatics/btw354
20
GeistlingerL.da SilvaV. H.CesarA. S. M.TiziotoP. C.WaldronL.ZimmerR.et al (2018). Widespread modulation of gene expression by copy number variation in skeletal muscle.Sci. Rep.8:1399. 10.1038/s41598-018-19782-4
21
GonçalvesT. M.de Almeida RegitanoL. C.KoltesJ. E.CesarA. S. M.da Silva AndradeS. C.MourãoG. B.et al (2018). Gene Co-expression analysis indicates potential pathways and regulators of beef tenderness in nellore cattle.Front. Genet.9:441. 10.3389/fgene.2018.00441
22
GraugnardD. E.PiantoniP.BionazM.BergerL. L.FaulknerD. B.LoorJ. J. (2009). Adipogenic and energy metabolism gene networks in longissimus lumborum during rapid post-weaning growth in Angus and Angus x Simmental cattle fed high-starch or low-starch diets.BMC Genomics10:142. 10.1186/1471-2164-10-142
23
HardieD. G. (2011). AMP-activated protein kinase–an energy sensor that regulates all aspects of cell function.Genes Dev.251895–1908. 10.1101/gad.17420111
24
HayS. M.McardleH. J.HayesH. E.StevensV. J.ReesW. D. (2016). The effect of iron deficiency on the temporal changes in the expression of genes associated with fat metabolism in the pregnant rat.Physiol. Rep.41–9. 10.14814/phy2.12908
25
HoffmanN. J.PenqueB. A.HabeggerK. M.SeallsW.TackettL.ElmendorfJ. S. (2015). Chromium enhances insulin responsiveness via AMPK.J. Nutr. Biochem.25565–572. 10.1016/j.jnutbio.2014.01.007.Chromium
26
JeW.KimM.ParkH.SikH.JaeM. (2006). AMPK activation increases fatty acid oxidation in skeletal muscle by activating PPARalpha and PGC-1.Biochem. Biophys. Res. Commun.340291–295. 10.1016/j.bbrc.2005.12.011
27
KneeskernS. G.DilgerA. C.LoerchS. C.ShikeD. W.FelixT. L. (2016). Effects of chromium supplementation to feedlot steers on growth performance, insulin sensitivity, and carcass characteristics.J. Anim. Sci.94217–226. 10.2527/jas.20159517
28
KoohmaraieM.KentM. P.ShackelfordS. D.VeisethE.WheelerT. L. (2002). Meat tenderness and muscle growth: Is there any relationship?Meat Sci.62345–352. 10.1016/S0309-1740(02)00127-4
29
LangfelderP.HorvathS. (2007). Eigengene networks for studying the relationships between co-expression modules.BMC Syst. Biol.1:54. 10.1186/1752-0509-1-54
30
LangfelderP.HorvathS. (2008). WGCNA: an R package for weighted correlation network analysis.BMC Bioinformatics9:559. 10.1186/1471-2105-9-559
31
LangfelderP.ZhangB.HorvathS. (2008). Defining clusters from a hierarchical cluster tree: the dynamic tree cut package for R.Bioinformatics24719–720. 10.1093/bioinformatics/btm563
32
LiJ.ZhouD.QiuW.ShiY.YangJ. J.ChenS.et al (2018). Application of weighted gene co-expression network analysis for data from paired design.Sci. Rep.81–8. 10.1038/s41598-017-18705-z
33
LiuY.DuM.LiX.ChenL.ShenQ.TianJ.et al (2016). Role of the ubiquitin-proteasome pathway on proteolytic activity in postmortem proteolysis and tenderisation of sheep skeletal muscle.Int. J. Food Sci. Technol.512353–2359. 10.1111/ijfs.13214
34
LoorJ. J. (2010). “Transcriptomics of muscle and fat development in growing cattle,” inEnergy and Protein Metabolism and Nutrition, ed.CrovettoG. M. (Italy: Wageningen Academic Publishers), 59–68.
35
MateescuR. (2014). “It is possible to genetically change the nutrient profile of beef,” inProceedings of the Beef improvement Federation Research Symposium & Annual Meeting, ed. Beef Improvement Federation (Lincoln, NE: Oklahoma State University), 87–92.
36
MateescuR. G.GarrickD. J.ReecyJ. M. (2017). Network analysis reveals putative genes affecting meat quality in Angus cattle.Front. Genet8:171. 10.3389/fgene.2017.00171
37
McCormickR. J. (2009). “Collagen,” inApplied Muscle Biology and Meat Science, edsDuM.McCormickR. J. (Boca Raton, FL: CRC Press), 360.
38
MenendezJ. A.LupuR. (2007). Fatty acid synthase and the lipogenic phenotype in cancer pathogenesis.Nat. Rev. Cancer7763–777. 10.1038/nrc2222
39
MihaylovaM. M.ShawR. J. (2011). The AMPK signalling pathway coordinates cell growth, autophagy and metabolism.Nat. Cell Biol.131016–1023. 10.1038/ncb2329
40
NajafpanahM. J.SadeghiM.ZaliA.Moradi-shahrebabakH.MousapourH. (2014). Chromium downregulates the expression of Acetyl CoA Carboxylase 1 gene in lipogenic tissues of domestic goats: a potential strategy for meat quality improvement.Gene543253–258. 10.1016/j.gene.2014.04.006
41
NCBI (2018). ARS-UCD1.2. Available at: https://www.ncbi.nlm.nih.gov/genome/annotation_euk/Bos_taurus/106/
42
NishimuraT. (2015). Role of extracellular matrix in development of skeletal muscle and postmortem aging of meat.Meat Sci.10948–55. 10.1016/j.meatsci.2015.05.015
43
OhD.LeeY.-S.LaB.LeeJ.-Y.ParkY.LeeJ.et al (2014). Identification of exonic nucleotide variants of the thyroid hormone responsive protein gene associated with carcass traits and Fatty Acid composition in korean cattle.Asian Australas. J. Anim. Sci.271373–1380. 10.5713/ajas.2014.14101
44
OmarB.Zmuda-TrzebiatowskaE.ManganielloV.GöranssonO.DegermanE. (2009). Regulation of AMP-activated protein kinase by cAMP in adipocytes: roles for phosphodiesterases, protein kinase B, protein kinase A. Epac and lipolysis.Cell. Signal.21760–766. 10.1016/j.cellsig.2009.01.015
45
PameiG.GnanarajP.SivakumarT.JawaharK.MuthuramalingamT.PothiappanP. (2014). Influence of dietary supplementation of chromium on the carcass traits of crossbred pigs.J. Adv. Vet. Anim. Res.1125–129. 10.5455/javar.2014.a22
46
PannierL.PethickD. W.BoyceM. D.BallA. J.JacobR. H.GardnerG. E. (2014). Associations of genetic and non-genetic factors with concentrations of iron and zinc in the longissimus muscle of lamb.Meat Sci.961111–1119. 10.1016/j.meatsci.2013.08.013
47
PighinD.PazosA.ChamorroV.PaschettaF.CunzoloS.GodoyF.et al (2016). A contribution of beef to human health: a review of the role of the animal production systems.Sci. World J. 2016: 8681491. 10.1155/2016/8681491
48
PilkisS. J.GrannerD. K. (1992). Molecular physiology of the regulation of hepatic gluconeogenesis and glycolysis.Annu. Rev. Physiol.54885–909. 10.1146/annurev.ph.54.030192.004321
49
PonsuksiliS.DuY.HadlichF.SiengdeeP.MuraniE.SchwerinM.et al (2013). Correlated mRNAs and miRNAs from co-expression and regulatory networks affect porcine muscle and finally meat properties.BMC Genomics14:533. 10.1186/1471-2164-14-533
50
ReverterA.HudsonN. J.WangY.TanS.-H.BarrisW.ByrneK. A.et al (2006). A gene coexpression network for bovine skeletal muscle inferred from microarray data.Physiol. Genomics2876–83. 10.1152/physiolgenomics.00105.2006
51
RitchieH.RoserM. (2018). Micronutrient Deficiency. Available at: https://ourworldindata.org/micronutrient-deficiency [accessed August 22, 2018].
52
RitchieM. E.PhipsonB.WuD.HuY.LawC. W.ShiW.et al (2015). limma powers differential expression analyses for RNA-sequencing and microarray studies.Nucleic Acids Res.43:e47. 10.1093/nar/gkv007
53
RobinsonM. D.McCarthyD. J.SmythG. K. (2010). edgeR: a Bioconductor package for differential expression analysis of digital gene expression data.Bioinformatics26139–140. 10.1093/bioinformatics/btp616
54
Ropka-molikK.KnapikJ.PieszkaM.SzmatołaT. (2017). Nutritional modification of SCD, ACACA and LPL gene expressions in different ovine tissues.Arch. Anim. Breed.60243–250. 10.5194/aab-60-243-2017
55
SadeghiM.JavadM.PanahN.RezaM.EmamiA. (2015). Transcription analysis of genes involved in lipid metabolism reveals the role of chromium in reducing body fat in animal models.J. Trace Elem. Med. Biol.3245–51. 10.1016/j.jtemb.2015.05.004
56
ShannonP.MarkielA.OzierO.BaligaN. S.WangJ. T.RamageD.et al (2003). Cytoscape: a software environment for integrated models of biomolecular interaction networks.Genome Res.132498–2504. 10.1101/gr.1239303.metabolite
57
SilvaV. H.RegitanoL. C.GeistlingerL.PértilleF.GiachettoP. F.BrassalotiR. A.et al (2016). Genome-wide detection of CNVs and their association with meat tenderness in Nelore cattle.PLoS One11:e0157711. 10.1371/journal.pone.0157711
58
Silva-VignatoB.CoutinhoL. L.CesarA. S. M.PoletiM. D.RegitanoL. C. A.BalieiroJ. C. C. (2017). Comparative muscle transcriptome associated with carcass traits of Nellore cattle.BMC Genomics18:506. 10.1186/s12864-017-3897-x
59
SmithL. R.MeyerG.LieberR. L. (2013). Systems analysis of biological networks in skeletal muscle function.Wiley Interdiscip. Rev. Syst. Biol. Med.555–71. 10.1002/wsbm.1197
60
StaniekH.WójciakR. W. (2018). The combined effect of supplementary Cr(III) propionate complex and iron deficiency on the chromium and iron status in female rats.J. Trace Elem. Med. Biol.45142–149. 10.1016/j.jtemb.2017.10.010
61
SuravajhalaP.KogelmanL. J. A.KadarmideenH. N. (2016). Multi-omic data integration and analysis using systems genomics approaches: methods and applications in animal production, health and welfare.Genet. Sel. Evol.48:38. 10.1186/s12711-016-0217-x
62
TajimaT.IijimaK.WatanabeT.YamaguchiH. (1981). The influence of calcium ions on the synthesis of collagen and glycosaminoglycans in human diploid cells in culture.Exp. Pathol.19219–225. 10.1016/S0232-1513(81)80067-9
63
TarazonaS.Furió-TaríP.TurràD.PietroA.Di NuedaM. J.FerrerA.et al (2015). Data quality aware analysis of differential expression in RNA-seq with NOISeq R/Bioc package.Nucleic Acids Res.43:e140. 10.1093/nar/gkv711
64
TiziotoP. C.DeckerJ. E.TaylorJ. F.SchnabelR. D.MudaduM. A.SilvaF. L.et al (2013). Genome scan for meat quality traits in Nelore beef cattle.Physiol. Genomics451012–1020. 10.1152/physiolgenomics.00066.2013
65
TiziotoP. C.GromboniC. F.NogueiraA. R.de SouzaM. M.MudaduM. D. A.TholonP.et al (2014). Calcium and potassium content in beef: influences on tenderness and associations with molecular markers in Nellore cattle.Meat Sci.96436–440. 10.1016/j.meatsci.2013.08.001
66
TiziotoP. C.TaylorJ. F.DeckerJ. E.GromboniC. F.MudaduM. A.SchnabelR. D.et al (2015). Detection of quantitative trait loci for mineral content of Nelore longissimus dorsi muscle.Genet. Sel. Evol.47:15. 10.1186/s12711-014-0083-3
67
VerhagenA. M.CoulsonE. J.VauxD. L. (2001). Inhibitor of apoptosis proteins and their relatives: IAPs and other BIRPs.Genome Biol.2 reviews3009.1–reviews3009.10. 10.1186/gb-2001-2-7-reviews3009
68
WangL.MichoelT. (2016). “Detection of Regulator Genes and eQTLs in Gene Networks,” inSystems Biology in Animal Production and HealthVol. 1ed.KadarmideenH. (Cham: Springer International Publishing), 1–23. 10.1007/978-3-319-43335-6_1
69
XuJ.JiJ.YanX. (2012). Cross-Talk between AMPK and mTOR in regulating energy balance.Crit. Rev. Food Sci. Nutr.52373–381. 10.1080/10408398.2010.500245
70
ZhbannikovI. Y.HunterS. S.FosterJ. A.SettlesM. L. (2017). “SeqyClean: a pipeline for high-throughput sequence data preprocessing,” inProceedings of the 8th ACM International Conference on Bioinformatics, Computational Biology, and Health Informatics ACM-BCB 17, (New York, NY: ACM), 407–416. 10.1145/3107411.3107446
Summary
Keywords
AMPK pathway, co-expression analysis, intramuscular fat, RNA sequencing, tenderness
Citation
Diniz WJS, Mazzoni G, Coutinho LL, Banerjee P, Geistlinger L, Cesar ASM, Bertolini F, Afonso J, de Oliveira PSN, Tizioto PC, Kadarmideen HN and Regitano LCA (2019) Detection of Co-expressed Pathway Modules Associated With Mineral Concentration and Meat Quality in Nelore Cattle. Front. Genet. 10:210. doi: 10.3389/fgene.2019.00210
Received
03 December 2018
Accepted
26 February 2019
Published
13 March 2019
Volume
10 - 2019
Edited by
Robert J. Schaefer, University of Minnesota Twin Cities, United States
Reviewed by
Martin Johnsson, Swedish University of Agricultural Sciences, Sweden; Fabyano Fonseca Silva, Universidade Federal de Viçosa, Brazil
Updates

Check for updates
Copyright
© 2019 Diniz, Mazzoni, Coutinho, Banerjee, Geistlinger, Cesar, Bertolini, Afonso, de Oliveira, Tizioto, Kadarmideen and Regitano.
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Luciana C. A. Regitano, luciana.regitano@embrapa.br
This article was submitted to Livestock Genomics, a section of the journal Frontiers in Genetics
Disclaimer
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.