Skip to main content

Advertisement

Comparative muscle transcriptome associated with carcass traits of Nellore cattle

Abstract

Background

Commercial cuts yield is an important trait for beef production, which affects the final value of the products, but its direct determination is a challenging procedure to be implemented in practice. The measurement of ribeye area (REA) and backfat thickness (BFT) can be used as indirect measures of meat yield. REA and BFT are important traits studied in beef cattle due to their strong implication in technological (carcass yield) and nutritional characteristics of meat products, like the degree of muscularity and total body fat. Thus, the aim of this work was to study the Longissimus dorsi muscle transcriptome of Nellore cattle, associated with REA and BFT, to find differentially expressed (DE) genes, metabolic pathways, and biological processes that may regulate these traits.

Results

By comparing the gene expression level between groups with extreme genomic estimated breeding values (GEBV), 101 DE genes for REA and 18 for BFT (false discovery rate, FDR 10%) were identified. Functional enrichment analysis for REA identified two KEGG pathways, MAPK (Mitogen-Activated Protein Kinase) signaling pathway and endocytosis pathway, and three biological processes, response to endoplasmic reticulum stress, cellular protein modification process, and macromolecule modification. The MAPK pathway is responsible for fundamental cellular processes, such as growth, differentiation, and hypertrophy. For BFT, 18 biological processes were found to be altered and grouped into 8 clusters of semantically similar terms. The DE genes identified in the biological processes for BFT were ACHE, SRD5A1, RSAD2 and RSPO3. RSAD2 has been previously shown to be associated with lipid droplet content and lipid biosynthesis.

Conclusion

In this study, we identified genes, metabolic pathways, and biological processes, involved in differentiation, proliferation, protein turnover, hypertrophy, as well as adipogenesis and lipid biosynthesis related to REA and BFT. These results enlighten some of the molecular processes involved in muscle and fat deposition, which are economically important carcass traits for beef production.

Background

Meat is the most important source of animal protein for the human diet; it consists mainly of skeletal muscle, and of varying amounts of connective tissue, implicated on its qualitative and quantitative characteristics, as well as small amounts of epithelial and nervous tissues. Meat represents the edible portion of the carcass, in other words, the part that will be destined for the final consumers and can be represented by the yield of commercial cuts [1, 2].

Commercial cuts yield is economically important since it affects the final value of the products due to the proportion of fat, muscle, and bone in the carcasses. The direct determination of meat yield is difficult in practice, therefore the measures of ribeye area (REA) and backfat thickness (BFT), sections of the Longissimus dorsi muscle, are often used as indirect measures of this trait [3,4,5].

REA and BFT are well studied traits in beef cattle due to their implication in technological and nutritional characteristics of meat products. The ribeye area is used as an indicator of degree of muscularity, edible mass of carcass and yield of cuts with high commercial value. This measure can also be associated with the length and weight of the carcass (hot carcass weight) [3, 6, 7].

The amount of BFT deposited on the carcass is related to the total body fat and plays a major role in beef’s flavor and juiciness, which is directly associated with production costs. In the meat industry, an adequate layer of fat acts as a thermal insulator during carcass cooling process, avoiding problems such as cold shortening [8, 9]. Also, the layer of fat is an important source of essential fatty acids and acts in the transport of fat-soluble vitamins, constituting a source of energy and insulation for the body of the animal [10].

Selection based on body composition, particularly on the relative proportion of muscle and fat in the carcass, is critical in meat-producing animals [5, 11]. Most carcass traits have moderate to high heritability, indicating that the selection may result in significant genetic progress [5]. According to Costa et al. [12] and Clímaco et al. [10], feedlot finished zebu breeds may present the same proportion of edible portion as other genotypes (crosses with taurine breeds), and even greater muscularity and higher carcass yield.

Several tools have been developed to improve the accuracy of animal selection and thus improve economically important traits in beef cattle, such as large-scale genotyping platforms, high-density panels of single nucleotide polymorphisms (SNP), and genome-wide association studies (GWAS). Besides these, many studies have used RNA Sequencing (RNA-Seq) to unravel complex traits in production animals. This high-throughput technology has been successfully employed in beef cattle for traits such as muscle development, intramuscular fat, and fatty acid profile, with interesting results of the phenotypic differences within and between populations [13,14,15,16,17,18,19,20].

Meat quality and carcass traits are influenced by a complex network of gene interactions in the muscle [21]. Therefore, elucidating the relationships between genes and how these genes, in turn, influence the carcass traits is critical for understanding the development of the animals, as well as the biological processes (BP) and metabolic pathways that may influence the final amount of fat and muscle in the carcasses.

Tizioto et al. [22] working with this same population of Nellore steers, identified six QTL (quantitative trait loci) that individually explained 0.8% of the additive genetic variance of REA, and a QTL that explained 0.36% of the variation in BFT. Gomes et al. [23] reported that SNPs in genes related to protein turnover, like genes regulating the ubiquitin-proteasome system, may be associated with growth and carcass traits in bovine. Junior et al. [24] found SNP-windows located on chromosomes 5, 6, 7, 8, 10, 12, 13, 16, 17, 20, 24 and 29 that together explained 8.72% and 11.38% of the genetic variance for REA and BFT in a population of Nellore cattle.

Despite those studies, there are still gaps to be filled about the molecular mechanisms that regulate carcass traits in cattle. Thus, the aim of this work was to study the Longissimus dorsi muscle transcriptome of Nellore cattle, associated with ribeye area and backfat thickness, to find differentially expressed genes, metabolic pathways, and biological processes that may regulate these traits. The results will improve our understanding of the molecular processes involved in muscle development and fat deposition of ruminants.

Results

Phenotypes and sequencing data

The phenotypic values of REA (cm2) and BFT (mm), animal identification, GEBVs (genomic estimated breeding values), the number of raw reads, and number and percentage of reads mapped against the Bos taurus UMD3.1 reference genome are shown in Tables 1 and 2. The heritability values for REA and BFT were 0.22 and 0.20, respectively [22]. There was no difference between the REA and BFT groups in regarding the intramuscular fat content, as well as the animals selected for REA were not significantly different for BFT and vice-versa (Additional file 1: Table S1). The correlation between REA and BFT in the sample of animals with contrasting GEBV (n = 22) tended to be low (r = −0.14).

Table 1 Phenotypic data, GEBV, the number of raw-reads, number of reads after cleaning, number and percentage of mapped reads for High and Low groups of ribeye area (REA)
Table 2 Phenotypic data, GEBV, the number of raw-reads, number of reads after cleaning, number and percentage of mapped reads for High and Low groups of backfat thickness (BFT)

The choice of GEBV to select animals within extreme groups was made following Meuwissen et al. [25] and Sosnicki and Newman [26], who emphasized the importance of choosing genomic values as a vehicle to incorporate molecular information into selection programs. Also, the correlation between the GEBV and REA phenotypic values was high, r = 0.93. The same occurred for BFT, with GEBV and BFT correlation value of r = 0.90.

On average 76.34% of total paired reads aligned against the reference genome. After filtering, 18,468 and 18,411 genes were used for differential expression analysis, for REA and BFT, respectively.

