Skip to main content

MiRNAs differentially expressed in skeletal muscle of animals with divergent estimated breeding values for beef tenderness



MicroRNAs (miRNAs) are small noncoding RNAs of approximately 22 nucleotides, highly conserved among species, which modulate gene expression by cleaving messenger RNA target or inhibiting translation. MiRNAs are involved in the regulation of many processes including cell proliferation, differentiation, neurogenesis, angiogenesis, and apoptosis. Beef tenderness is an organoleptic characteristic of great influence in the acceptance of meat by consumers. Previous studies have shown that collagen level, marbling, apoptosis and proteolysis are among the many factors that affect beef tenderness. Considering that miRNAs can modulate gene expression, this study was designed to identify differentially expressed miRNAs that could be modulating biological processes involved with beef tenderness.


Deep sequence analysis of miRNA libraries from longissimus thoracis muscle allowed the identification of 42 novel and 308 known miRNAs. Among the known miRNAs, seven were specifically expressed in skeletal muscle. Differential expression analysis between animals with high (H) and low (L) estimated breeding values for shear force (EBVSF) revealed bta-mir-182 and bta-mir-183 are up-regulated (q value < 0.05) in animals with L EBVSF, and bta-mir-338 is up-regulated in animals with H EBVSF. The number of bovine predicted targets for bta-mir-182, bta-mir-183 and bta-mir-338 were 811, 281 and 222, respectively, which correspond to 1204 unique target genes. Among these, four of them, MEF2C, MAP3K2, MTDH and TNRC6B were common targets of the three differentially expressed miRNAs. The functional analysis identified important pathways related to tenderness such as apoptosis and the calpain–calpastatin system.


The results obtained indicate the importance of miRNAs in the regulatory mechanisms that influence muscle proteolysis and meat tenderness and contribute to our better understanding of the role of miRNAs in biological processes associated with beef tenderness.


MicroRNAs (miRNAs) are small 22 nucleotides endogenous non-coding ribonucleic acids (RNAs) [1] that negatively modulate the expression of genes in plants, animals and virus at a post-transcriptional level through cleavage or translational inhibition [2, 3]. MiRNA sequences are highly conserved among species, from nematode to cattle and humans, a reason why they are of central importance to biology and developmental decisions [3,4,5,6,7]. Increasing evidence indicates that miRNAs play an important regulatory role in several biological processes such as cell proliferation, differentiation, neurogenesis, angiogenesis, and apoptosis as well as epigenetic changes [2, 8].

In animals, they were previously reported to be related to embryonic development and function of skeletal muscle, adipose, mammary and immune tissues [4, 8]. For example, miR-1 and miR-133 are muscle-specific and are involved in the modulation of muscle proliferation [9]; miR-133 increases proliferation of C2C12 myoblasts [9]; miR-486 is an inducer of myoblast differentiation [10], and miR-26a is induced during skeletal muscle regeneration [11]. However, little is known about the role of miRNAs in beef tenderness.

Among the traits of economic value in livestock species, meat quality, specifically beef tenderness, is considered the primary attribute of sensory satisfaction of the beef consumers [12,13,14]. It is a complex trait with economical importance to the beef industry and has been a major focus of many studies [15, 16].

The present investigation was undertaken to identify differentially expressed miRNAs and functional pathways associated with beef tenderness in Nelore (Bos indicus species) cattle. We hypothesized that variation in shear force at 14 days of aging could be associated with the difference in miRNA expression in skeletal muscle. Thus, we sequenced miRNAs from longissimus thoracis (LT) muscle of animals with high (H) and low (L) EBV for shear force (SF) values to detect differential expressed miRNAs and to identify putative biological processes associated with beef tenderness.


Animals and phenotype

The population used in this study was previously described in detail by Tizioto et al. [17]. Briefly, a total of 390 Nelore steers, offspring of 34 sires unrelated were used to obtain phenotypic data. The animals were raised at pasture until approximately 23 months of age when they were moved to a feedlot with identical nutrition and handling conditions. The animals were slaughtered at an average age of 25 months and an endpoint of 5 mm of backfat thickness (BFT). Immediately after exsanguination, samples were collected from the longissimus thoracis (LT) muscle between the 12th and 13th ribs and frozen in liquid nitrogen until RNA extraction. Measurements of meat tenderness were determined by the Warner–Bratzler shear force (WBSF) in 2.54 cm thick steaks obtained from the same muscle after aging at the 2 °C cold chamber for 24 h, at 7 and 14 days postmortem as described into detail by Carvalho et al. [18]. The WBSF values were calculated as the average of eight cores.

For this study, samples were ranked on estimated breeding values for shear force at 14 days of aging (EBVSF14) calculated from the previous study of our group [19], and we selected 34 animals with either the highest (H, n = 15) or lowest (L, n = 19) EBVSF14 to form the groups that were tested for miRNAs differential expression analysis.

RNA extraction and small RNAs libraries construction

The total RNA was extracted from 100 mg of frozen LT muscle using the TRIzol reagent (Life Technologies, Carlsbad, CA). RNA integrity (RIN) was verified by Bioanalyzer 2100 (Agilent, Santa Clara, CA, USA) and a minimum threshold of RIN seven was used for library construction. Small RNAs libraries were constructed from 1 μg of total RNA from each of the 34 samples using the Illumina TruSeq small RNA Sample Prep Kit (Illumina Inc., San Diego, CA, USA) according to the manufacturer’s protocol. PCR amplification was performed for 15 cycles. Library quality was determined using the High Sensitivity DNA Chip and an Agilent 2100 Bioanalyzer (Agilent Technologies) and quantified with qPCR with the KAPA Library Quantification kit (KAPA Biosystems, Foster City, CA, USA). The individual libraries were adjusted to 20 pM concentrations; sequencing was performed using a Miseq Reagent Kit v3 for 150 cycles in an Illumina Miseq Sequencing System (Illumina Inc., San Diego, CA, USA). This kit allows the generation of 25 million sequences reads per flow cell.