Differential expression analysis

Differential gene expression analysis between High and Low groups was conducted with DESeq2 software from R. DESeq2 uses statistical models based on a negative binomial distribution and is widely used to analyze RNA-Seq data since it allows more flexibility in assigning variations between samples [27]. One hundred and one differentially expressed (DE) genes were identified (false discovery rate, FDR 10%) between HighREA and LowREA groups, being 72 down-regulated and 29 up-regulated in the LowREA group. For BFT, 18 DE genes (FDR 10%) were identified, from which 13 were up-regulated and 5 were down-regulated in the LowBFT group. Figures 1 and 2 shows a Volcano plot of log2 foldChange (x-axis) vesus -log10 p value (FDR-corrected, y-axis) for REA and BFT, respectively. The gene annotation, log2foldChange, adjusted p value and p value of down- and up-regulated genes of REA and BFT can be found in Additional files 2: Table S2 and Additional file 3: Table S3, consecutively.

Fig. 1
figure1

Volcano plot of log2FoldChange (x-axis) versus –log10 p value (FDR-corrected, y-axis) of high and low genomic breeding value groups for ribeye area in Nellore steers with FDR 10%

Fig. 2
figure2

Volcano plot of log2FoldChange (x-axis) versus –log10 p value (FDR-corrected, y-axis) of high and low genomic breeding value groups for backfat thickness in Nellore steers with FDR 10%

Functional enrichment analysis

The functional enrichment analysis performed by DAVID (Database for Annotation, Visualization and Integrated Discovery) software identified two KEGG (Kyoto Encyclopedia of Genes and Genomes) pathways (p value <0.1) for REA: MAPK (Mitogen-Activated Protein Kinase) signaling pathway (bta04010) and endocytosis pathway (bta04144). The DE genes enriched for MAPK pathway were: MAX and PPM1B down-regulated, and ARRB2, PTPRR and STMN1 up-regulated in the LowREA group. For the endocytosis pathway, the enriched genes were: NEDD4 and NEDD4L down-regulated, CHMP4A and ARRB2 up-regulated in the LowREA group. In the enrichment analysis performed by BINGO (Biological Networks Gene Ontology) software, three significant biological processes (FDR 5%) were identified for REA: response to endoplasmic reticulum stress (GO: 0034976), cellular protein modification process (GO: 0006464) and macromolecule modification (GO: 0043412). These BP can be seen in Table 3. Redundant terms were not found by REVIGO (Reduce + Visualize Gene Ontology).

Table 3 Significant biological processes (FDR 5%) identified by BINGO comparing high and low genomic breeding value groups for ribeye area

For BFT, the functional enrichment analysis performed by DAVID identified five biological processes (p value <0.1) (Additional file 4: Table S4). The DE genes identified in these BP were IDO1 and ACHE, down and up-regulated in the LowBFT group, respectively. The second enrichment analysis, performed by BINGO software, identified 18 significant biological processes (FDR 5%), grouped into eight clusters of semantically similar Gene Ontology (GO) terms by REVIGO (Table 4). The DE genes identified in the BP were ACHE and SRD5A1 down-regulated, RSAD2 and RSPO3 up-regulated in the LowBFT group.

Table 4 Significant biological processes (FDR 5%) identified by BINGO comparing high and low genomic breeding value groups for backfat thickness

Discussion

Understanding how growth and development work may provide elements for increasing the profit and quality of meat production [28, 29]. Growth in livestock occurs mainly as a function of deposition of muscle and adipose tissue in the animal’s body [23]. As mentioned before, the ribeye area is a direct indicator of animal’s muscular development and has been used to predict the amount of lean meat in the carcass. In the other hand, the backfat thickness is used as an indirect indicator of meat in the carcass, and is very important to predict animal’s total body fat [3,4,5, 30, 31].

There are several biological processes involved in animal and muscle growth, such as the coordinated expression of many transcription factors (myogenic regulatory factors), genes and metabolic pathways, from the embryonic and fetal development until animals approach maturity. Most of the change in muscle weight during embryonic and fetal development is due to hyperplasia, the increase in number of muscle fibers. The postnatal stage of muscle growth (hypertrophy) consists in the increase in size of existing fibers. Both processes can be regulated by genetic factors, growth factors (insulin-like growth factors), hormones, and even environmental factors (mainly nutrition) acting as a positive or negative regulator of animal’s growth [29, 32, 33].

Furthermore, the processes of protein synthesis and degradation, also called protein turnover, affect muscle growth rates and can consequently alter carcass traits in beef cattle [23, 29, 34, 35]. During muscle hypertrophy, there is a balance between protein synthesis and degradation that may result in protein deposition, and therefore muscle growth [32]. Altogether, these processes will lead to differences in muscle and fat deposition, and hence animals with different proportions of REA and BFT.

Ribeye area

The enrichment analysis performed by DAVID identified two pathways (KEGG). The first one was MAPK pathway, which is responsible for transduction of extracellular signals to their intracellular targets in various cell types, including skeletal muscle cells. This pathway acts in the control of fundamental cellular processes, such as proliferation, growth, migration, differentiation, apoptosis, and more specifically to muscle cells, hypertrophy [36,37,38]. According to Noordman, Jansen and Hendriks [39], the MAPK pathway is the main mechanism used by growth factors in processes such as cell proliferation and differentiation.

When activated, MAPKs phosphorylate several intracellular targets, which include numerous transcription factors, resulting in the reprogramming of gene expression and cell division [40]. The activity is regulated by autophosphorylation or by phosphorylation of other kinases. On the other hand, the inactivation occurs by the process of dephosphorylation, which can be initiated by protein tyrosine phosphatases (PTPs) and metal-dependent protein phosphatases (PPMs) [39, 41, 42].

The PTPRR gene (protein tyrosine phosphatase, receptor type R), up-regulated in the LowREA group, can act by regulating the dephosphorylation of MAPKs and inhibiting cellular processes of proliferation and differentiation [43, 44]. Li et al. [42], working with PTPRR expression in mice hippocampus, verified that a greater expression of this gene led to an increase in MAPK dephosphorylation and consequently neuronal apoptosis and a decrease in cellular proliferation, showing that this gene may be acting on inhibition of the MAPK pathway in the LowREA group.

The PPM1B (protein tyrosine phosphatase, Mg2+/Mn2+ dependent 1B) gene, down-regulated in LowREA group, encodes a protein of the PPM family and acts in MAPK pathway dephosphorylation. Wei and Liang [41] identified a negative correlation between PPM1B and muscle atrophy, that is, PPM1B expression gradually decreased when muscle atrophy increased.

The second pathway found in the present study was endocytosis, which is fundamental for eukaryotic cells and is highly conserved between species and cell types. Endocytosis acts on the regulation of several processes, like cell adhesion and migration, extracellular signal transduction, cell growth and differentiation [45]. Junior et al. [24] also found genes involved in cell cycle regulation and transportation of cellular substances associated with REA in Nellore cattle.

Four genes were enriched in the endocytosis pathway: CHMP4A, ARRB2, NEDD4 and NEDD4L. Within them, NEDD4 (neural precursor cell expressed, developmentally down-regulated 4, E3 ubiquitin protein ligase) and NEDD4L (neural precursor cell expressed, developmentally down-regulated 4-like, E3 ubiquitin protein ligase), also known by NEDD4–2, encode ubiquitin proteins ligases belonging to the Nedd4 family. Among their functions, they may aid protein internalization in the cells [46,47,48].

In addition to the protein internalization function, NEDD4 is required for cell surface expression of the IGF-1R (insulin-like growth factor, type 1 receptor) and insulin receptor, and is a positive regulator of IGF-1 (insulin-like growth factor, type 1) and insulin signaling [46, 47]. In mammals, the Insulin-like growth factors (IGF) axis is the largest fetal and postnatal growth regulator and is strongly related to muscle differentiation [32, 49, 50]. Studies with knockout mice for the NEDD4 gene showed that loss of NEDD4 reduced IGF-1 and insulin signaling, delayed embryonic development, and reduced growth and body weight [46]. Junior et al. [24] in a GWAS study associated with REA, BFT and hot carcass weight found this gene enriched for the GO terms “cellular protein metabolic process” and “protein metabolic process” related to protein turnover and, consequently animal growth and development. In the present study, NEDD4 and NEDD4L were down-regulated in LowREA group, emphasizing their importance in regulating muscle growth.

ARRB2 (arrestin β-2), identified in both pathways – endocytosis and MAPK –, was up-regulated in the LowREA group. β-arrestins are multifunctional signaling molecules ubiquitously expressed that act as endocytosis regulators in different types of cell surface receptors [51, 52]. According to Luttrell and Lefkowitz [51], β-arrestins can serve as scaffold proteins for MAPK pathway proteins. Additionally, Yan et al. [53] show the involvement of β-arrestin 2 in the activation of MAPK pathway.

Analysis with BINGO software ascertained three biological processes: response to endoplasmic reticulum stress (GO: 0034976), cellular protein modification process (GO: 0006464) and macromolecule modification (GO: 0043412) (Table 3). Among the genes identified in these BP, AMFR (autocrine motility factor), a down-regulated gene in LowREA group, appears in all of them. AMFR – also known as gp78 – encodes a RING (Really Interesting New Gene) class E3 ubiquitin protein ligase that is involved in the mechanism of protein quality control, eliminating misfolded proteins from the endoplasmic reticulum of eukaryotic cells [54,55,56].

The endoplasmic reticulum (ER) is a ubiquitous multifunctional organelle, which ensures the correct protein formation, and plays a key role in lipids and sterols synthesis, and in intracellular calcium maintenance [57]. ER stress can occur due to perturbations in its homeostasis, such as chemical damage, gene mutations, nutrient insufficiency, cell differentiation, oxidative stress, and fluctuation in calcium concentrations, leading to changes in protein structure, resulting in the accumulation of misfolded proteins in the ER lumen [58, 59]. ER stress can alter gene expression and cause post-transcriptional modifications, change cell physiology, and even induce cell apoptosis [57, 60, 61].

According to Nakanishi, Sudo and Morishima [61] and Nakanishi, Dohmae and Morishima [62] the ER response to stress related to the induction of apoptosis may be favorable for myogenesis. Nakanishi, Dohmae, and Morishima [62] working with mouse myoblast cells (C2C12) demonstrated that apoptosis induced by ER stress controls the differentiation of myoblasts, so only cells resistant to apoptosis undergo terminal differentiation to muscle tissue formation, improving the myoblasts quality.

The other two biological processes identified in the present study - cellular protein and macromolecule modification processes - are intrinsically related, since proteins can be classified as macromolecules. According to Cantin and Yates III [63], most proteins need to undergo modifications to carry out their activities or become biologically active, and these changes are called post-translational modifications (PTM).

PTMs are chemical changes that modify the protein structure reversibly or irreversibly, through proteolytic cleavage or covalent modifications in specific amino acid residues [64,65,66]. According to Blom et al. [65] and Zou and Blank [67], phosphorylation is the primary protein modifier and is considered a key event in several transductional signaling cascades, such as the MAPK pathway. As discussed previously, PTPRR and PPM1B, identified in the MAPK pathway, were also identified in the macromolecular modification process. These two genes encode protein phosphatases that can act by dephosphorylating and thus decreasing the activity of MAPK proteins [39, 41, 42].

Another gene found in the macromolecular modification process that encode a protein phosphatase is EYA2 (EYA transcriptional coactivator and phosphatase 2). Unlike PTPRR and PPM1B, this gene acts as transcription factor inducing myogenic regulatory factors, such MEF3 and MYOG, as well has important roles in differentiated muscle cells [68, 69]. In our finds, EYA2 was down-regulated in LowREA group, showing its importance as a positive regulator of muscle growth.

Another common PTM is ubiquitination, which acts in directing short-life proteins to the proteasome degradation pathway [34, 54,55,56]. This ubiquitin-dependent proteolysis ensures protein turnover that is essential to cell survival [70]. In addition, the ubiquitination process also functions on cellular processes like signal transduction, enzymatic activation, endocytosis, molecular trafficking, chromatin rearrangement and DNA repair [71]. Among the genes identified in the cellular protein modification and macromolecular modification processes that encode ubiquitin proteins, NEDD4 and NEDD4L already had their functions in muscle discussed here. The UBE4A (ubiquitin conjugation factor E4 A), down-regulated in the LowREA group, is an important and ubiquitously expressed gene for the ubiquitination process, that may also participate in growth and differentiation processes [72]. Gomes et al. [23] have already reported that genes related to protein turnover were associated with growth and carcass traits in Nellore cattle. Latterly, Junior et al. [24] also found GO terms related to protein turnover in an association study with REA, BFT and hot carcass weight in Nellore. These findings show us the important role of ubiquitination and ubiquitin proteins for muscle growth, and consequently to improve REA in the animals.

Backfat thickness

The enrichment analysis performed by BINGO software identified 18 biological processes (FDR 5%), which were grouped into eight clusters of semantically similar GO terms by REVIGO (Table 4).

Acetylcholinesterase gene (ACHE), present in four similarity clusters and up-regulated in the LowBFT group, is an essential component of the neuromuscular junction (NMJ). This enzyme is highly conserved in mammals and appears in multiple molecular forms, which originate from the alternative splicing of ACHE gene [73, 74]. ACHE expression in muscle is regulated by transcriptional and post-transcriptional events, with an increased expression in the early stages of myogenic differentiation, and reaching a plateau when the myotubes are mature [74, 75].

Mouisel et al. [76] reported a loss of muscular weight in the hind limb muscles of ACHE knockout mice. Soysal et al. [77] concluded that acetylcholinesterase inhibitors could cause weight loss and change muscle mass index in elderly people. Despite not showing a direct relation with BFT, it is clear that this gene has a role during animal’s muscle growth.