MiRNA sequencing data analysis

After sequencing, the Illumina CASAVA v1.8 pipeline was used to generate and de-multiplex the raw fastq sequences. The quality of Illumina deep sequencing data was determined by using the FastQC program (version 0.9.5) [20]. Adapters and low quality reads were trimmed using Cutadapt (version 1.2.1) [21] with the following parameters: − b AATCTCGTATGCCGTCTTCTGCTTGC-O 3-m 17-f fastq-q 24, where − b is the Illumina sequence adapter, -O indicates the minimum number of matching bases necessary to recognize the adapter, -m represents the minimum length of sequence and –q stands for the sequences quality.

Filtered reads were then processed following mirDeep2 analysis pipeline [22]. Sequences were aligned to the UMD3.1 Bos taurus taurus reference genome (available at the Ensembl database []) using the module. Only alignments with 0 mismatches in the seed region (first 18 nt of a read sequence) of a read mapped to the genome were retained.

Differentially expressed miRNAs

Initially, miRNAs with zero counts for all samples were removed. Next, the miRNAs presenting a count different from zero in at least 1/5 of the samples were maintained. Read count data was normalized to account for differences in starting RNA quantity and other library effects using the Upper Quartile method [23].

The QuasiSeq package [24] developed in R [25] was used to detect miRNAs differentially expressed. QuasiSeq uses a quasi-likelihood method through which is incorporated uncertainty in the estimated variances during the test for differential expression [24] and also provides a self-tuning approach to shrinking miRNAs dispersion estimates. The QuasiSeq method is built from a generalized linear model and allow fixed effects to be added in the model. The PROC MIXED procedure by SAS [26] was performed to determine significant (p < 0.05) fixed effects and covariates that were used in miRNA differential expression analysis by QuasiSeq package. In the QuasiSeq model, the final age and contemporary group (animals of the same slaughter group, origin and birth year) were fitted as a covariate and fixed effect, respectively. High and Low groups were defined for EBVSF at 14 days of aging and used to select animals, while the phenotype for QuasiSeq analysis was WBSF at 14 days of aging, as adopted by Gonçalves et al. [19]. Thus, it is important to highlight that differentially expressed miRNA were identified as a function of WBSF at 14 days of aging.

We used the q-value method to control the false discovery rate (FDR) at 5% [27], to correct for multiple tests (the total number of miRNA tested in the model).

Target gene prediction and functional analysis