The SRD5A1 (steroid-5-alpha-reductase, alpha polypeptide 1), also known as 5α-reductase type 1, was identified in the androgen biosynthetic process (GO:0006702). Androgens, such as testosterone, play a critical role in muscle, increasing protein synthesis and energy metabolism, and promoting growth and muscle strength increase [78, 79]. Ferrando et al. [79] hypothesized that testosterone might stimulate IGF-1 release in muscle tissue. Although skeletal muscle can synthesize and metabolize testosterone, this action on target organs often requires its metabolic conversion to one or more active products [80], such as DHT (dihydrotestosterone) metabolized by the 5α-reductase enzyme. DHT is one of the most potent natural androgens because of its high affinity for androgen receptors; it has several physiological effects on skeletal muscle, like activation of signaling pathways and anabolic action in protein synthesis, as well as the maintenance of muscle homeostasis [81,82,83,84].

Several studies found an association of SRD5A1 and its corresponding protein (5α-reductase) with muscle weight and strength [83,84,85,86]. Sato et al. [86] identified a positive correlation of 5α-reductase protein with the cross-sectional area of quadriceps femoralis muscle in humans. Although not DE for REA, SRD5A1 was up-regulated in the LowBFT group, that is, the group that presented a higher proportion of muscle mass represented by higher values of REA (Table 2). In contrast, Sun et al. [87] studying putative target genes for miRn25 and n26, highly expressed miRNAs in bovine backfat thickness, identified SRD5A1 related to lipid synthesis in adult animals.

The R-spondin 3 gene (RSPO3), down-regulated in the LowBFT group, encodes a member of a protein family widely recognized as an agonist of the canonical Wnt signaling pathway (or Wnt/β-catenin pathway) [88,89,90,91], one of the biological processes found here enriched for this gene (GO:0060070). This pathway plays an essential role during embryonic muscle development and in skeletal muscle homeostasis during adulthood [89, 90, 92]. This pathway also is an important regulator of adipocyte differentiation [93, 94].

Han et al. [88] studied the role of RSPOs (r-spondins) during myogenic differentiation using primary satellite cells and C2C12 cells from mouse myoblast. The authors observed that silencing RSPO2 and RSPO3 significantly affected the Myf5 expression, the rate of myogenic differentiation and the myotubes formation in mice muscle cells. The authors also found that RSPOs can act via canonical Wnt signaling pathway in the positive regulation of myogenesis in skeletal muscle.

Li et al. [93] studying adipose-derived mesenchymal stem cells in porcine found that the activation of Wnt signaling pathway suppressed mRNA and protein expression of the adipocyte-specific genes C/EBPα (CCAAT/enhancer-binding protein-α) and PPARγ (peroxisome proliferator-activated receptor-γ), inhibiting adipogenesis in these cells. Chen et al. [94] also found that Wnt signaling pathway may inhibit adipogenic differentiation in porcine intramuscular preadipocytes. So, even identifying this gene down-regulated in the LowBFT group, it is likely that it was acting as a negative regulator of BFT in the animals.

RSAD2 gene (radical domain of S-adenosyl methionine containing 2), found in the defense response to virus process (GO:0051607), is a type I interferon response gene, which has been used in the clinical prediction of some diseases, also related to skeletal muscle myopathies caused by inflammatory cytokines [95,96,97]. Wei et al. [98] working with growing pigs fed a linseed-enriched diet, found RSAD2 up-regulated in the Longissimus dorsi muscle of treated animals. Dogan et al. [99], studying the structure and composition of mice fat, muscle, and liver, reported that RSAD2 may act as a modulator of lipid droplet content and lipid biosynthesis in adipose tissue. These findings coincide with this work, in which the RSAD2 gene was down-regulated in the LowBFT group.

Conclusions

Our results emphasize the complexity of gene regulation in the Longissimus dorsi muscle of Nellore cattle associated with REA and BFT. We identified 101 DE genes in the extreme GEBV groups for REA. These genes were enriched for metabolic pathways and biological processes mostly involved in differentiation, the proliferation of muscle cells, protein turnover and hypertrophy, such as the MAPK pathway, cellular protein, and macromolecule modification processes. For BFT, we identified 18 DE genes involved in biological processes that may regulate positively or negatively adipogenesis, lipid biosynthesis and muscle growth. These results might help us to enlighten the molecular processes involved in muscle and fat deposition, which are economically important carcass traits for beef production.

Methods

Animals, samples, and phenotypes

Three hundred eighty-five (385) Nellore steers from Embrapa (Brazilian Agricultural Research Corporation) breeding herd, raised between 2009 and 2011, were included in this study. To breed this herd, 34 unrelated bulls were selected representing the main breeding lineages used in Brazil based on the information of the National Summary of Nellore produced by the Brazilian Association of Zebu Breeders (ABCZ) and the National Research Center for Beef Cattle.

The animals were raised in grazing systems, under the same conditions of handling and nutrition until 21 months of age when they were taken to feedlots. All animals were slaughtered at an average age of 25 months. The slaughter was carried out in a commercial slaughterhouse located in the city of Bariri (São Paulo), under the supervision of the Federal Inspection Service (SIF) and within the standards established by the Brazilian Ministry of Agriculture, Livestock and Food Supply (MAPA), for more details see [100]. At the time of slaughter, approximately 5 g of the Longissimus dorsi muscle was collected between the 12th and 13th ribs (right half carcass) and were stored in liquid nitrogen. Twenty-four hour after slaughtering, steaks corresponding to a cross section of the Longissimus dorsi muscle between the 12th and 13th ribs (left half carcass) were sampled with bone and transported to the laboratory of Embrapa Pecuária Sudeste (São Carlos, SP), where REA and BFT were measured. REA was measured with a grid and the BFT with a graduated ruler.

The genomic estimated breeding values (GEBV) were obtained by the GenSel program [101], which uses Bayesian methodology. The a priori values of genetic and residual variance were obtained from the Bayes C analysis, in which the a priori genetic and residual variance was equal to 1 [102]. Using the estimated a priori values, a new Bayes C analysis was performed to obtain GEBVs for each animal. The SNP markers information was obtained as described by Cesar et al. [103] using BovineHD 770 k BeadChip (Infinium BeadChip, Illumina, San Diego, CA, USA). For BFT, 384 animals were used in the GEBV estimate. The animals were separated into two groups of six animals each (High and Low), based on the extreme values of GEBVs for each of the two traits. Of the 12 animals selected for each trait, two of them were in common for both traits (Additional file 5: Table S5). A Student’s t-test was performed to verify the difference in REA and BFT level between High and Low groups. The phenotypic values of intramuscular fat content were also included in this test to ascertain these animals were not significantly different for this trait. The phenotypic correlation between REA and BFT using the selected animals (n = 22) was estimated with the R software.

RNA extraction, quality analysis, library preparation and sequencing