The prediction of target genes of differentially expressed miRNAs was made using the TargetScan website ( and the B. taurus species. Only the expressed genes in muscle tissue were maintained, which were measured in a previous RNAseq study of our research group [28].

Gene Ontology analysis and functional annotation of target genes were done to identify significative canonical pathways (p-value < 0.1) from the target genes and also to visualize their networks and by IPA Core Analysis from Ingenuity® Pathway Analysis software (IPA®, QIAGEN Redwood City, CA).


Animals, phenotype and sequencing data

The animals with extreme values of EBVSF14 were selected to compose the high and low groups for miRNA-sequencing analysis (see “Methods”). The average values of SF at 14 days aging and EBVSF14 for the group with high EBVSF14 (H group; n = 15) and the group with low EBVSF14 (L group; n = 19) were 6.50 kgf/cm2 and 0.63; and 2.65 kgf/cm2 and − 0.62, respectively. A statistical test of means was performed between H and L groups, and the results indicate that the samples selected were divergent for SF at 24 h, seven and 14 days of aging, as well as EBVSF14, but were not different for intramuscular fat and ribeye area (Table 1).

Table 1 Test of means (t-test) between groups for shear force (kgf/cm2) at 24 h, seven and 14 days of aging, estimated breeding values at 14 days of aging (EBV), intramuscular fat (IMF, %) and ribeye area (REA, cm2)

A total of 40,513,215 reads were obtained from miRNA sequencing analysis. After removing low quality reads, adaptor’s sequences and insufficiently tagged, a total of 33,222,144 reads were ultimately maintained. The distribution of the reads length was from 20 to 25 nucleotides (nt), which is the typical length of the small RNA processed by Dicer.

The alignment of the reads against the Bos taurus UMD3.1 reference genome showed that ~ 87% of the reads were correctly mapped by the Bowtie2 algorithm using the miRDeep2 software [22]. The read mapping statistic is shown in Additional file 1.

Identification of known and novel miRNAs in LT muscle

Novel miRNAs in the bovine LT muscle were identified by a deep sequence of 34 small RNA libraries. Known miRNAs were identified by comparison between the precursor and mature miRNAs annotated in the miRBase database (Version 21). In total, 308 mature miRNAs were identified as known miRNAs, and among them, seven were previously described in the literature as skeletal muscle specific (Table 2). This was considered when the miRNA was expressed at least 20 times more than the average of its expression in the other tissues [29].

Table 2 Specific skeletal muscle miRNAs identified in samples of both phenotypic extremes

To identify novel miRNAs, we compared the seed sequence from miRNAs found in this study with miRNAs already annotated in miRBase for others species close phylogenetically to B. indicus. The seed sequence corresponds to the 2–8 nt region of the 5′ end of the miRNA, and is very important for recognition and silencing of mRNA [30, 31]. The B. taurus, H. sapiens, M. musculus, S. scrofa, and E. caballus species were chosen according to the database of phylogenetic trees of animal genes, TreeFam ( A total of 42 novel miRNAs was identified in the LT muscle from B. indicus (Additional file 2).

Identification of differentially expressed miRNAs

In order to identify miRNAs differentially expressed in LT muscle of animals with divergent estimated breeding values for beef tenderness the QuasiSeq R package was used. Three miRNAs were differentially expressed between the groups (q value < 0.05). Two miRNAs (bta-mir-182 and bta-mir-183) were up-regulated in the L group compared to the H group, and one miRNA (bta-mir-338) was up-regulated in H group compared to the L group (Table 3).

Table 3 Differentially expressed miRNAs between animals of the L and H groups based in estimated breeding values for shear force at 14 days of aging

Target genes prediction and functional pathways

To identify the potential functions of the differentially expressed miRNAs, bovine target genes of these miRNAs were predicted by TargetScan 6.2 ( From the list of target genes (Additional file 3), we selected only the genes that are expressed in LT muscle to perform the functional analysis (Fig. 1). For this selection, RNA sequencing (RNA-Seq) data of the same set of samples [28] of this study was used (Additional file 4).

Fig. 1

Summary of target genes predicted for bta-mir-182, bta-mir-183 and bta-mir-338 in all tissues and in Longissimus thoracis muscle specifically

The number of predicted targets expressed in LT muscle for bta-mir-182, bta-mir-183 and bta-mir-338 was 811, 281 and 222, respectively. These numbers correspond to 1204 target genes, after removing duplicates. Four genes (MEF2C, MAP3K2, MTDH, and TNRC6B) were found as common targets of the three differentially expressed miRNAs (Fig. 2).

Fig. 2

Venn chart of the common and specific target genes of differentially expressed miRNAs

The Ingenuity® Pathway Analysis software conducted the functional enrichment study. First, all predicted target genes were mapped to the IPA Knowledge Base, allowing us to reveal which molecules are being encoded by these genes (Table 4).

Table 4 Summary of molecules encoded by target genes predicted of differentially expressed miRNAs

Core analysis (functional analysis) revealed that gene expression, cellular development, cell morphology, cellular assembly and organization and cellular function and maintenance were ranked in the top of significant molecular and cellular functions (p < 0.05) (Fig. 3).

Fig. 3

Molecular and cellular functions of total predicted target genes of differentially expressed miRNAs. The likelihood association among the genes and biological functions is represented as − log(p-value), larger bars are more significant than shorter bars considering p-value < 0.05 as cutoff for significance (representing in the threshold vertical line)

The most significant biological processes associated with the target genes for the three differentially expressed miRNAs were cancer, post-translational modification, cell morphology, small molecule biochemistry, nucleic acid metabolism and vitamin and mineral metabolism (Table 5).

Table 5 Biological processes related to predicted targets genes

We identified 161, 112 and 61 significative (p < 0.1) canonical pathways from the list of target genes of the bta-mir-182, bta-mir-183, and bta-mir-338, respectively (Additional file 5). From these lists, we identified five canonical pathways associated with meat tenderness regulation: apoptosis signaling, glutathione biosynthesis, regulation of cellular mechanics by calpain protease and calcium signaling and transport. Table 6 shows the miRNA and targets genes involved in these pathways.

Table 6 Characterization of few candidate canonical pathways for meat tenderness


Increasing evidence indicate that miRNAs play an important regulatory role in several biological processes such as cell proliferation, differentiation, neurogenesis, angiogenesis, and apoptosis as well as epigenetic changes [2, 8], which could promote phenotypic variation among individuals. Among the phenotypic traits of interest for many researchers, beef tenderness is in evidence due to economical importance to the beef industry. Thus, the present investigation was undertaken to characterize the miRNAs expressed in skeletal muscle of Nelore cattle and to identify differentially expressed miRNAs and functional pathways associated with beef tenderness in Nelore (Bos indicus species) cattle.

Among the miRNAs preferentially expressed in LT muscle, miR-208b and miR-499 are located within introns of myosin genes [32]; the pairs miR-1-1/133a-2, miR-1-2/133a-1, and miR-206/133b are encoded by bicistronic transcripts on different chromosomes and have been shown to play roles in the control of muscle growth and differentiation [32]. The serum response factor (SRF) and myocyte enhancer factor-2 (MEF2) control the expression of miR-1-1/133a-2 and miR1-1–2/133a-1 [32]. These miRNAs exert opposing effects in the process mentioned above, with miR-1 playing a pro-apoptotic role and miR-133 playing an anti-apoptotic role in cardiomyocyte apoptosis [32, 33]. MiR-1 and miR-133 are evolutionary conserved and are found in most animal species, from Drosophila to human [9].

MiR-486 has been described as an inducer of myoblast differentiation through its negative regulation of PAX7, a transcription factor required for the biogenesis of muscle satellite cells and the specification of myogenic lineage [10]. MiR-206 also promotes muscle differentiation, as a previous study reported that the inhibition of this miRNA by antisense oligonucleotide inhibits cell cycle withdrawal and differentiation [34].

The differentially expressed miRNAs (bta-mir-182, bta-mir-183, and bta-mir-338) are predicted to modulate the expression of many genes that are involved in several pathways. Among the pathways identified, we chose to discuss the ones more likely to be associated with muscle proteolysis and meat tenderness.

It is widely accepted that the mechanism of meat tenderizing is an enzymatic process involving multiple proteolytic systems and that apoptosis is the first step in the conversion of muscle into meat [35]. There is large variability in meat tenderness, and this variability comes mainly from the biochemical reactions taking place during this conversion, immediately after slaughter [36].

Apoptosis is a physiological process of cell death. It is very important for development and tissue homeostasis and is mediated by a particular group of cysteine peptidases called caspases [35]. These proteins are divided according to their location in the apoptosis pathway in apoptosis initiator caspases (caspase 8, 9, 10 and 12) and effector caspases (caspases 3, 6 and 7) [37]. These enzymes play a vital role in the induction, amplification, and transduction of intracellular apoptotic signals [38] and it has been suggested that there is an overexpression of these proteases in tender meat [39].

The apoptosis pathway signaling is presented in Fig. 4, where the relationship between the protein caspase and the anti-apoptotic protein BCL2 can be observed. Previous evidence report that an increase in the intracellular concentration of calcium in apoptosis results in the activation of calpain, and with this, blocking of the BCL2 protein family [40]. It has also been shown that disturbance in calcium homeostasis in the endoplasmic reticulum can lead to activation of caspase 12 by calpain and this last one can cleave the anti-apoptotic protein BCL-XL (a member of Bcl-2 family) turning it into a pro-apoptotic protein [40].

Fig. 4

Targets molecules of bta-mir-182 involved in the apoptosis signaling pathway. Solid lines indicate a direct connection while broken lines indicate an indirect relationship

We, therefore, hypothesized that up-regulation of bta-mir-182 in the low EBVSF14 animals could be down-regulating the anti-apoptotic BCL2 protein level, and thus promoting apoptosis that could contribute to muscle proteolysis and tenderness. Conversely, it is important to note that bta-mir-182 could be downregulating CAPN5, CASP2, and CASP9, and this could reduce the apoptotic process.

In the apoptosis signaling pathway, it is also possible to observe the presence of ROCK-1, one of the two isoforms of Rho-kinase (ROCK). This protein was reported as a direct substrate of caspase-3 protein in driving the apoptosis of myocytes [41], and it also promotes actin–myosin-mediated contractile force generation by phosphorylating its target proteins. The myosin binding subunit of myosin light chain (MLC) phosphatase 1 (MYPT1) myosin light chain 2 (MLC2) and LIM kinases are downstream substrates of ROCK, modulating the organization of actin cytoskeleton [42].

Another important pathway identified in the study was the glutathione pathway. Glutathione (GSH) is a cysteine-containing tripeptide (glutamine, glycine, and cysteine) involved in the antioxidation system and intracellular redox state. It can be present in two forms, the reduced (GSH) and the oxidized glutathione (GSSG), and the ratio of the two forms allow the characterization of the oxidative stress in cells [43].

D’Alessandro and colleagues [44], in a study of Chianina cattle, described an accumulation of oxidative stress in tender meat samples, when GSSG/GSH ratios were higher than in tough ones. Oxidative stress might be related to meat tenderness due to the promotion of reactive oxygen species (ROS) inducing protein fragmentation [44]. The majority of GSH is found in the cytosol, but a small percentage is located in the mitochondria, contributing to the protection of this organelle from reactive oxygen species (ROS) [45]. Thus, the decrease of GSH can increase the reactive oxygen species (ROS) or accelerate the mitochondrial damage [44].

The up-regulation of bta-mir-183 in the L group could be downregulating the target gene Glutamate-cysteine ligase (GCLM), a rate-limiting enzyme of glutathione synthesis. The decrease in GSH can increase ROS, stimulate apoptosis and consequently contribute to tenderness.

In this context, it is important to note that in a previous study using the same animals and phenotypes used here, a QTL for SF at 24 h was found located on BTA23 at 24 Mb, regions containing the glutathione S-transferase alpha gene family [17].

Two important target genes found in our study are calpain and calpastatin. Calpain is a calcium-activated protease and injection of calcium in muscles accelerates postmortem proteolysis and tenderization [46]. In post-mortem muscle, calcium concentration in the cytoplasm increases gradually during rigor mortis while the sarcoplasmic reticulum is emptied [13]. This translocation of calcium results in different processes affecting the permeability of the sarcoplasmic reticulum membrane as binding of pro-apoptotic Bcl2 members [40]. Calpastatin (CAST) is a calpain proteolytic enzyme inhibitor. Calpain’s role in post mortem transformation of muscle in meat is extensively studied and is widely accepted that the proteolytic activity of this protein contributes to tenderness [47, 48]. Previous studies showed that the correlations between different tenderness rates in beef, pork, and lamb are inversely related to the calpain: calpastatin ratio [49].

Calpastatin is a predicted target for bta-mir-182, suggesting that high expression levels of this miRNA would be reducing the translation of the CAST gene, resulting in a lower inhibitory effect on calpain, higher post-mortem proteolytic activity and consequently greater tenderness (Fig. 4).

Proteins of the heat shock protein family (HSP27 and HSP70) were found in this study as targets of bta-mir-338 differentially expressed miRNA. Our results corroborate previous research from our group with the same animals of this work, where both proteins (HSP70 e HSP27) appeared as down-regulated in animals with lower values for shear force [18]. Others authors also reported low levels of HSP27 and HSP70-1A/B associated with animals with more tender meat [50, 51]. This evidenced negative relationship between certain HSP levels and meat tenderness could be linked to the anti-apoptotic activity of these proteins [47].

Regarding the common targets identified for the three differentially expressed miRNAs, we highlight the myocyte-specific enhancer factor 2C (MEF2C), mitogen-activated protein kinase kinase kinase 2 (MAP3K2), and metadherin (MTDH). The MEF2C transcription factor is restricted to skeletal muscle, brain, and spleen, playing a crucial role in the morphogenesis and myogenesis. A previous study identified genetic variants of this gene in seven different breeds of cattle, such as, Aberdeen Angus, Charolais, Hereford, Limousin, Simmental, Polish Friesian and Polish Red, which could constitute as potential genetic markers for the characteristics of carcass and meat quality in cattle [52]. Besides that, the MEF2C bovine gene has been mapped on chromosome 7, which contains quantitative trait loci (QTL) responsible for the average daily gain, body and carcass weight [53] as well as the fat thickness in the Longissimus muscle [54]. Both the biological importance and the chromosomal location suggests that the bovine MEF2C gene could be a promising functional and positional candidate gene responsible for carcass and meat quality traits in cattle [55].

The MAP3K2 gene encodes a member of the serine/threonine protein kinase family that preferentially mediates the activation of other kinases involved in MAP kinase signaling pathway [56]. MAP3K2 is involved in cell proliferation, differentiation, and cell migration [57]. Previous studies have reported that MAP3K2 can promote cell proliferation in different types of cancer [58,59,60] and that some miRNAs can suppress the tumor by targeting MAP3K2 [60,61,62]. On the other hands, the overexpression of MAP3K2 inhibited cell proliferation in chickens but did not induce apoptosis [63]. The association of the MAP3K2 markers with loin muscle area (LMA) was previously identified in a Duroc pig population, evaluating QTL for carcass merit and meat quality traits [64], suggesting that this gene could be considered suitable candidates for future studies of growth traits and meat production in domestic animals including cattle.

The MTDH gene was also identified in a previous study, as a target gene for microRNA bta-mir-885, which was exclusively expressed in the semitendinosus muscle (STD) from Japanese black cattle when compared to masseter muscle (MS). The functional annotation of MTDH gene revealed its possible relationship in skeletal system development and regulation of transcription, respectively [65].


The results obtained indicate the importance of miRNAs in the regulatory mechanisms that influence muscle proteolysis and meat tenderness and contribute to our better understanding of the role of miRNAs in biological processes associated with beef tenderness. Further studies are necessary to explore the implementation of these miRNAs as biomarkers in cattle breeding contributing to the selection of animals with improved meat tenderness.



estimated breeding value


estimated breeding value for shear force at 14 days


false discovery rate




ingenuity pathway analysis




longissimus thoracis




quantitative trait loci


  1. 1.

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

    CAS  PubMed  Article  Google Scholar 

  2. 2.

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

    CAS  Article  Google Scholar 

  3. 3.

    McDaneld TG. MicroRNA: mechanism of gene regulation and application to livestock. J Anim Sci. 2009;87(14 Suppl):E21–8.

    CAS  PubMed  Article  Google Scholar 

  4. 4.

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

    CAS  PubMed  Article  Google Scholar 

  5. 5.

    Kedde M, Agami R. Interplay between microRNAs and RNA-binding proteins determines developmental processes. Cell Cycle. 2008;7(7):899–903.

    CAS  PubMed  Article  Google Scholar 

  6. 6.

    Ibanez-Ventoso C, Vora M, Driscoll M. Sequence relationships among C. elegans, D. melanogaster and human microRNAs highlight the extensive conservation of microRNAs in biology. PLoS ONE. 2008;3(7):e2818.

    PubMed  PubMed Central  Article  Google Scholar 

  7. 7.

    Zhao Y, Srivastava DA. A developmental view of microRNA function. Trends Biochem Sci. 2007;32(4):189–97.

    CAS  PubMed  Article  Google Scholar 

  8. 8.

    Cupp AS, Matthews J, Huff-Lonergan E, Spurlock DM, McLean D. Cell biology symposium: the role of microRNA in cell function. J Anim Sci. 2009;87(14 Suppl):E19–20.

    CAS  PubMed  Article  Google Scholar 

  9. 9.

    Chen JF, Mandel EM, Thomson JM, Wu Q, Callis TE, Hammond SM, Conlon FL, Wang DZ. The role of microRNA-1 and microRNA-133 in skeletal muscle proliferation and differentiation. Nat Genet. 2006;38(2):228–33.

    CAS  PubMed  Article  Google Scholar 

  10. 10.

    Dey BK, Gagan J, Dutta A. miR-206 and -486 induce myoblast differentiation by downregulating Pax7. Mol Cell Biol. 2011;31:203–14.

    CAS  PubMed  Article  Google Scholar 

  11. 11.

    Dey BK, Gagan J, Yan Z, Dutta A. miR-26a is required for skeletal muscle differentiation and regeneration in mice. Genes Dev. 2012;26(19):2180–91.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  12. 12.

    Ouali A. Meat tenderization: possible causes and mechanisms: a review. J Muscle Foods. 1990;1(2):129–65.

    Article  Google Scholar 

  13. 13.

    Koohmaraie M. The role of Ca(2+)-dependent proteases (calpains) in post mortem proteolysis and meat tenderness. Biochimie. 1992;74(3):239–45.

    CAS  PubMed  Article  Google Scholar 

  14. 14.

    Miller MF, Carr MA, Ramsey CB, Crockett KL, Hoover LC. Consumer thresholds for establishing the value of beef tenderness. J Anim Sci. 2001;79(12):3062–8.

    CAS  PubMed  Article  Google Scholar 

  15. 15.

    Thompson J. Managing meat tenderness. Meat Sci. 2002;62(3):295–308.

    PubMed  Article  Google Scholar 

  16. 16.

    Hocquette JF, Wezemael LV, Chriki S, Legrand I, Verbeke W, Farmer L, Scollan ND, Polkinghorne R, Rodbotten R, Allen P, Pethick DW. Modelling of beef sensory quality for a better prediction of palatability. Meat Sci. 2014;97(3):316–22.

    PubMed  Article  Google Scholar 

  17. 17.

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

    CAS  PubMed  Article  Google Scholar 

  18. 18.

    Carvalho ME, Gasparin G, Poleti MD, Rosa AF, Balieiro JC, Labate CA, Nassu RT, Tullio RR, Regitano LC, Mourão GB, et al. Heat shock and structural proteins associated with meat tenderness in Nellore beef cattle, a Bos indicus breed. Meat Sci. 2014;96(3):1318–24.

    CAS  PubMed  Article  Google Scholar 

  19. 19.

    Gonçalves TM, Regitano LCA, Koltes JE, Cesar ASM, Andrade SCS, Mourão GB, Gasparin G, Moreira GCM, Fritz-Waters E, Reecy JM, Coutinho LL. Gene Co-expression analysis indicates potential pathways and regulators of beef tenderness in Nellore cattle. Front Genet. 2018;9:441.

    PubMed  PubMed Central  Article  Google Scholar 

  20. 20.

    Andrews S. FastQC: a quality control tool for high throughput sequence data. Cambridge: Babraham Institute; 2010.

    Google Scholar 

  21. 21.

    Martin M. Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnet J. 2011;17:10–2.

    Article  Google Scholar 

  22. 22.

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

    PubMed  Article  Google Scholar 

  23. 23.

    Bullard JH, Purdom E, Hansen KD, Dudoit S. Evaluation of statistical methods for normalization and differential expression in mRNA-Seq experiments. BMC Bioinformatics. 2010;11(1):94.

    PubMed  PubMed Central  Article  Google Scholar 

  24. 24.

    Lund SP, Nettleton D, McCarthy DJ, Smyth GK. Detecting differential expression in RNA-sequence data using quasi-likelihood with shrunken dispersion estimates. Stat Appl Genet Mol Biol. 2012;11:8.

    Article  Google Scholar 

  25. 25.

    Team RDC. R: a language and environment for statistical computing. Vienna: The R Foundation for Statistical Computing; 2011.

    Google Scholar 

  26. 26.

    Saxton AM. Genetic analysis of complex traits using SAS. Cary, NC: SAS Institute; 2004. p. 312.

    Google Scholar 

  27. 27.

    Storey JD, Tibshirani R. Statistical significance for genomewide studies. Proc Natl Acad Sci USA. 2003;100:9440–5.

    CAS  PubMed  Article  Google Scholar 

  28. 28.

    Cesar AS, Regitano LC, Koltes JE, Fritz-Waters ER, Lanna DP, Gasparin G, Mourão GB, Oliveira PS, Reecy JM, Coutinho LL. Putative regulatory factors associated with intramuscular fat content. PLoS ONE. 2015;10(6):e0128350.

    PubMed  PubMed Central  Article  Google Scholar 

  29. 29.

    Lee EJ, Baek M, Gusev Y, Brackett DJ, Nuovo GJ, Schmittgen TD. Systematic evaluation of microRNA processing patterns in tissues, cell lines, and tumors. RNA. 2008;14:35–42.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  30. 30.

    Lewis BP, Burge CB, Bartel DP. Conserved seed pairing, often flanked by adenosines, indicates that thousands of human genes are microRNA targets. Cell. 2005;120:15–20.

    CAS  PubMed  Article  Google Scholar 

  31. 31.

    Lewis BP, Shih IH, Jones-Rhoades MW, Bartel DP, Burge CB. Prediction of mammalian microRNA targets. Cell. 2003;115(7):787–98.

    CAS  PubMed  Article  Google Scholar 

  32. 32.

    van Rooij E, Liu N, Olson EN. MicroRNAs flex their muscles. Trends Genet. 2008;24(4):159–66.

    PubMed  Article  Google Scholar 

  33. 33.

    Liu N, Williams AH, Kim Y, McAnally J, Bezprozvannaya S, Sutherland LB, Richardson JA, Bassel-Duby R, Olson EN. An intragenic MEF2-dependent enhancer directs muscle-specific expression of microRNAs 1 and 133. Proc Natl Acad Sci USA. 2007;104(52):20844–9.

    CAS  PubMed  Article  Google Scholar 

  34. 34.

    Kim HK, Lee YS, Sivaprasad U, Malhotra A, Dutta A. Muscle-specific microRNA miR-206 promotes muscle differentiation. J Cell Biol. 2006;174:677–87.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  35. 35.

    Ouali A, Herrera-Mendez CH, Coulis G, Becila S, Boudjellal A, Aubry L, Sentandreu MA. Revisiting the conversion of muscle into meat and the underlying mechanisms. Meat Sci. 2006;74(1):44–58.

    PubMed  Article  Google Scholar 

  36. 36.

    Boudida Y, Gagaoua M, Becila S, Picard B, Boudjellal A, Herrera-Mendez CH, Sentandreu M, Ouali A. Serine protease inhibitors as good predictors of meat tenderness: which are they and what are their functions? Crit Rev Food Sci Nutr. 2016;56(6):957–72.

    CAS  PubMed  Article  Google Scholar 

  37. 37.

    Earnshaw WC, Martins LM, Kaufmann SH. Mammalian caspases: structure, activation, substrates, and functions during apoptosis. Annu Rev Biochem. 1999;68:383–424.

    CAS  PubMed  Article  Google Scholar 

  38. 38.

    Fan TJ, Han LH, Cong RS, Liang J. Caspase family proteases and apoptosis. Acta Biochim Biophys Sin (Shanghai). 2005;37(11):719–27.

    CAS  Article  Google Scholar 

  39. 39.

    Ouali A, Gagaoua M, Boudida Y, Becila S, Boudjellal A, Herrera-Mendez CH, Sentandreu MA. Biomarkers of meat tenderness: present knowledge and perspectives in regards to our current understanding of the mechanisms involved. Meat Sci. 2013;95(4):854–70.

    CAS  PubMed  Article  Google Scholar 

  40. 40.

    Gil-Parrado S, Fernandez-Montalvan A, Assfalg-Machleidt I, Popp O, Bestvater F, Holloschi A, Knoch TA, Auerswald EA, Welsh K, Reed JC, et al. Ionomycin-activated calpain triggers apoptosis. A probable role for Bcl-2 family members. J Biol Chem. 2002;277(30):27217–26.

    CAS  PubMed  Article  Google Scholar 

  41. 41.

    Chang J, Xie M, Shah VR, Schneider MD, Entman ML, Wei L, Schwartz RJ. Activation of Rho-associated coiled-coil protein kinase 1 (ROCK-1) by caspase-3 cleavage plays an essential role in cardiac myocyte apoptosis. Proc Natl Acad Sci USA. 2006;103(39):14495–500.

    CAS  PubMed  Article  Google Scholar 

  42. 42.

    Shi J, Wu X, Surma M, Vemula S, Zhang L, Yang Y, Kapur R, Wei L. Distinct roles for ROCK1 and ROCK2 in the regulation of cell detachment. Cell Death Dis. 2013;4:e483.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  43. 43.

    Huang J, Jia Y, Li Q, Son K, Hamilton C, Burris WR, Bridges PJ, Stromberg AJ, Matthews JC. Gluathione content and expression of proteins involved with glutathione metabolism differs in longissimus dorsi, subcutaneous adipose, and liver tissues of finished vs growing beer steers. J Anim Sci. 2018;96:5152–65.

    PubMed  Google Scholar 

  44. 44.

    D’Alessandro A, Marrocco C, Rinalducci S, Mirasole C, Failla S, Zolla L. Chianina beef tenderness investigated through integrated Omics. J Proteomics. 2012;75(14):4381–98.

    PubMed  Article  Google Scholar 

  45. 45.

    Fernandez-Checa JC, Kaplowitz N, Garcia-Ruiz C, Colell A, Miranda M, Mari M, Ardite E, Morales A. GSH transport in mitochondria: defense against TNF-induced oxidative stress and alcohol-induced defect. Am J Physiol. 1997;273:G7–17.

    CAS  PubMed  Article  Google Scholar 

  46. 46.

    Lian T, Wang L, Liu Y. A new insight into the role of calpains in post-mortem meat tenderization in domestic animals: a review. Asian-Australas J Anim Sci. 2013;26:443–54.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  47. 47.

    Beere HM. Death versus survival: functional interaction between the apoptotic and stress-inducible heat shock protein pathways. J Clin Invest. 2005;115(10):2633–9.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  48. 48.

    Pearson G, Robinson F, Beers Gibson T, Xu BE, Karandikar M, Berman K, Cobb MH. Mitogen-activated protein (MAP) kinase pathways: regulation and physiological functions. Endocr Rev. 2001;22(2):153–83.

    CAS  PubMed  Google Scholar 

  49. 49.

    Koohmaraie M, Whipple G, Kretchmar DH, Crouse JD, Mersmann HJ. Postmortem proteolysis in longissimus muscle from beef, lamb and pork carcasses. J Anim Sci. 1991;69(2):617–24.

    CAS  PubMed  Article  Google Scholar 

  50. 50.

    Bouley J, Chambon C, Picard B. Proteome analysis applied to the study of muscle development and sensorial qualities of bovine meat. Sci. Aliment. 2003;23:75–8.

    CAS  Article  Google Scholar 

  51. 51.

    Picard B, Gagaoua M, Micol D, Cassar-Malek I, Hocquette JF, Terlouw CEM. Inverse relationships between biomarkers and beef tenderness according to contractile and metabolic properties of the muscle. J Agric Food Chem. 2014;62(40):9808–18.

    CAS  PubMed  Article  Google Scholar 

  52. 52.

    Flisikowski EJK, Msscp PÁ, Promoter ÁRÁ. Nucleotide sequence and variations of the bovine myocyte enhancer factor 2C (MEF2C) gene promoter in Bos Taurus cattle. Mol Biol Rep. 2011;38:1269–76.

    PubMed  Article  Google Scholar 

  53. 53.

    Casas E, Shackelford M, Smith TPL, Stone RT. Detection of quantitative trait loci for growth and carcass composition in cattle. J Anim Sci. 2003;81:2976–83.

    CAS  PubMed  Article  Google Scholar 

  54. 54.

    Ferraz JB, Pinto LF, Meirelles FV, Eler JP, De Rezende FM, Oliveira EC, Almeida HB, Woodward B, Nkrumah D. Association of single nucleotide polymorphisms with carcass traits in Nellore cattle. Genet Mol Res. 2009;8:1360–6.

    CAS  PubMed  Article  Google Scholar 

  55. 55.

    Juszczuk-Kubiak E, Wicinska K, Starzynski RR. Postnatal expression patterns and polymorphism analysis of the bovine myocyte enhancer factor 2C (Mef2C) gene. Meat Sci. 2014;94:753–8.

    Article  Google Scholar 

  56. 56.

    Greenblatt MB, Shin DY, Oh H, Lee K-Y, Zhai B, Gygi SP, Lotinum S, Baron R, Liu D, Su B, Glimcher LH, Shim J-H. MEKK2 mediates an alternative β-catenin pathway that promotes bone formation. PNAS. 2016;113(9):E1226–35.

    CAS  PubMed  Article  Google Scholar 

  57. 57.

    Yu J, Tan Q, Deng B, Fang C, Qi D, Wang R. The microRNA-520a-3p inhibits proliferation, apoptosis and metastasis by targeting MAP3K2 in non-small cell lung cancer. Am J Cancer Res. 2015;5(2):802–11.

    CAS  PubMed  PubMed Central  Google Scholar 

  58. 58.

    Cazares LH, Troyer D, Mendrinos S, Lance RA, Nyalwidhe JO, Beydoun HA, Clements MA, Drake RR, Semmes OJ. Imaging mass spectrometry of a specific fragment of mitogen-activated protein kinase/extracellular signal-regulated kinase kinase kinase 2 discriminates cancer from uninvolved prostate tissue. Clin Cancer Res. 2009;15(17):5541–51.

    CAS  PubMed  Article  Google Scholar 

  59. 59.

    Mirza AA, Kahle MP, Ameka M, Campbell EM, Cuevas BD. MEKK2 regulates focal adhesion stability and motility in invasive breast cancer cells. Biochim Biophys Acta. 2014;1843(5):945–54.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  60. 60.

    Zhang W, Kong G, Zhang J, Wang T, Ye L, Zhang X. MicroRNA-520b inhibits growth of hepatoma cells by targeting MEKK2 and cyclin D1. PLoS ONE. 2012;7:e314502012.

    Google Scholar 

  61. 61.

    Huang T, She K, Peng G, Wang W, Huang J, Li J, Wang Z, He J. MicroRNA-186 suppresses cell proliferation and metastasis through targeting MAP3K2 in non-small cell lung cancer. Int J Oncol. 2016;49:1437–44.

    CAS  PubMed  Article  Google Scholar 

  62. 62.

    Wang LL, Zhang M. miR-582-5p is a potential prognostic marker in human non-small cell lung cancer and functions as a tumor suppressor by targeting MAP3K2. Eur Rev Med Pharmacol Sci. 2018;22:7760–7.

    PubMed  Google Scholar 

  63. 63.

    Zhang X, Song H, Qiao S, Liu J, Xing T, Yan X, Li H, Wang N. MiR-17-5p and miR-20a promote chicken cell proliferation at least in part by upregulation of c-Myc via MAP3K2 targeting. Sci Rep. 2017;7(1):15852.

    PubMed  PubMed Central  Article  Google Scholar 

  64. 64.

    Choi I, Bates RO, Raney NE, Steibel JP, Ernst CE. Evaluation of QTL for carcass merit and meat quality traits in a US commercial Duroc population. Meat Sci. 2012;92:132–8.

    PubMed  Article  Google Scholar 

  65. 65.

    Muroya S, Taniguchi M, Shibata M, Oe M, Ojima K, Nakajima I, Chikuni K. Profiling of differentially expressed microRNA and the bioinformatic target gene analyses in bovine fast- and slow-type muscles by massively parallel sequencing. J Anim Sci. 2013;91:90–103.

    CAS  PubMed  Article  Google Scholar 

Download references

Authors’ contributions

BIGK, LCAR, and LLC conceived the idea of this research. BIGK, GG and MDP participated in the experimental design and miRNA-sequencing. BIGK and ASMC performed data analysis. BIGK and MDP drafted the manuscript. BIGK, MDP, ASMC, LCAR, GG, GCMM and LLC collaborated with interpretation and discussion of the results. LCAR and LLC provided the experimental environment, phenotype and data analysis support. All authors read and approved the final manuscript.


We thank the Brazilian Agricultural Research Corporation (Embrapa), the National Council for Scientific and Technological Development (CNPq) and the São Paulo State Research Foundation (FAPESP) for financial support. Berna I. G. Kappeler, Mirele D. Poleti and Aline S. M. Cesar are FAPESP fellows (Grants 2014/00943-5, 2013/21017-9 and 2014/11871-5). This study was financed in part by the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior—Brasil (CAPES)—Finance Code 001, and Conselho Nacional de Desenvolvimento Científico e Tecnológico—Brasil (CNPq).

Competing interests

The authors declare that they have no competing interests.

Consent for publication

Not applicable.

Data availability

The dataset supporting the conclusions of this article is available in the European Nucleotide Archive (ENA) repository (EMBL-EBI), under accession PRJEB18514 [].

Ethics approval and consent to participate

All experimental procedures involving steers were approved by the Institutional Animal Care and Use Committee Guidelines from Brazilian Agricultural Research Corporation—EMBRAPA (process number: Macroprograma 1, 01/2005) and sanctioned by the president Dr. Rui Machado.


The design of this study, and the phenotype collection, miRNA-sequencing, analysis and interpretation of data were conducted with funding from EMBRAPA (Macroprograma 1, 01/2005), FAPESP (processes number 2012/23638-8), and CAPES (Finance Code 001). The writing of the manuscript and fellowships were supported by FAPESP (2014/00943-5, 2014/11871-5, 2013/21017-9) and CNPq (150914/2017-2).

Publisher’s Note

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

Author information



Corresponding author

Correspondence to Luiz L. Coutinho.

Additional files

Additional file 1.

The number of raw-reads, number of reads after cleaning (filtered), number and percentage of mapped reads for High and Low groups based on estimated breeding values for shear force.

Additional file 2.

Novel miRNAs.

Additional file 3.

Bovine target genes of bta-mir-182 (sheet 1), bta-mir-183 (sheet 2), and bta-mir-338 (sheet 3).

Additional file 4.

List of Longissimus thoracis muscle expressed genes in Nelore.

Additional file 5.

Significative (p < 0.1) canonical pathways from the list of target genes of bta-mir-182 (Sheet 1), bta-mir-183 (Sheet 2) and bta-mir-338 (Sheet 3).

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (, 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 ( applies to the data made available in this article, unless otherwise stated.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Kappeler, B.I.G., Regitano, L.C.A., Poleti, M.D. et al. MiRNAs differentially expressed in skeletal muscle of animals with divergent estimated breeding values for beef tenderness. BMC Molecular Biol 20, 1 (2019).

Download citation


  • Beef
  • Bos indicus
  • bta-miR
  • MicroRNA
  • Shear force