The total RNA was extracted from 100 mg of frozen Longissimus dorsi muscle collected at the slaughter using the TRIzol reagent (Life Technologies, Carlsbad, CA, USA), following the manufacturer’s instructions. At the end of the extraction process, RNA integrity was verified by the Bioanalyzer 2100 (Agilent, Santa Clara, CA, USA). The mean RIN (RNA integrity number) of all samples was 7.75. For library preparation, 2 μg of RNA from each sample was used, according to the protocol described in the TruSeq RNA Sample Preparation kit v2 guide (Illumina, San Diego, CA, USA). The libraries were quantified by quantitative PCR using the KAPA Library Quantification kit (KAPA Biosystems, Foster City, CA, USA) and the average library size was estimated using the Bioanalyzer 2100 (Agilent, Santa Clara, CA, USA). After quantification, the samples were diluted and pooled into three pools out of six samples each. Three lanes of a sequencing flowcell were clustered, using the TruSeq PE Cluster kit v3-cBot-HS (Illumina, San Diego, CA, USA). They were sequenced on the HiSeq2500 ultra-high throughput sequencing system (Illumina, San Diego, CA, USA) using the TruSeq SBS kit v3-HS (200 cycles), according to the manufacturer’s instructions. All sequencing analyses were performed at the Genomics Center at the College of Agriculture “Luiz de Queiroz” of the University of São Paulo.

Quality control and read alignment

The adapter sequences and low complexity reads were removed in an initial data-filtering step using SeqClean software (https://sourceforge.net/projects/seqclean/files/). The FastQC software was used to analyze the quality of raw reads (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/). The Tophat version 2.1.0 software [104] was used to map the reads against to the UMD3.1 Bos taurus reference genome (http://www.ensembl.org/Bos_taurus/Info/Index/). Read counts (mRNA abundance) for all mapped genes were calculated using the HTSeq version 0.6.1 software (http://www-huber.embl.de/HTSeq) [105]. Only read sequences that uniquely mapped to known chromosomes were used in this study.

Identification and annotation of differentially expressed genes

Differentially expressed genes were identified using the DESeq2 software of R [106]. Before the statistical analysis was performed, the read count data was filtered based on previous studies [16, 17], as follows: i) genes with zero counts were removed (not expressed); ii) genes with less than one read per sample on average were removed (very low expression level); iii) genes that did not appear in at least three samples were removed (rarely expressed). After filtering, a total of 18,468 genes for REA and 18,411 for BFT were analyzed for differential expression employing the “nbinomWaldTest” function of the DESeq2 that assumes the level of gene expression as a negative binomial distribution. Exploratory plots were made to check the dispersion estimates (Additional file 6: Figures S1 and S2, and Additional file 7: Figures S3 and S4). The Benjamini and Hochberg [107] methodology was used to control the false discovery rate (FDR) at 10%. The DE genes were annotated by the online tool BioMart (http://www.ensembl.org/biomart) from Ensembl. Genes that lacked annotation information were annotated using the NCBI (National Center for Biotechnology Information - http://www.ncbi.nlm.nih.gov/) and Panther (http: /www.pantherdb.org) databases.

Functional enrichment analysis

The functional enrichment analysis of DE genes (FDR 10%) for KEGG (Kyoto Encyclopedia of Genes and Genomes) pathways was carried out with the online tool DAVID (Database for Annotation, Visualization and Integrated Discovery) version 6.7 [108].

Also, another enrichment analysis was performed to identify biological processes related to the DE genes, using BINGO (Biological Networks Gene Ontology) version 3.0.3 [109], a Cytoscape [110] version 3.4.0 app. BINGO is a free-use tool that determines GO terms that are over-represented in a set of genes using the “Hypergeometric test” as a statistical test. BPs that presents FDR 5% [107] were considered significant. Lastly, REVIGO (Reduce + Visualize Gene Ontology), an algorithm that summarizes long lists of GO terms was used to remove redundant GO terms. REVIGO performs a simple clustering procedure finding a representative subset of GO terms that is based on semantic similarity measures [111].

Abbreviations

BFT:

Backfat thickness

BINGO:

Biological Networks Gene Ontology

BP:

Biological process

DAVID:

Database for Annotation, Visualization and Integrated Discovery

DE:

Differentially expressed

DHT:

Dihydrotestosterone

ER:

Endoplasmic reticulum

FDR:

False discovery rate

GEBV:

Genomic estimated breeding value

GO:

Gene ontology

GWAS:

Genome-wide association study

IGF:

Insulin-like growth factor

KEGG:

Kyoto Encyclopedia of Genes and Genomes

MAPK:

Mitogen-Activated Protein Kinase

NMJ:

Neuromuscular junction.

PPM:

Metal-dependent protein phosphatase.

PTM:

Post-translational modification.

PTP:

Protein tyrosine phosphatase.

QTL:

Quantitative trait loci.

REA:

Ribeye area.

REVIGO:

Reduce + Visualize Gene Ontology.

RSPOs:

r-spondins.

SNP:

Single nucleotide polymorphism.

References

  1. 1.

    DMd LJ, AHdN R, Urbano SA, MdV M, LPdA A. Alguns aspectos qualitativos da carne bovina: uma revisão. Acta Vet Bras. 2011;5(4):351–8.

  2. 2.

    Alves DD, Mancio AB. Maciez da carne bovina - Uma revisão. Revista da FZVA. 2007:14(1).

  3. 3.

    Bergen R, McKinnon JJ, Christensen D, Kohle N. Prediction of lean yield in yearling bulls using real-time ultrasound. Can Vet J. 1996;76(3):305–11.

  4. 4.

    Moser DW, Bertrand JK, Misztal I, Kriese LA, Benyshek LL. Genetic parameter estimates for carcass and yearling ultrasound measurements in Brangus cattle. J Anim Sci. 1998;76(10):2542–8.

  5. 5.

    Burrow HM, Moore SS, Johnston DJ, Barendse W, Bindon BM. Quantitative and molecular genetic influences on properties of beef: a review. Aust J Exp Agric. 2001;41(7):893–919.

  6. 6.

    Bianchini W, Silveira AC, Jorge AM, Arrigoni MDB, Martins CL, Rodrigues É, et al. Effect of genetic group on carcass traits and fresh and aged beef tenderness from young cattle. R Bras Zootec. 2007;36(6):2109–17.

  7. 7.

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

  8. 8.

    Wilson DE. Application of ultrasound for genetic improvement. J Anim Sci. 1992;70(3):973–83.

  9. 9.

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

  10. 10.

    Clímaco SM, ELdA R, Mizubuti IY, LdDFd S, MAAdF B, Bridi AM. Desempenho e características de carcaça de bovino de corte de quatro grupos genéticos terminados em confinamento. Rev Bras Zootec. 2011;40(7):1562–7.

  11. 11.

    Andersson L. Genetic dissection of phenotypic diversity in farm animals. Nat Rev Genet. 2001;2(2):130–8.

  12. 12.

    Costa D, Abreu JBR, RdC M, JCGd S, Rodrigues VC, JCDd S, et al. Características de Carcaça de Novilhos Inteiros Nelore e F1 Nelore x Holandês. Ciência Animal Bras. 2007;8(4):685–94.

  13. 13.

    Baldwin RL, Li RW, Li CJ, Thomson JM, Bequette BJ. Characterization of the longissimus lumborum transcriptome response to adding propionate to the diet of growing Angus beef steers. Physiol Genomics. 2012;44(10):543–50.

  14. 14.

    Djari A, Esquerré D, Weiss B, Martins F, Meersseman C, Boussaha M, et al. Gene-based single nucleotide polymorphism discovery in bovine muscle using next-generation transcriptomic sequencing. BMC Genomics. 2013;14:307.

  15. 15.

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

  16. 16.

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

  17. 17.

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

  18. 18.

    Tizioto PC, Coutinho LL, Decker JE, Schnabel RD, Rosa KO, Oliveira PS, et al. Global liver gene expression differences in Nelore steers with divergent residual feed intake phenotypes. BMC Genomics. 2015;16:242.

  19. 19.

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

  20. 20.

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

  21. 21.

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

  22. 22.

    Tizioto PC, Decker JE, Taylor JF, Schnabel RD, Mudadu MA, Silva FL, et al. Genome scan for meat quality traits in Nelore beef cattle. Physiol Genomics. 2013;45(21):1012–20.

  23. 23.

    Gomes RC, Silva SL, Carvalho ME, Rezende FM, Pinto LF, Santana MH, et al. Protein synthesis and degradation gene SNPs related to feed intake, feed efficiency, growth, and ultrasound carcass traits in Nellore cattle. Genet Mol Res. 2013;12(3):2923–36.

  24. 24.

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

  25. 25.

    Meuwissen TH, Hayes BJ, Goddard ME. Prediction of total genetic value using genome-wide dense marker maps. Genetics. 2001;157(4):1819–29.

  26. 26.

    Sosnicki AA, Newman S. The support of meat value chains by genetic technologies. Meat Sci. 2010;86(1):129–37.

  27. 27.

    Ching T, Huang S, Garmire LX. Power analysis and sample size estimation for RNA-Seq differential expression. RNA. 2014;20(11):1684–96.

  28. 28.

    Rehfeldt C, Fieldler I, Dietl G, Ender K. Myogenesis and postnatal skeletal muscle cell growth as influenced by selection. Livest Prod Sci. 2000;66(2):177–88.

  29. 29.

    Maltin C, Balcerzak D, Tilley R, Delday M. Determinants of meat quality: tenderness. Proc Nutr Soc. 2003;62(2):337–47.

  30. 30.

    Yokoo MJ, Albuquerque LG, Lôbo RB, Bezerra LAF, Araujo FRC, Silva JAV, et al. Genetic and environmental factors affecting ultrasound measures of longissimus muscle area and backfat thickness in Nelore cattle. Livest Sci. 2008;117(Issues 2–3):147–54.

  31. 31.

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

  32. 32.

    Oksbjerg N, Gondret F, Vestergaard M. Basic principles of muscle development and growth in meat-producing mammals as affected by the insulin-like growth factor (IGF) system. Domest Anim Endocrinol. 2004;27(3):219–40.

  33. 33.

    Aberle ED, Forrest JC, Gerrard DE, Mills EW. Principles of meat science. 5th ed. Dubuque, Iowa, USA: Kendall Hunt; 2012.

  34. 34.

    Fanzani A, Conraads VM, Penna F, Martinet W. Molecular and cellular mechanisms of skeletal muscle atrophy: an update. J Cachexia Sarcopenia Muscle. 2012;3(3):163–79.

  35. 35.

    Gumucio JP, Mendias CL. Atrogin-1, MuRF-1, and sarcopenia. Endocrine. 2013;43(1):12–21.

  36. 36.

    Martineau LC, Gardiner PF. Insight into skeletal muscle mechanotransduction: MAPK activation is quantitatively related to tension. J Appl Physiol (1985). 2001;91(2):693–702.

  37. 37.

    Meister M, Tomasovic A, Banning A, Tikkanen R. Mitogen-activated protein (MAP) Kinase scaffolding proteins: a recount. Int J Mol Sci. 2013;14(3):4854–84.

  38. 38.

    Chakraborty C, Sharma AR, Patra BC, Bhattacharya M, Sharma G, Lee SS. MicroRNAs mediated regulation of MAPK signaling pathways in chronic myeloid leukemia. Oncotarget. 2016;7(27):42683–97.

  39. 39.

    Noordman YE, Jansen PA, Hendriks WJ. Tyrosine-specific MAPK phosphatases and the control of ERK signaling in PC12 cells. J Mol Signal. 2006;1:4.

  40. 40.

    Liang Q, Molkentin JD. Redefining the roles of p38 and JNK signaling in cardiac hypertrophy: dichotomy between cultured myocytes and animal models. J Mol Cell Cardiol. 2003;35(12):1385–94.

  41. 41.

    Wei J, Liang B-S. PPM1B and P-IKKβ expression levels correlated inversely with rat gastrocnemius atrophy after denervation. Braz J Med Biol Res. 2012;45(8):711–5.

  42. 42.

    Li L, Gao G, Shankar J, Joshi B, Foster LJ, Nabi IR. p38 MAP kinase-dependent phosphorylation of the Gp78 E3 ubiquitin ligase controls ER-mitochondria association and mitochondria motility. Mol Biol Cell. 2015;26(21):3828–40.

  43. 43.

    Tabernero L, Aricescu AR, Jones EY, Szedlacsek SE. Protein tyrosine phosphatases: structure-function relationships. FEBS J. 2008;275(5):867–82.

  44. 44.

    Su PH, Lin YW, Huang RL, Liao YP, Lee HY, Wang HC, et al. Epigenetic silencing of PTPRR activates MAPK signaling, promotes metastasis and serves as a biomarker of invasive cervical cancer. Oncogene. 2013;32(1):15–26.

  45. 45.

    Doherty GJ, McMahon HT. Mechanisms of endocytosis. Annu Rev Biochem. 2009;78:857–902.

  46. 46.

    Cao XR, Lill NL, Boase N, Shi PP, Croucher DR, Shan H, et al. Nedd4 controls animal growth by regulating IGF-1 signaling. Sci Signal. 2008;1(38):ra5.

  47. 47.

    Monami G, Emiliozzi V, Morrione A. Grb10/Nedd4-mediated multiubiquitination of the insulin-like growth factor receptor regulates receptor internalization. J Cell Physiol. 2008;216(2):426–37.

  48. 48.

    Goel P, Manning JA, Kumar S. NEDD4-2 (NEDD4L): the ubiquitin ligase for multiple membrane proteins. Gene. 2015;557(1):1–10.

  49. 49.

    Barker J, Liu J-P, Robertson EJ, Efstratiadis A. Role of insulin-like growth factors in embryonic and postnatal growth. Cell. 1993;75(1):73–82.

  50. 50.

    Wu Z, Woodring PJ, Bhakta KS, Tamura K, Wen F, Feramisco JR, et al. p38 and extracellular signal-regulated kinases regulate the myogenic program at multiple steps. Mol Cell Biol. 2000;20(11):3951–64.

  51. 51.

    Luttrell LM, Lefkowitz RJ. The role of beta-arrestins in the termination and transduction of G-protein-coupled receptor signals. J Cell Sci. 2002;115(Pt 3):455–65.

  52. 52.

    Yu MC, Su LL, Zou L, Liu Y, Wu N, Kong L, et al. An essential function for beta-arrestin 2 in the inhibitory signaling of natural killer cells. Nat Immunol. 2008;9(8):898–907.

  53. 53.

    Yan H, Li H, Denney J, Daniels C, Singh K, Chua B, et al. β-arrestin 2 attenuates cardiac dysfunction in polymicrobial sepsis through gp130 and p38. Biochem Biophys Rep. 2016;7:130–7.

  54. 54.

    Fang S, Ferrone M, Yang C, Jensen JP, Tiwari S, Weissman AM. The tumor autocrine motility factor receptor, gp78, is a ubiquitin protein ligase implicated in degradation from the endoplasmic reticulum. Proc Natl Acad Sci U S A. 2001;98(25):14422–7.

  55. 55.

    Chen Z, Ballar P, Fu Y, Luo J, Du S, Fang S. The E3 ubiquitin ligase gp78 protects against ER stress in zebrafish liver. J Genet Genomics. 2014;41(7):357–68.

  56. 56.

    Zhang T, Xu Y, Liu Y, Ye Y. gp78 functions downstream of Hrd1 to promote degradation of misfolded proteins of the endoplasmic reticulum. Mol Biol Cell. 2015;26(24):4438–50.

  57. 57.

    Kaufman RJ. Stress signaling from the lumen of the endoplasmic reticulum: coordination of gene transcriptional and translational controls. Genes Dev. 1999;13(10):1211–33.

  58. 58.

    Cimellaro A, Perticone M, Fiorentino TV, Sciacqua A, Hribal ML. Role of endoplasmic reticulum stress in endothelial dysfunction. Nutr Metab Cardiovasc Dis. 2016;26(10):863–71.

  59. 59.

    Nakanishi K, Kakiguchi K, Yonemura S, Nakano A, Morishima N. Transient Ca2+ depletion from the endoplasmic reticulum is critical for skeletal myoblast differentiation. FASEB J. 2015;29(5):2137–49.

  60. 60.

    Ron D: Translational control in the endoplasmic reticulum stress response. In: J Clin Invest vol 110; 2002: 1383-1388.

  61. 61.

    Nakanishi K, Sudo T, Morishima N. Endoplasmic reticulum stress signaling transmitted by ATF6 mediates apoptosis during muscle development. J Cell Biol. 2005;169(4):555–60.

  62. 62.

    Nakanishi K, Dohmae N, Morishima N. Endoplasmic reticulum stress increases myofiber formation in vitro. FASEB J. 2007;21(11):2994–3003.

  63. 63.

    Cantin GT, Yates JR 3rd. Strategies for shotgun identification of post-translational modifications by mass spectrometry. J Chromatogr A. 2004;1053(1–2):7–14.

  64. 64.

    Mann M, Jensen ON. Proteomic analysis of post-translational modifications. Nat Biotechnol. 2003;21(3):255–61.

  65. 65.

    Blom N, Sicheritz-Ponten T, Gupta R, Gammeltoft S, Brunak S. Prediction of post-translational glycosylation and phosphorylation of proteins from the amino acid sequence. Proteomics. 2004;4(6):1633–49.

  66. 66.

    Jensen ON. Modification-specific proteomics: characterization of post-translational modifications by mass spectrometry. Curr Opin Chem Biol. 2004;8(1):33–41.

  67. 67.

    Zou X, Blank M. Targeting p38 MAP kinase signaling in cancer through post-translational modifications. Cancer Lett. 2017;384:19–26.

  68. 68.

    Ridgeway AG, Skerjanc IS. Pax3 is essential for skeletal myogenesis and the expression of Six1 and Eya2. J Biol Chem. 2001;276(22):19033–9.

  69. 69.

    Hudson NJ, Lyons RE, Reverter A, Greenwood PL, Dalrymple BP. Inferring the in vivo cellular program of developing bovine skeletal muscle from expression data. Gene Expr Patterns. 2013;13(3–4):109–25.

  70. 70.

    Ciechanover A, Orian A, Schwartz AL. The ubiquitin-mediated proteolytic pathway: mode of action and clinical implications. J Cell Biochem Suppl. 2000;34:40–51.

  71. 71.

    Nguyen LK, Kolch W, Kholodenko BN. When ubiquitination meets phosphorylation: a systems biology perspective of EGFR/MAPK signalling. Cell Commun Signal. 2013;11:52.

  72. 72.

    Contino G, Amati F, Pucci S, Pontieri E, Pichiorri F, Novelli A, et al. Expression analysis of the gene encoding for the U-box-type ubiquitin ligase UBE4A in human tissues. Gene. 2004;328:69–74.

  73. 73.

    Massoulie J. The origin of the molecular diversity and functional anchoring of cholinesterases. Neurosignals. 2002;11(3):130–43.

  74. 74.

    Tsim KW, Choi RC, Xie HQ, Zhu JT, Leung KW, Lau FT, et al. Transcriptional control of different subunits of AChE in muscles: signals triggered by the motor nerve-derived factors. Chem Biol Interact. 2008;175(1–3):58–63.

  75. 75.

    Angus LM, Chan RY, Jasmin BJ. Role of intronic E- and N-box motifs in the transcriptional induction of the acetylcholinesterase gene during myogenic differentiation. J Biol Chem. 2001;276(20):17603–9.

  76. 76.

    Mouisel E, Blondet B, Escourrou P, Chatonnet A, Molgo J, Ferry A. Outcome of acetylcholinesterase deficiency for neuromuscular functioning. Neurosci Res. 2006;55(4):389–96.

  77. 77.

    Soysal P, Isik AT, Stubbs B, Solmi M, Volpe M, Luchini C, et al. Acetylcholinesterase inhibitors are associated with weight loss in older people with dementia: a systematic review and meta-analysis. J Neurol Neurosurg Psychiatry. 2016;87(12):1368–74.

  78. 78.

    Ramamani A, Aruldhas MM, Govindarajulu P. Impact of testosterone and oestradiol on region specificity of skeletal muscle-ATP, creatine phosphokinase and myokinase in male and female Wistar rats. Acta Physiol Scand. 1999;166(2):91–7.

  79. 79.

    Ferrando AA, Sheffield-Moore M, Yeckel CW, Gilkison C, Jiang J, Achacosa A, et al. Testosterone administration to older men improves muscle function: molecular and physiological mechanisms. Am J Physiol Endocrinol Metab. 2002;282(3):E601–7.

  80. 80.

    Fuxjager MJ, Schuppe ER, Hoang J, Chew J, Shah M, Schlinger BA. Expression of 5alpha- and 5beta-reductase in spinal cord and muscle of birds with different courtship repertoires. Front Zool. 2016;13:25.

  81. 81.

    Yoshioka M, Boivin A, Ye P, Labrie F, St-Amand J. Effects of dihydrotestosterone on skeletal muscle transcriptome in mice measured by serial analysis of gene expression. J Mol Endocrinol. 2006;36(2):247–59.

  82. 82.

    Sato K, Iemitsu M, Aizawa K, Ajisaka R. Testosterone and DHEA activate the glucose metabolism-related signaling pathway in skeletal muscle. Am J Physiol Endocrinol Metab. 2008;294(5):E961–8.

  83. 83.

    Aizawa K, Iemitsu M, Maeda S, Mesaki N, Ushida T, Akimoto T. Endurance exercise training enhances local sex steroidogenesis in skeletal muscle. Med Sci Sports Exerc. 2011;43(11):2072–80.

  84. 84.

    Windahl SH, Andersson N, Borjesson AE, Swanson C, Svensson J, Moverare-Skrtic S, et al. Reduced bone mass and muscle strength in male 5alpha-reductase type 1 inactivated mice. PLoS One. 2011;6(6):e21402.

  85. 85.

    Aizawa K, Iemitsu M, Maeda S, Otsuki T, Sato K, Ushida T, et al. Acute exercise activates local bioactive androgen metabolism in skeletal muscle. Steroids. 2010;75(3):219–23.

  86. 86.

    Sato K, Iemitsu M, Matsutani K, Kurihara T, Hamaoka T, Fujita S. Resistance training restores muscle sex steroid hormone steroidogenesis in older men. FASEB J. 2014;28(4):1891–7.

  87. 87.

    Sun J, Zhou Y, Cai H, Lan X, Lei C, Zhao X, et al. Discovery of novel and differentially expressed MicroRNAs between fetal and adult Backfat in cattle. PLoS One. 2014:9(2).

  88. 88.

    Han XH, Jin YR, Seto M, Yoon JK. A WNT/beta-catenin signaling activator, R-spondin, plays positive regulatory roles during skeletal myogenesis. J Biol Chem. 2011;286(12):10649–59.

  89. 89.

    von Maltzahn J, Chang NC, Bentzinger CF, Rudnicki MA. Wnt signaling in myogenesis. Trends Cell Biol. 2012;22(11):602–9.

  90. 90.

    Yoon JK, Lee JS. Cellular signaling and biological functions of R-spondins. Cell Signal. 2012;24(2):369–77.

  91. 91.

    Shi GX, Mao WW, Zheng XF, Jiang LS. The role of R-spondins and their receptors in bone metabolism. Prog Biophys Mol Biol. 2016;122(2):93–100.

  92. 92.

    Buckingham M, Montarras D. Skeletal muscle stem cells. Curr Opin Genet Dev. 2008;18(4):330–6.

  93. 93.

    Li HX, Luo X, Liu RX, Yang YJ, Yang GS. Roles of Wnt/beta-catenin signaling in adipogenic differentiation potential of adipose-derived mesenchymal stem cells. Mol Cell Endocrinol. 2008;291(1–2):116–24.

  94. 94.

    Chen X, Luo Y, Jia G, Liu G, Zhao H, Huang Z. The effect of arginine on the Wnt/beta-catenin signaling pathway during porcine intramuscular preadipocyte differentiation. Food Funct. 2017;8(1):381–6.

  95. 95.

    Filkova M, Senolt L, Vencovsky J. The role of resistin in inflammatory myopathies. Curr Rheumatol Rep. 2013;15(6):336.

  96. 96.

    He P, Zhang Z, Liao W, Xu D, Fu M, Kang Y. Screening of gene signatures for rheumatoid arthritis and osteoarthritis based on bioinformatics analysis. Mol Med Rep. 2016;14:1587–93.

  97. 97.

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

  98. 98.

    Wei H, Zhou Y, Jiang S, Huang F, Peng J. Transcriptional response of porcine skeletal muscle to feeding a linseed-enriched diet to growing pigs. J Anim Sci Biotechnol. 2016;7:6.

  99. 99.

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

  100. 100.

    Tizioto PC, Meirelles SL, Veneroni GB, Tullio RR, Rosa AN, Alencar MM, et al. A SNP in ASAP1 gene is associated with meat quality and production traits in Nelore breed. Meat Sci. 2012;92(4):855–7.

  101. 101.

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

  102. 102.

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

  103. 103.

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

  104. 104.

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

  105. 105.

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

  106. 106.

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

  107. 107.

    Benjamini Y, Hochberg Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J R Stat Soc Ser B Methodol. 1995;57(1):289–300.

  108. 108.

    Huang DW, Sherman BT, Tan Q, Kir J, Liu D, Bryant D, et al. DAVID bioinformatics resources: expanded annotation database and novel algorithms to better extract biology from large gene lists. Nucleic Acids Res. 2007;35(Web Server issue):W169–75.

  109. 109.

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

  110. 110.

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

  111. 111.

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

Download references

Acknowledgements

We thank Embrapa and University of São Paulo for the collaborative efforts and CAPES for the scholarship to BSV.

Funding

This study was conducted with funding from Embrapa (Macroprograma 1, 01/2005) and FAPESP (process number 2012/23638–8). LCAR and LLC were granted CNPq fellowships.

Availability of data and materials

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

Author information

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

Correspondence to Bárbara Silva-Vignato.

Ethics declarations

Ethics approval and consent to participate

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

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Publisher’s Note

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

Additional files

Additional file 1: Table S1.

Test of means (t-test) of intramuscular fat (IMF), backfat thickness (BFT) and ribeye area (REA) between groups with High and Low REA and BFT in the Longissimus dorsi muscle of Nellore steers. (XLS 33 kb)

Additional file 2: Table S2.

Differentially expressed genes obtained between high and low genomic breeding value groups for ribeye area in Nellore steers with FDR 10%. The table contains the Ensembl gene identification, gene symbol, mean normalized counts, Log2 fold Change from low to high group, p value and p value adjusted for a false discovery rate of 10%. (XLS 45 kb)

Additional file 3: Table S3.

Differentially expressed genes obtained between high and low genomic breeding value groups for backfat thickness in Nellore steers with FDR 10%. The table contains the Ensembl gene identification, gene symbol, mean normalized counts, Log2 fold Change from low to high group, p value and p value adjusted for a false discovery rate of 10%. (XLS 31 kb)

Additional file 4: Table S4.

Functional enrichment analysis performed by DAVID (p value <0.1) comparing high and low genomic breeding value groups for backfat thickness. The table contains Gene Ontology category, identification and description, p value, p value adjusted for a false discovery rate of 10% and the DE genes for each category. (XLS 29 kb)

Additional file 5: Table S5.

Correspondence between sample names and European Nucleotide Archive (ENA) identifier (accession PRJEB13188 and PRJEB19421). (XLS 29 kb)

Additional file 6: Figures S1 and S2.

Plots of dispersion estimates of all ribeye area and backfat thickness genes. The black dots are the dispersion estimates of the empirical values of each gene; the red line represents the trend line; the blue dots represent the genes estimates regressed through the trend line used in the hypothesis test; and the blue circles above the “cloud” of points are genes which have high gene-wise dispersion estimates which are labelled as dispersion outliers, and will not be used in the hypothesis test. (PDF 110 kb)

Additional file 7: Figures S3 and S4.

Histogram of p-values from RNA-Seq data of Longissimus dorsi muscle of Nellore steers by DESeq2 software, for ribeye area and backfat thickness. Y-axis represents the frequency of p-values and x-axis represents the residual of p-values. (PDF 108 kb)

Rights and permissions

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

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Keywords

  • Backfat thickness
  • Bos taurus indicus
  • Ribeye area
  • RNA-Seq