Gene Co-expression Analysis Indicates Potential Pathways and Regulators of Beef Tenderness in Nellore Cattle

Beef tenderness, a complex trait affected by many factors, is economically important to beef quality, industry, and consumer’s palatability. In this study, RNA-Seq was used in network analysis to better understand the biological processes that lead to differences in beef tenderness. Skeletal muscle transcriptional profiles from 24 Nellore steers, selected by extreme estimated breeding values (EBVs) for shear force after 14 days of aging, were analyzed and 22 differentially expressed transcripts were identified. Among these were genes encoding ribosomal proteins, glutathione transporter ATP-binding cassette, sub-family C (CFTR/MRP), member 4 (ABCC4), and synaptotagmin IV (SYT4). Complementary co-expression analyses using Partial Correlation with Information Theory (PCIT), Phenotypic Impact Factor (PIF) and the Regulatory Impact Factor (RIF) methods identified candidate regulators and related pathways. The PCIT analysis identified ubiquitin specific peptidase 2 (USP2), growth factor receptor-bound protein 10 (GBR10), anoctamin 1 (ANO1), and transmembrane BAX inhibitor motif containing 4 (TMBIM4) as the most differentially hubbed (DH) transcripts. The transcripts that had a significant correlation with USP2, GBR10, ANO1, and TMBIM4 enriched for proteasome KEGG pathway. RIF analysis identified microRNAs as candidate regulators of variation in tenderness, including bta-mir-133a-2 and bta-mir-22. Both microRNAs have target genes present in the calcium signaling pathway and apoptosis. PIF analysis identified myoglobin (MB), enolase 3 (ENO3), and carbonic anhydrase 3 (CA3) as potentially having fundamental roles in tenderness. Pathways identified in our study impacted in beef tenderness included: calcium signaling, apoptosis, and proteolysis. These findings underscore some of the complex molecular mechanisms that control beef tenderness in Nellore cattle.

Introduction

Meat tenderness is an important beef palatability trait for consumers, and most of the biological mechanisms involved in this trait are not completely understood (Ouali et al., 2006). Beef tenderness is commonly measured by SF and is influenced by genetic factors, environment and genotype-environment interaction (Koch et al., 1982; Wheeler et al., 2005). The average heritability for beef tenderness is about 0.3, corroborating that genetics have a reasonable contribution for this trait (Koch et al., 1982). Unfortunately, tenderness and other phenotypes associated with meat quality such as taste and color, are difficult to measure directly, and can only be determined after an animal is slaughtered (Gonçalves, 2015; Mateescu et al., 2017).

For complex traits such as beef tenderness, next-generation sequencing (NGS) technologies have been applied in transcriptomic analyses to identify DE genes and biological pathways potentially involved with meat quality phenotypes (Metzker, 2010). However, biological interpretation of RNA-Seq remains a challenge (Barabasi and Oltvai, 2004; de la Fuente, 2010). Co-expression analysis provides a powerful approach to integrate the transcripts relationship into different pathways based on phenotypic variation (Reverter et al., 2006; Reverter and Chan, 2008). The PCIT, RIF, and PIF algorithms determine the changes in correlation or wiring of a network between treatments or extreme levels of phenotypes. PCIT identifies significant partial correlations among pairs of genes and can be used to quantify differential connectivity as well as the strength of a connection to a hub gene (Reverter et al., 2006; Reverter and Chan, 2008). PIF and RIF analyses are based on the approach of differential wiring or differential correlation, which is the difference in co-expression between a certain pair of genes under two different situations. More specifically, PIF is the average expression of two conditions combined, multiplied by the DE value of genes. PIF scores allow the variability and differential expression of a gene to be prioritized in determining its importance in relation to differences in phenotypic variation (Reverter et al., 2010). RIF is calculated from the cumulative differential wiring of each regulator compared to the target DE genes and weighted for PIF (Hudson et al., 2009). The RIF algorithm can identify genes that may not show differential expression but still alter the molecular “wiring” of a pathway or act as biomarkers of phenotypic variation. For instance, myostatin was identified as a regulator in a study of gene expression was conducted between two breeds of cattle, Piedmontese; which is double-muscled due to a mutation in the myostatin gene, and Wagyu, which do not contain the mutation. The key regulator of muscle growth, myostatin, was not DE in animals carrying the mutation, but this gene was differentially wired based on the RIF analysis (Hudson et al., 2009).

Previous gene co-expression networks analysis allowed the integration of regulatory processes, which affected gene expression level of the studied system (Reverter et al., 2006). Complex networks of gene interactions influence meat quality, including beef tenderness (Mateescu et al., 2017). Understanding the processes involved in postmortem traits impacting meat quality is critical to elucidate the mechanisms involved in the conversion of muscle to meat (Ponsuksili et al., 2013). The discovery of potential regulators by gene co-expression analysis includes small RNAs. MiRNAs are small non-coding RNAs that are highly conserved among mammals and can regulate translation by binding to the 3′-UTR of target mRNAs. They play a role in many biological processes including myogenesis and hypertrophy (Wang et al., 2014). In muscle, miRNAs can regulate myoblast proliferation and differentiation and may control gene expression related to the variation of beef tenderness during stress (Zhao et al., 2012a). The identification of miRNAs enriched in skeletal muscle can enhance the understanding of its biology and function (Wang et al., 2014). Thus, gene co-expression analysis can be used to possibly identify miRNAs that regulate the skeletal muscle development.

The objective of this study was to identify DE genes associated with beef tenderness in skeletal muscle (LD) from Nellore cattle selected for extreme estimated breeding values (EBVs) for SF after 14 days of aging. A secondary goal was to unveil possible mechanisms associated with biological processes and variation in beef tenderness at the pathway level. The final goal was to explore regulators involved in beef tenderness variation using co-expression methods (PCIT, PIF, and RIF) in order to increase our understanding of genes and pathways involved in the variation of beef tenderness.

Materials and Methods

Animals and Phenotypic Data

The Institutional Animal Care and Use Committee (IACUC) Guidelines from the Brazilian Agricultural Research Corporation – EMBRAPA approved all experimental procedures involving steers prior to conducting this study (Macroprograma 1, 01/2005). The animals selected for RNA-Seq analysis came from a larger project where meat quality traits were evaluated (Tizioto et al., 2013a). A total of 310 Nellore steers sired by 33 unrelated Nellore bulls were studied. The sires were selected based on the national summary of Nellore cattle to represent the main breeding lineages used in Brazil. Steers were raised in feedlots with identical nutritional rations and handling conditions until slaughter at an average age of 25 months as previously reported (Cesar et al., 2014). Samples from skeletal muscle (LD) located between the 12th and 13th ribs were collected at two-time points: at slaughter for RNA sequencing analysis, and 24 h after slaughter for SF analysis and aging.

The beef tenderness was represented by the SF measured from the steaks with 2.54 cm of thickness. The SF was obtained by using a TA XT2i (Stable Micro Systems Ltd., Surrey, United Kingdom) texture analyzer coupled to a Warner-Bratzler blade, at 1.016 mm thickness (Wheeler et al., 1997). Measurements of SF were made after aging in a 2°C chamber for 24 h at 7 and 14 days after slaughter (Tizioto et al., 2013a). The measurements at 14 days of aging confirmed the results of tenderness identified at 7 days of aging and were used as the SF phenotype.

Statistical Analysis of Shear Force

A goal of the larger overarching meat quality project in Nellore was to determine the impact of additive genetic control on meat quality traits and to develop genetic breeding tools to improve these traits in Nellore cattle. For this reason, animals were ranked based on the EBV, instead of the phenotypic value. The EBV represents the animal’s predicted genetic merit, half of which is inherited by its progeny to impact their performance. The benefits of using EBV are the ability (1) to make full use of information from all relatives; (2) to account for fixed environmental effects so that animals can be compared across groups, giving wider scope for selection by focusing on genetic effects (Sosnicki and Newman, 2010). Therefore, EBVs were used to select the extreme groups of animals for RNA-Seq, PCIT, PIF, and RIF analyses.

Prior to RNA-Seq analysis, EBVSF were obtained based on the SF performance of the half-sib data using standard best linear unbiased prediction (BLUP) procedures under an animal model performed using the mixed procedure of SAS (Henderson, 1973; Saxton and SAS Institute, 2004; Mrode, 2014). Animals were ranked based on the EBVSF at 14 days of aging in order to adjust for fixed effects that could interfere with phenotypic measurements. Selected animals were divided into two groups with extreme values: highest (H) values of EBVSF (n = 11) and lowest (L) values of EBVSF (n = 13). The mixed model used to estimate breeding values contained a numerator coefficient matrix of relationship (A) containing 284 animals. The basic model equation in matrix notation was as follows:

$$\mathrm{\gamma }=X\mathrm{\beta }+\mathit{Zu}+\mathrm{\epsilon }$$

Where the γ (dependent variable) is the phenotype SF phenotype in kgf/cm², and X is the incidence matrix of fixed effects associated with the vector of fixed effectsβ. Fixed effects included contemporary group (based on slaughter group and birth season) and age, which was used as a linear covariate. The incidence matrix of random animal effects in u, where u ∼ N (0, Aσ_a²) for SF at 14 days of aging was represented by Z. The random animal effect was included to estimate the additive genetic effect. The ε is the vector of random residuals, normally and independently distributed as ε∼ (NID) (0, Iσ_e²). Variance component estimation was conducted by restricted maximum likelihood (REML) methodology (Patterson and Thompson, 1971).

RNA-Seq Analysis

Total RNA of 24 LD muscle samples (100 mg) were extracted with Trizol reagent (Thermo Fisher Scientific, Carlsbad, CA, United States). The RNA was quantified using a spectrophotometer (NanoDrop 200 - Thermo Fisher Scientific, Wilmington, DE, United States). RNA Integrity Number (RIN) was verified using a 1% agarose gel and a Bioanalyzer 2100 (Agilent Technologies, Santa Clara, CA, United States). All samples had RIN values greater than seven.

Libraries were prepared using the TruSeq SBS Kit v3-HS (200 cycles). The average cDNA library size was estimated using an Agilent Bioanalyzer 2100 (Agilent Technologies, Santa Clara, CA, United States) and quantified using quantitative PCR with the KAPA Library Quantification kit (KAPA Biosystems, Foster City, CA, United States). Samples were diluted to 17 pM and clustered on the sequencing flowcell by cBOT (Illumina, San Diego, CA, United States) using the TruSeq paired-end (PE) Cluster Kit v3-cBOT-HS kit. Libraries were barcoded, multiplexed, and sequenced in seven lanes on a HiScanSQ (Illumina, San Diego, CA, United States) to generate 100 nt-long paired-end reads. These procedures were performed at the Center for Functional Genomics (USP, Campus “Luiz de Queiroz”, SP, Brazil).

Sequence quality was evaluated with FastQC (Andrews, 2010). Phred quality score of the reads was assessed using the software Seqyclean (Zhbannikov et al., 2013). Reads with nucleotide quality score < 26 and reads length < 65 bp were discarded. Potential vector and adaptor sequences were compared to the UniVec database to remove possible artificial sequence contaminants. Sequence reads were aligned to the UMD3.1 Bos taurus masked reference genome (Ensembl release 80) using TopHat 2.0.10 (Trapnell et al., 2009; Zimin et al., 2009). Reads with more than one mismatch or non-unique mapping were discarded. A reference-guided assembly of transcripts was performed using Cufflinks v2.1.1 to identify novel transcripts (Trapnell et al., 2010). Read counts for each transcript were quantified with HTSeq v0.5.4p2, and used for differential expression analysis (Anders and Huber, 2010).

Differential Expression Analysis

Prior to the DE analysis, raw count data was filtered for zero (non-expressed) and low expression to ensure transcripts were expressed across the samples (Lund et al., 2012). First, all transcripts with zero counts were removed from all samples as these genes were assumed to be non-expressed in skeletal muscle. Second, all transcripts were required to have at least two reads to remove very lowly expressed genes. Finally, transcripts were required to have non-zero read counts in 1/5 of all samples to remove transcripts that may be quantified in error. In this study, transcripts with no reads in six or more samples were dropped from the analysis. Read count data was normalized for differences in library size using the 75% Quartile method to account for differences in RNA quantity and other library effects (Bullard et al., 2010).

The R package QuasiSeq, based on quasi-likelihood methods added to the edgeR software, was used to estimate the read count dispersion for each gene (Lund et al., 2012). The use of QuasiSeq for DE analysis allowed fixed effects to be added into the original generalized linear model as well as a variable to account for over-dispersion to reduce the error in modeling differences in gene expression. An ANOVA model was conducted in SAS to determine what variables might be important to the differential expression analysis. Final age nested within year, contemporary group (slaughter group, origin, and birth season) and lane of sequencing within flow cell were included as fixed effects in the QuasiSeq gene expression model since they were significant in the SAS analysis (P-value < 0.05). The QLSpline effect fitted a quasi-negative binomial distribution to account for over-dispersion in the count data (Lund et al., 2012). This approach identified DE transcripts as a function of SF at 14 days of aging. High and low groups were defined for EBVSF at 14 days of aging and used to select the animals, while the phenotype for QuasiSeq analysis was SF at 14 days of aging. The groups H and L for EBVSF at 14 days of aging were used to discuss the association of transcripts with the SF phenotype. The direction of expression will be up or down-regulated based on the b values, and the average of the normalized counts are enough to represent the direction of expression. The q-value methodology was used to control the false discovery rate (FDR) at 10% (Storey and Tibshirani, 2003).

Annotation, Functional Analysis, and miRNA Prediction

The Ensembl 80 annotation file used herein contains entries representing both protein coding genes and types of RNAs such as rRNA, snoRNA, snRNA, miRNA, lincRNA, and pseudogenes. When annotation was unavailable for genes identified as DE, several approaches were used to obtain annotations in this study. Novel transcripts identified by Cufflinks were annotated using the Genome2Seq and GOanna tools at AgBase (McCarthy et al., 2010; Trapnell et al., 2010). Uncharacterized transcripts were annotated by searching for orthologous genes at BioMart using Ensembl v80 database (Cunningham et al., 2014). CUFFs transcripts (short transcripts assembled by Cufflinks or transfrags in GTF format) were annotated by the best sequence match generated by the NCBI Basic Local Alignment Search Tool (BLAST) (Altschul et al., 1997).

Functional enrichment analysis is a very effective way to provide biological knowledge to a large dataset. The functional analysis of curated gene ontology terms was conducted with the Database for Annotation, Visualization and Integrated Discovery (DAVID) v6.7 tool to identify Gene Ontology (GO) terms and KEGG pathways with significant enrichment scores (Huang et al., 2008). The Benjamini and Hochberg correction was used to account for multiple testing (FDR < 10%) in this process (Benjamini and Hochberg, 1995).

Standard RNA-Seq protocols that capture poly-A RNAs will, as an added benefit, capture the expression levels of miRNA precursors (pre-miRNAs) as they are initially polyadenylated. It was not the main objective of this study to identify miRNAs, but since it was possible to retain longer pre-miRNAs, we were able to subsequently identify predicted target genes and functions associated with the targets. Pre-miRNAs identified in this study exhibited transcript lengths higher than 65 bp and were reported by the name of the annotated miRNA according to Ensembl database and miRBase. A miRNA target prediction method was used to allow the exploration of miRNA targets potentially regulated by a miRNA (Betel et al., 2008). Transcript targets were predicted with Miranda v3.3a software for miRNAs of interest using the Target Scanning Algorithm (Enright et al., 2003; Betel et al., 2008). Predicted transcript targets of each miRNA were analyzed for over-representation of pathway annotations using Ingenuity Pathway Analysis (IPA) (Enright et al., 2003).

PCIT, PIF, and RIF Co-expression Analysis

A modified version of PCIT was used to calculate the direct and partial correlations of all possible pair-wise transcript interactions (Reverter and Chan, 2008; Koesterke et al., 2013). The significance of these pair-wise interactions was determined using an information theory approach (Reverter and Chan, 2008). In order to account for gene length to compare transcripts to one another, gene expression values were transformed to fragment per kilobase per million mapped reads (FPKM) values. The FPKM values were corrected for the same covariates used for the QuasiSeq model to create residual FPKMs for each transcript used in the PCIT analysis. The differential hubbing or differential connectivity is the difference in the significant connections (i.e., correlations) of a gene between the two experimental conditions under study. This is computed by subtracting the significant connections a gene has in the two conditions (Hudson et al., 2009). Only significant partial correlation values ≥ 0.90 were used in the DH analysis to determine differential connectivity of hub genes. Thus, the DH analysis considered only the hub genes. DH was computed as the difference of significant connections of a transcript between the H versus the L group (Reverter and Chan, 2008; Hudson et al., 2009; Koesterke et al., 2013). For each DH transcript, a list of all positively and all negatively correlated transcripts was generated. The transcripts with a significant correlation with DH transcripts are called DH transcript targets, and they were tested by gene set enrichment analysis to determine if DH targets were biologically meaningful based on the pathways impacted. These data were analyzed using DAVID to identify over-represented GO terms and KEGG pathways, where the FDR was controlled using Benjamini and Hochberg methodology (FDR < 10%) (Benjamini and Hochberg, 1995). The PCIT results were visualized by BioLayout software (Enright and Ouzounis, 2001).

In addition, a modified RIF analysis was used to identify changes in the correlation of a network between the H and L SF groups (Yang et al., 2015). Instead of limiting possible regulators to transcription factors, all transcripts (including those that mapped on miRNA entries present in Ensembl v80 annotation file) were tested as possible regulators in this analysis. To determine the importance of each DE transcript with respect to phenotype, PIF was also calculated for each transcript. RIF1 and RIF2 scores were calculated as described by Reverter et al. (2010) and modified to consider all transcripts as regulators as described by Yang et al. (2015) using residual FPKM gene expression values, correlation values from PCIT (≥0.90) and PIF values for all transcripts (Reverter et al., 2010; Yang et al., 2015).

Results

Differential expression analysis identified 22 DE transcripts (q-value < 0.10) for SF at 14 days of aging. The top DH transcripts identified by PCIT analysis and the functional analyses performed on the correlated transcripts were associated with the proteasome KEGG pathway. Two miRNAs bta-mir-133a-2 and bta-mir-22 were identified among the top regulatory candidates by RIF analysis. The myoglobin (MB), enolase 3 (ENO3), and carbonic anhydrase 3 (CA3) were top transcripts identified by PIF analysis.

Phenotypic and Sequencing Results

Previously, as part of the larger phenotypic study of this population, the average SF at 24 h, seven and 14 days were estimated at 8.70, 5.93, and 4.56 kgf/cm², respectively (Tizioto et al., 2013a). In order to focus on the genetic mechanisms that influence variation in SF at 14 days of aging, animals with high and low EBVSF were selected for RNA-Seq analysis. The extreme animals selected for this study had average SF values of 6.52 kgf/cm² for the H group and 2.73 kgf/cm² for the L group at 14 days of aging. The selected animals were among the 5% most extreme for SF at 14 days of aging relative to the full population. Table 1 provides the phenotypes for each group of individuals selected for RNA-Seq, the average measurements of SF at different levels of aging, the EBV values, and the EBV accuracies. The phenotypes were consistent with a normal aging process, with a reduction in SF from 24 h to 14 days for both the H and L groups (Table 1).

TABLE 1

TABLE 1. Average shear force (kgf/cm²), EBVSF values and corresponding accuracies calculated for the low (L) group and high (H) groups at seven and 14 days of aging.

On average, 31 million sequence reads were generated per sample from sequencing, and Seqyclean removed 25% of these reads. Approximately 19 million reads were aligned to the 3.1 UMD Bos taurus genome per sample. On average, 84% of the reads were mapped, and of these, 77% aligned uniquely (reads that aligned once in the genome). Read mapping statistics, including the total number of reads mapped and the percentage of reads mapped uniquely to the genome, are provided in Supplementary Table A.

Identification of Differentially Expressed Transcripts

A total of 15,693 annotated transcripts were identified after data filtering, including 1,743 novel transcripts. From the 22 transcripts identified as DE (q-value < 0.10) by QuasiSeq, eight had higher expression in the H group, and 14 transcripts had higher expression in the L group. All the DE transcripts were annotated, except for one where no annotation could be identified (Table 2). A complete list of the relative expression levels is presented in the Supplementary Table B.

TABLE 2

TABLE 2. Differentially expressed transcripts obtained between H and L groups for EBVSF at 14 days of aging in Nellore steers with q-value < 0.10.

The low number of DE transcripts may be due to the complex statistical model used to perform the analysis, including many variables. For this reason, there was no significantly enriched GO term and KEGG pathway from the DE transcripts based on DAVID analysis (FDR adjusted P-value < 0.1).

Hub Transcripts Identified With PCIT

PCIT analysis identified transcripts that may be acting as regulators of pathways impacting tenderness. The DH analysis identified correlated transcripts with an absolute correlation ≥ 0.9. Table 3 reports the top 10 negative and positive DH transcripts when comparing the H group relative to L group.

TABLE 3

TABLE 3. Top 10 negatively and positively differentially hubbed (DH) transcripts.

Functional enrichment analysis was performed on the targets of DH transcripts to determine their biological relevance. Each of the top DH transcripts, independent of their DH value, exhibited positively and negatively correlated transcript targets that were enriched for GO terms by DAVID tools (FDR adjusted P-value < 0.1). The Supplementary Tables C–V present the most significant (FDR adjusted P-value < 0.1) enrichment results, including GO terms of biological process (BP), molecular function (MF), cellular component (CC), and KEGG pathway information for all negatively and all positively correlated transcript targets of the DH transcripts presented in Table 3.

Targets of top DH transcripts enriched for terms associated with proteasome GO terms or KEGG pathways. Four transcripts, including: tripartite motif containing 45 (TRIM45), ubiquitin specific peptidase 2 (USP2), growth factor receptor-bound protein 10 (GRB10), and HECT domain-containing e3 ubiquitin protein ligase 4 (HECTD4) exhibited negative DH between the High (H) and Low (L) groups for EBVSF at 14 days of aging and were enriched for proteasome (Figure 1). The transcripts transmembrane BAX inhibitor motif containing 4 (TMBIM4), transmembrane protein 150a (TMEM150A), eukaryotic translation initiation factor 2b, subunit 1 alpha, 26kda (EIF2B1), and CDK5 regulatory subunit associated protein 3 (CDK5RAP3) exhibited positive DH between the high (H) and Low (L) groups for EBVSF at 14 days of aging, and were enriched for terms related to proteasome function (Figure 2). Figures 1, 2 were constructed using the BioLayout software to visualize DH between the H and L groups for EBVSF at 14 days of aging. The figures selected represented the top DH transcripts related to proteasome function and had the most different networks.

FIGURE 1

FIGURE 1. Negative DH between the H and L groups for EBVSF at 14 days of aging. The center spot represents a hub transcript with a high value of DH, red lines represent positive DH and blue lines represent negative DH. All lines represent significant partial correlations ≥ 0.90 at p-value < 0.05 based on the information theory approach, which accounts for all comparisons in the analysis.

FIGURE 2

FIGURE 2. Positive DH between the H and L groups for EBVSF at 14 days of aging. The center spot represents a hub transcript with a high value of DH, the red lines represent positive DH and blue lines represent negative DH. All lines represent significant partial correlations ≥ 0.90 at p-value < 0.05 based on the information theory approach, which accounts for all comparisons in the analysis.

Potential Regulatory Transcripts Identified by RIF and PIF

Co-expression analyses were performed to identify transcripts associated with differences in the correlation wiring of a network (RIF1) and the magnitude of changes in expression of transcripts associated with tenderness (RIF2). The transcripts identified by RIF1 and RIF2 analyses are shown in Tables 4, 5. A total of 13,485 transcripts identified by RIF1 had negative values, and 52 transcripts had positive values. For RIF2, a total of 50 transcripts identified by RIF2 had negative values, and 13,487 transcripts had positive values.

TABLE 4

TABLE 4. Top negative and positive transcripts identified by Regulatory Impact Factor 1 (RIF1) Z-scores.

Functional enrichment analysis was performed by DAVID using the list of the top 3,000 negative and positive for both RIF1 and RIF2. The significant (FDR adjusted P-value < 0.05) enriched GO terms are presented in the Supplementary Table W. There were no significant enriched GO terms for transcripts with positive RIF1values and negative RIF2 values. Enriched categories for negative RIF1 include the GO BP terms carboxylic acid metabolic process, cellular ketone metabolic process, and lipid metabolic process, while the enriched categories for positive RIF2 include the GO BP terms cellular metabolic process, intracellular organelle, and membrane-bounded organelle (Supplementary Table W).

PIF analysis identified 995 significant transcripts (FDR adjusted P-value < 0.05) that could impact variation in tenderness. Among the top 10 transcripts identified by PIF, potential biomarkers of tenderness variation included MB, ENO3 and CA3 (Table 6). Functional analysis of significant PIF transcripts identified GO terms involved in the cellular metabolic process, nitrogen compound metabolic process, biosynthetic process and primary metabolic process (FDR adjusted P-value < 0.1), and are presented in Supplementary Table X.

Some potential regulators identified by both RIF and PIF analyses showed association with mitochondria activity, such as the positive RIF1 transcript mitochondrially encoded cytochrome C oxidase I (MT-CO1) (Table 4), the positive RIF2 transcripts ATPase family, AAA domain containing 3A (ATAD3A) and histidine triad nucleotide binding protein 2 (HINT2) (Table 5), and two transcripts from cytochromes c oxidase complex identified by PIF: cytochrome C oxidase subunit VIIa polypedtide 1 (muscle) (COX71A) and cytochrome C oxidase subunit VIIb (COX7B) (Table 6).

TABLE 5

TABLE 5. Top negative and positive transcripts identified by Regulatory Impact Factor 2 (RIF2) Z-scores^∗.

TABLE 6

TABLE 6. Top 10 transcripts identified by Phenotypic Impact factor (PIF) analysis with FDR adjusted P-value < 0.05.

Transcripts and Pathways Regulated by miRNAs

The miRNAs identified in this study had sequences higher than 65bp, and for this reason they were considered miRNA precursors. Bta-mir-22 (ENSBTAT00000042309; miRBase; Acc: MI0011372) had 95 bp of transcript length (including UTRs and CDS), and bta-mir-133a-2 (ENSBTAT00000042285; miRBase; Acc: MI0009732) contained 87 bp of transcript length (including UTRs and CDS) (Tables 4, 5). Since both bta-mir-133a-2 and bta-mir-22 were identified as potential regulators by RIF analysis, the target transcripts of this miRNA were identified and analyzed for enriched pathways. The target transcripts of miRNAs bta-mir-133a-2 and bta-mir-22 were enriched for 251 and 308 significant IPA pathways, respectively. The complete list of significant IPA pathways is presented in the Supplementary Tables Y, Z. Both miRNAs identified by the co-expression analysis share the same targets with lethal (let-7) miRNA, and the bta-mir-22 is also a target of bta-mir-133a-2. The two miRNAs targeted a variety of structural protein transcripts. The potential regulator bta-mir-133a-2, identified by both positive RIF1 and negative RIF2 as a candidate regulator, had 23 myosin transcript targets, including 10 from the myosin heavy chain class, two from the myosin light chain class, and 11 unconventional myosins. The canonical pathway actin cytoskeleton signaling was enriched for both bta-mir-133a-2 and bta-mir-22 (P-value = 4.90E-06 and 7.24E-08, respectively). The bta-mir-133a-2 and bta-mir-22 had some overlapping targets, such as the eukaryotic translation elongation factor 1 alpha 2 (EEF1A2), which is a paralog of the eukaryotic translation elongation factor 1 alpha 1 (EEF1A), identified as a DE transcript. Another DE transcript targeted by bta-mir-133a-2 was the apolipoprotein L3 (APOL3). DH transcript ubiquitin specific peptidase 2 (USP2) was a target of bta-mir-22, while tripartite motif containing 45 (TRIM45) was a transcript target of bta-mir-133a-2. Syntaxin-binding protein 6 (amisyn) (STXBP6), alpha-kinase 3 (ALPK3), and anoctamin 1, calcium activated chloride channel (ANO1) were also transcript targets of both miRNAs.

Some of the predicted transcript targets identified by differential expression and co-expression analysis, such as the DE transcript synaptotagmin IV (SYT4) and many DH transcripts, were related to calcium signaling. Calpastatin (CAST) and two genes from the Calpain family were targets of bta-mir-22, including: calpain 1, (mu/I) large subunit (CAPN1) and calpain 11 (CAPN11). Two other genes from the Calpain family were transcript targets of both miRNAs: calpain 3, (p94) (CAPN3) and the calpain 7 (CAPN7). The canonical calcium-signaling pathway was enriched for transcript targets of both bta-mir-133a-2 (P-value = 8.71E-07) and bta-mir-22 (P-value = 4.79E-04). The bta-mir-133a-2 targets also enriched for calcium-induced T lymphocyte apoptosis (P-value = 7.76E-04) and calcium transport I (P-value = 5.62E-03) pathways.

Multiple proteins related to mitochondrial function were identified by co-expression analyses. The bta-mir-133a-2 and bta-mir-22 targets enriched for the mitochondrial L-carnitine shuttle pathway (P-value = 3.16E-02 and 1.12E-02, respectively). In addition, while glutathione reductase (GSR) was a target of bta-mir-22, glutathione synthetase (GSS) was a target of bta-mir-133a-2. Bta-mir-133a-2 had glutathione biosynthesis (P-value = 6.03E-02) as a significant enriched pathway. Both miRNAs had glutathione S-transferase mu 4 (GSTM4) and glutathione S-transferase, alpha 4 (GSTA4) as common targets.

Discussion

The Relationship of Differentially Expressed Transcripts With Beef Tenderness

This study integrates DE genes and co-expression analysis associated with SF in Nellore cattle. It also contrasts differences in EBVs to focus on changes in gene expression associated with additive genetic differences in SF phenotypic value. Genetic variation in complex traits can be divided into many components like additive, dominance, and interaction effects of genes (Hill et al., 2008). Complex models fit to identify additive variation are capable of detecting effects that are not exclusively additive, such as epistatic effects, when estimating EBVs. Increasing knowledge about biological pathways and gene networks implies that gene-gene interactions (i.e., epistasis) might be important in trait variation (Hill et al., 2008). Nowadays, there are methods for integrating the “omics” and genetic information, creating systems genetics approaches as potential tools for use in animal breeding to study complex traits (Kadarmideen et al., 2006). The aim of these methodologies is to complement, add, and integrate quantitative genetics and transcriptomics to improve model accuracy, as the RNA-Seq provides an unprecedented level of accuracy and precision for measurements of gene expression profiles (Kadarmideen et al., 2006). In a study of boar taint in non-castrated pigs, the authors used summarized EBV to select animals for transcriptomics analysis (Drag et al., 2017). The co-expression analysis revealed modules that were associated with the summarized EBVs, indicating that EBVs were reliable indicators of phenotype variation (Drag et al., 2017). In addition, other studies successfully linked differences in additive genetic merit (EBV) with RNA-Seq based DE genes, pathways, and network analyses (Kogelman et al., 2014). Thus, these reports corroborate that this approach is reliable for differential expression studies.

The identification of 22 DE genes by QuasiSeq is due to the more conservative approach of this methodology, which adjusts for over-dispersion and provides the additional advantage of sharing information across genes when estimating dispersions (Lund et al., 2012). In a study comparing methodologies for RNA-Seq analysis, QuasiSeq was the top ranked method (Burden et al., 2014). Beef tenderness is a complex trait that can be influenced by many factors (i.e., breed, muscle, nutrition) (Bongiorni et al., 2016). As such, it is not unexpected that the number of DE genes identified as a function of beef tenderness would be prone to large variability based on the influence of these many contributing factors. An RNA-Seq study of beef tenderness was performed between two different Italian breeds, Chianina and Maremmana, and within each breed using LD muscle (Bongiorni et al., 2016). These two breeds were selected for different reasons: Maremmana cattle represent robustness and Chianina cattle represent conformation and fast growth. Regardless of breed, the group tender versus tough showed 494 DE genes (Bongiorni et al., 2016). Within Chianina breed, 607 genes were DE, while 163 genes were identified within the Maremmana breed, showing that the number of DE genes for tenderness can vary in different breeds. Most of the DE genes that differed by breed were involved in the regulation of cellular and biological processes (Bongiorni et al., 2016). Another transcriptomic study of LD muscle in Angus steers using microarray and RT-PCR analysis identified 53 DE genes for beef tenderness (Zhao et al., 2012b). Many of the genes identified in this study were involved in skeletal muscle contraction and lipid metabolism (Zhao et al., 2012b). Finally, an RNA-Seq study of tenderness in Nellore cattle identified 40 DE genes, however, just 17 had known functions (Fonseca et al., 2017). This number of DE genes is consistent with our results, and some of the genes with known functions were also related to ubiquitin metabolism and transport of calcium (Fonseca et al., 2017). A GWAS analysis was previously performed in the same population of Nellore cattle, and many QTLs were identified for meat quality traits. However, just a few QTLs were associated with tenderness, and most with small effects (Tizioto et al., 2013b). A proteomic study was also performed in the same population evaluating Longissimus thoracis muscle to identify DE proteins related to meat tenderness in Nellore cattle and only 17 DE spots were identified (Carvalho et al., 2014). These proteins were involved in muscle structure, cell organization, and metabolism (Carvalho et al., 2014). The functional roles of many of the transcripts identified in our study are consistent with the biological roles of genes and proteins previously associated with tenderness. DE transcripts identified included the following: components of the ribosome, serine proteases, and regulators of calcium and lipid transport (Gonçalves, 2015).

Postmortem processes result in changes in structural protein genes, such as shortening of the muscle due to rigor mortis, which in turn influences beef tenderness (Lonergan et al., 2010). In this study, protein synthesis may be a factor impacting muscle SF because of the presence of five ribosomal proteins identified as DE transcripts: 60S ribosomal protein L9 (RPL9), 60S ribosomal protein L31 (RPL31), 60S ribosomal protein L36 (RPL36), 40S ribosomal protein S14 (RPS14), and 40S ribosomal protein S15 (RPS15) (Table 2). These five ribosomal component transcripts plus EEF1A1 are expressed at higher levels in the L group, which represents the group with more tender meat (lower values of EBVSF at 14 days of aging) (Supplementary Table B). Ribosomal proteins interact with elongation factor EEF1 during protein biosynthesis (Graifer and Karpova, 2015). A tenderness study in male and female Qinchuan cattle identified high expression of genes related to translational elongation and protein biosynthesis, such as EEF1A1, ribosomal protein S27a (RPS27A), and eukaryotic translation initiation factor 1 (EIF1) in tender meat (Zhang et al., 2011). Protein synthesis was not measured in our study; however, based on the weight gain observed during the feedlot stage, there is no indication of difference in protein accretion between the L and H groups. Thus, a better hypothesis might be that animals in the L group have a higher rate of protein turnover, which is an essential process for cell survival (Voet et al., 2013).

Serine, 2 (trypsin 2) (PRSS2) protein is a serine protease that is involved in many biological processes like collagen metabolic process, proteolysis, serine-type endopeptidase activity, apoptosis, and calcium ion binding (Mihalyi, 1953; Sentandreu et al., 2002; De Bruin et al., 2003; Vandenabeele et al., 2005; Abo-Ismail et al., 2013). Serine proteases are proteolytic enzymes responsible for digestion of polypeptides (Di Cera, 2009). Proteolytic enzymes like calpains are important for the meat aging process. Since calpains are known to be important in meat tenderness, it is plausible that changes in the expression or function of PRSS2 may impact SF in Nellore cattle (Calkins and Sullivan, 2007). However, although some authors found that serine proteases could make meat more tender and may participate in muscle proteolysis, in this study, PRSS2 had higher expression in the H group, which represents the group with tougher meat (higher values of EBVSF at 14 days of aging) (Uytterhaegen et al., 1994; Alarcon-Rojo and Dransfield, 1995; Boudida et al., 2016). A previous study performed in the same population identified a quantitative trait loci (QTL) associated with SF at 24 h after slaughter, and the associated SNP within this QTL is located near the serine proteinase inhibitor 2 (SERPIN2) on BTA2 (Tizioto et al., 2013b). All LD muscle samples used in this study were collected immediately after slaughter and frozen immediately for posterior RNA extraction, which suggests that at the time of the sample collection, PRSS2 transcript levels were highly abundant. One possible hypothesis is that PRSS2 might have been degraded by the presence of serpins before the aging process starts, making the meat tough. Thus, understanding the high expression level of PRSS2 in the H group might be important to determine the if this mechanism impacts beef tenderness differences between Bos taurus and Bos indicus cattle.

The family of synaptotagmins are transmembrane proteins that act as calcium regulators of exocytosis (Südhof, 2002; Johnson et al., 2010). SYT4 is a paralog of synaptotagmin 11 (SYT11), which is a known as a calcium sensor and plays a role in the regulation of the synaptic transmission (Glavan et al., 2009). Calcium is required to recover from injury, and synaptic exocytosis is a similar process as that used by muscle cells repair the membrane damage (Steinhardt et al., 1994; Lek et al., 2013). Calpains are calcium-binding proteins that can also participate in cell membrane repair (Mellgren et al., 2009; Lek et al., 2013). Calpains and synaptotagmins appear to interact and possibly function in tandem in some cell types (Fukuda and Itoh, 2004; Aganna et al., 2006). As the calpains, synaptotagmins may also be involved in meat tenderness. A SNP in the synaptotagmin 9 (SYT9) was described as the most significant SNP associated with slice SF in swine (Nonneman et al., 2013). In our study, SYT4 had higher expression in the H group, which may indicate that calcium concentration is not significantly reduced in the muscle before the aging process starts. Therefore, the consequent toughness of the meat might be related to understanding why the calcium is not being used by proteolytic enzymes, which is a common process for meat tenderization during the 14 days of aging.

Apolipoproteins are related to lipid transporter activity, and share characteristics with members of the apoptotic Bcl-2 family (Uzureau et al., 2016). Some members of this family like the apolipoprotein A5 (APOA5) and apolipoprotein C3 (APOC3) are regulators of triglyceride, lipoprotein metabolism, cholesterol, LDL, and HDL (O’Brien et al., 2005). Apolipoproteins were associated with meat tenderness in swine (Hui et al., 2013). SNPs in both APOA5 and APOC3 genes were associated with various meat quality traits like cooked weight percentage, drip loss rate, meat color value of LD muscle, and SF (Hui et al., 2013). In this study, the APOL3 transcript had higher expression in the H groups that may suggest that lipid transportation and metabolism can impact beef tenderness process.

The ATP-binding cassette, sub-family C (CFTR/MRP), member 4 (ABCC4), is also known as multidrug resistance-associated protein, which are ATP-dependent-export and anion transporters (Hammond et al., 2007). A substrate of ABCC4 is glutathione (GSH), which exhibits greater excretion in the presence of ABCC4 (Bai et al., 2004). The export of GSH by ABCC4 is observed in the cell during apoptosis. A decrease of GSH can lead to an increase in the presence of reactive oxygen species (ROS) or accelerate the mitochondrial damage. When GSH is not exported, cell survival increases (Hammond et al., 2007). In this study, ABCC4 had higher expression in the L group, indicating a potential increase in the exportation of GSH from the cell and higher levels of oxidation. In agreement with this hypothesis, another study in Chianina cattle indicated that the presence of ROS can induce protein fragmentation and consequently impact beef tenderness (D’Alessandro et al., 2012). Oxidative stress was also impacts glycolysis and heat shock proteins responses, which are mechanisms known to positively effect beef tenderness (D’Alessandro et al., 2012). In addition, a previous study performed in the same population identified a QTL that contained the glutathione S-transferase alpha gene family located on BTA23 at 24 Mb that accounted for 0.11% of the additive genetic variance for SF at 24 h of aging (Tizioto et al., 2013b).

Our study indicates that despite the low number of DE genes detected, the DE transcripts identified might trigger molecular mechanisms elucidating the differences in SF between the L and H group in Nellore cattle. The identification of molecular targets can effectively be applied in studies of animal breeding and genetics to understand the differences in the meat quality of Nellore cattle and other breeds.

Differential Hubbing Transcripts Affecting Beef Tenderness

The DH analysis identified transcripts related to proteasome and ubiquitin function, as well as calcium signaling, consistent with previous associations with SF (Zhao et al., 2012b). Fourteen of the top hub genes identified by PCIT were all associated with protein degradation, indicating a large re-wiring (i.e., change in the correlation) of the protein degradation pathways between the H and L groups. One of these transcripts is USP2, which is from the family of de-ubiquitinating enzymes and can be expressed as two distinct isoforms with opposing effects in muscle differentiation. While overexpression of USP2a stimulates differentiation, the overexpression of USP2b inhibits this process (Wing, 2013). USP2 was previously associated with beef tenderness in a transcriptome study between male and female Qinchuan cattle (Zhang et al., 2011). In addition, another gene from the same family, ubiquitin specific peptidase 32 (USP32), was associated with tenderness in Nellore cattle (Fonseca et al., 2017). TRIM45, from the TRIM family ubiquitin ligases, is a component of the ubiquitin-proteasome and is involved in a variety of functions like growth regulation, apoptosis, and may suppress transcriptional activities of mitogen-activated protein kinase (MAPK) signaling pathway targets (Shibata et al., 2012). In agreement with these results, TRIM45 and another gene from the same family, tripartite motif 32 (TRIM32), were respectively associated with beef tenderness in Chianina and Maremmana cattle (Bongiorni et al., 2016).

HECT domain containing E3 ubiquitin protein ligase 4 (HECTD4) is an E3 ubiquitin-protein, which was identified as DH. Negative correlates of HECTD4 enriched for proteasome function using the KEGG pathway database. The ubiquitin-proteasome pathway (UPP) controls the protein turnover and degradation of intracellular proteins (Glickman and Ciechanover, 2002; Lecker et al., 2006). The E3 ubiquitin-protein ligase is a component of the UPP and the main enzyme responsible for the process of degradation (Glickman and Ciechanover, 2002; Lecker et al., 2006). Eukaryotic translation initiation factor 2B, subunit 1 alpha (EIF2B1) encodes one of five subunits of eukaryotic initiation factor 2B (eIF2B), which is known to affect global rates of protein synthesis. If protein synthesis is impaired, then eIF2B expression also is reduced as part of the ubiquitin-proteasome system based on studies in rat muscle (Tuckow et al., 2011). Mex-3 RNA binding family member C (MEX3C) is also an E3 ubiquitin-protein ligase that plays a role in ubiquitination of proteins, immune response, growth, and adiposity (Jiao et al., 2012a,b; Han et al., 2014; Kuniyoshi et al., 2014). In addition, many other DH transcripts were enriched for ubiquitin-proteasome components or functions (Supplementary Tables C–V). Growth factor receptor-bound protein 10 (GBR10) is an intracellular adaptor protein that may act as a growth suppressor, and its deletion causes hypermuscularity in mice (Holt et al., 2012). In KEGG pathways, the negatively correlated transcript targets of GBR10 enriched for proteasome and pyrimidine metabolism, while positively correlated transcripts enriched for endocytosis and ubiquitin-mediated proteolysis. CDK5 regulatory subunit associated protein 3 (CDK5RAP3) plays a role in cell signaling pathways, including transcriptional regulation and DNA damage response. CDK5RAP3 might be more or less susceptible to the ubiquitin proteasome system-mediated protein degradation depending on its interaction with other proteins and regulators (Wu et al., 2010). The positively correlated transcript targets of CDK5RAP3 enriched for the proteasome pathway, containing ATP-dependent degradation of ubiquitinated, 20S proteasome core, and 19S regulatory complex (RC) components, which are required for ubiquitin-dependent proteolysis and peptide cleavage (Glickman and Ciechanover, 2002) (Supplementary Tables C–V).

Besides TRIM45, USP2, RTP4, GRB10, HECTD4, EIF2B1, MEX3C, and CDK5RAP3 whose targets were enriched for functions related to the proteasome, other target transcripts correlated with other DH transcripts were identified, including: ALPK3, transmembrane BAX inhibitor motif containing 4 (TMBIM4), receptor (chemosensory) transporter protein 4 (RTP4), syntaxin binding protein 6 (amisyn) (STXBP6), transmembrane protein 150a (TMEM150A), ectonucleotide pyrophosphatase/phosphodiesterase 4 (putative) (ENPP4), and solute carrier family 25, member 44 (SLC25A44) (Supplementary Tables C–V). The eight hub genes presented in Figures 1, 2 are protein-coding genes according to Biomart Ensembl annotation, and most of them are linked to proteasome in the literature, except for two (TMBIM4 and TMEM150A). Many of the DH transcripts correlated with the proteasome KEGG pathway were not previously identified as associated with beef tenderness and thus might be involved in the protein degradation process of tenderization.

Two hub genes were related to calcium-activated mechanisms: ANO1 and TMBIM4. ANO1 is a calcium-activated chloride channel, which is essential for numerous physiological functions (Caputo et al., 2008). This gene could be a potential regulator of beef tenderness as the infusion of calcium chloride in meat increases tenderness (Koohmaraie, 1994). Calcium chloride injections in muscle increased the degradation of filamin (a large actin-binding protein), stimulating proteolysis and tenderization (Uytterhaegen et al., 1994; Harris et al., 2001; Lonergan et al., 2010). TMBIM4 is an anti-apoptotic protein that modulates calcium entry, promotes calcium release from intracellular reserves, and plays a role in cell survival (Gubser et al., 2007). TMBIM4 is expressed in several tissues including skeletal muscle in humans and also encodes a protein that contains transmembrane domains (Lisak et al., 2015). This gene is known as Golgi anti-apoptotic proteins (GAAP) in humans and may activate calpain 2, promote cell adhesion, and migration by controlling calcium concentration (Saraiva et al., 2013).

Differential hubbing analysis indicated that many potential regulators of beef tenderness are involved in the protein degradation process, and a few DH transcripts are related to calcium-activated mechanisms. Most of these transcripts were not identified in previous studies of SF, and offer novel insights into future research on SF and meat quality.

Co-expression Analysis Indicates microRNAs Impact Beef Tenderness

In muscle, miRNAs have been suggested to play a role in control gene expression impacting variation in beef tenderness (Zhao et al., 2012a). Several intriguing candidate transcripts were identified as possible regulators of variation in SF in this study, including bta-mir-133a-2 and bta-mir-22. The methodology used for library preparation and sequencing was a standard protocol for mRNA analysis. Although this study was not designated to find miRNAs, the procedure allowed the identification of two miRNAs precursors. Most miRNAs are transcribed from long primary miRNA (pri-miRNA) transcripts, which contains the mature miRNA as part of a predicted RNA hairpin (Zeng et al., 2005). The pri-miRNAs are then cleaved into ∼70 nt hairpin precursor miRNAs that might be able to rescue the miRNA function (Zeng et al., 2005; Van Wynsberghe et al., 2011). In agreement with this finding, other authors have also identified miRNAs using standard library preparation and sequencing for RNA-Seq analysis (Armenta-Medina et al., 2017; Lepe-Soltero et al., 2017). In addition, one miRNA was identified in the PIF analysis in a study of gene co-expression analysis, using the same RNA-Seq protocol, showing that it is possible to identify miRNAs or miRNA precursors even if the experiment is not designed to sequence miRNAs (Cesar et al., 2015). The bta-mir-133a-2 and bta-mir-22 were only identified in to the RIF analysis, and were among the top 10 potential regulators (based on extreme RIF Z score). The bta-mir-133a-2 was identified as both an extreme positive RIF1 and extreme negative RIF2 candidate regulatory gene. For this reason, these two miRNAs might be crucial regulators, as they have many target genes and pathways that can be involved in the control of meat tenderness. Functional analysis of miRNAs targets indicated that pathways related to calcium signaling, proteolysis, and apoptosis are impacted by fluctuations in bta-mir-133a-2 and bta-mir-22. In addition, both RIF and PIF co-expression analysis corroborate with the DH transcripts identified, indicating that protein degradation and calcium activity are key mechanisms to understand the differences between the H and L SF groups.

RIF1 results identified ubiquitin-conjugating enzyme E2D1 (UBE2D1). This gene promotes the ubiquitination of the 26S proteasome in cells (Jacobson et al., 2014). UBE2D1 may regulate SF variability since this gene encodes for a protein involved in proteasome function. Other biologically interesting RIF1 candidates identified included: histone cluster 1, H2ac (HIST1H2AC); synuclein, gamma (breast cancer-specific protein 1) (SNCG); and mohawk homeobox (MKX). Histones are proteins that play a role in transcriptional regulation, chromatin organization, and DNA packaging (Keogh et al., 2017). Genes coding for histone proteins like HIST1H2AC were up-regulated in skeletal muscle tissue of cattle expressing compensatory growth in Holstein Friesian bulls, suggesting an increase in feed efficiency in cattle (Keogh et al., 2017). SNCG increases metastasis at the cellular level and regulates many pathways including cell motility (Miao et al., 2014). MKX is a potential transcriptional repressor and critical regulator of gene expression that can repress myogenic differentiation (MYOD) gene expression in mice and negatively regulate muscle differentiation (Anderson et al., 2012; Chuang et al., 2014). Reduced expression of MYOD in mice muscle might be associated with tenderness (Lian et al., 2013). Interesting RIF2 candidates identified included: ATPase family, AAA domain containing 3A (ATAD3A); interferon induced transmembrane protein 1 (IFITM1); and TAF5-like RNA polymerase II, p300/CBP-associated factor (PCAF)-associated factor, 65kDa (TAF5L). ATAD3 is a member of the ATPases family, a highly conserved ubiquitous protein associated with development and mitochondrial DNA integrity (Issop et al., 2015). IFITM1 impacts cell growth. A transcriptomic analysis combined with RIF was performed in Blond d’Aquitaine and Charolais during pre-natal development and identified interferon signaling molecules as upregulated in the Blond d’Aquitaine breed compared to Charolais (Cassar-Malek et al., 2017). Interferon signaling has been suggested to play a role in the formation of new muscle fibers (Cassar-Malek et al., 2017). TAF5L is involved in the regulation of myogenic transcription, and MYOD function is enhanced by PCAF and the transcriptional coactivators p300 (Kuninger et al., 2006). Thus, it is plausible that these genes may impact muscle growth and subsequently tenderness.

In other species, miR-133a can promote muscle regeneration in skeletal muscle and regulates cardiac and muscle remodeling (Roderburg et al., 2013). The miR-133a is associated with fibrosis in liver and its overexpression can decrease the expression of collagens, which might indicate similar function in muscle (Roderburg et al., 2013). Both fibrosis and collagen expression could impact muscle tenderness. MiR-133a is a muscle-specific miRNA, or myomiR, which plays an important role in skeletal muscle proliferation and regeneration (Horak et al., 2016). Among the roles of myomiRs, miR-133a-2 promotes myoblast proliferation, regulates alternative splicing, is pro-apoptotic, controls mitochondrial metabolism, and impacts muscle fiber type (Horak et al., 2016). In addition, miR-133 is negatively regulated by mechanistic target of rapamycin (serine/threonine kinase) (MTOR), which regulates apoptosis, plays a role in myogenesis, and might control cell proliferation and differentiation (Horak et al., 2016).

In humans, miR-22, also a myomiR, impacts cell proliferation and apoptosis (Bersani et al., 2016). The role of bta-mir-22 has not been well characterized in skeletal muscle, but it plays a role in calcium signaling in the heart (Gurha et al., 2012). Consistent with the known function of bta-mir-22, target transcripts identified in this study were involved in calcium signaling. Transcript targets for both bta-mir-22 and bta-mir-133a-2 included calcium transporters such as the calcium channel, voltage-dependent, N type, alpha 1B subunit (CACNA1B) and ATPase, calcium-transporting, plasma membrane 2 (ATP2B2). Transcripts activated by calcium or calcium-binding such as calcium/calmodulin-dependent protein kinase ID (CAMK1D), calcium/calmodulin-dependent protein kinase kinase 2, beta (CAMKK2), and calcium-dependent secretion activator 2 (CADPS2) were targets of both miRNAs. Calpain 5 (CAPN5) and calpain 15 (CAPN15) were targets of bta-mir-133a-2, while calpain 7 (CAPN7) was target of bta-mir-22. Interestingly, SYT4 identified in the differential expression analysis was also targeted by bta-mir-22.

Apoptosis is known to be an important process for beef tenderness and to be involved in the aging process of meat (Ouali et al., 2006). Both miRNAs had many transcript targets that are anti- or pro-apoptotic, such as Bcl2-like 1 (BCL2L1) that inhibits cell death, and Bcl2-like 13 (apoptosis facilitator) (BCL2L13), which promotes apoptosis (Kataoka et al., 2001). Other transcript targets identified in the RIF co-expression analysis included programmed cell death 6 interacting protein (PDCD6IP) and programmed cell death 10 (PDCD10). In agreement with these findings, a protein study of SF in Nellore cattle identified heat shock protein 27 (HSP27) and heat shock protein 70 (HSP70) as having anti-apoptotic functions (Carvalho et al., 2014). However, the two proteins had differing effects. While HSP27 may contribute to increased tenderness, HSP70 is suggested to negatively impact tenderness (Carvalho et al., 2014).

Among the top 10 PIF transcripts, myoglobin (MB) had the highest ranking, which implicates it as the most important transcript related to the SF variation based on this method. MB is a cytoplasmic hemoprotein, which can serve as an oxygen reservoir releasing oxygen during periods of hypoxia or anoxia (Ordway and Garry, 2004). It can keep the intracellular concentration of oxygen constant even in high levels of activity (Ordway and Garry, 2004). MB has a similar role to that of creatine phosphokinase, which increases ATP concentrations during muscle activity; however, MB increases the oxygen concentrations instead of ATP (Ordway and Garry, 2004). Oxidative energy metabolism is important to the transformation of muscle into meat. When cells are deprived of oxygen, ATP synthesis occurs via anaerobic conditions with glycolytic phosphorylation by creatine (Ouali et al., 2006). In a study of beef tenderness in Chianina cattle, expression of creatine was higher in the tender group (D’Alessandro et al., 2012). MB was also a target of bta-mir-22, present in high quantities in type I muscle fibers, and responsible for oxygen storage and transport (Picard et al., 2011). Type I fibers were associated with beef tenderness, as increased type I fibers results in more tender beef (Strydom et al., 2000; Renand et al., 2001; Choi and Kim, 2009). In addition, myoglobin is involved in iron ion binding. The presence of iron may impact beef tenderness as markers in the Cast and mu-calpain (CAPN4751) genes were associated with iron levels (Casas et al., 2014). Thus, changes in myoglobin may act as a surrogate to track differences in fiber type between the extreme tenderness samples used in this study.

High glycolytic metabolism and oxidative metabolism are dependent on muscle type and breed (Cassar-Malek and Picard, 2016). For example, enolase 3 (beta, muscle) (ENO3) was positively correlated with tenderness in glycolytic muscles, like semitendinosus, and in breeds that present glycolytic muscle properties like Limousin (Cassar-Malek and Picard, 2016). In another study, two isoforms of the glycolytic enzyme enolase (ENO1 and ENO3) were associated with tenderness in Charolais (Guillemin et al., 2011). However, ENO1 and ENO3 have different biological functions. ENO1 (alpha, muscle) is a transcriptional regulator involved in the immunoglobulin response, while ENO3 (beta, muscle) is involved in the development and regeneration of muscle. Also, ENO3 was more abundant than ENO1 in steers compared to young bulls (Guillemin et al., 2011). Among the top PIF transcripts, ENO3 had the second highest score, and carbonic anhydrase III (CA3) had the third highest score in the ranking. CA3 is a high-speed enzyme involved in proton homeostasis and converting carbohydrate into carbon dioxide. Changes in proton homeostasis decrease the electrostatic repulsion of myofibrillar proteins and facilitate the lateral shrinkage of muscle fibers (sarcomere length) (Pearce et al., 2011). The high presence of the glycolytic enzyme CA3 is also associated with increased meat tenderness (Laville et al., 2009; Choi et al., 2010; D’Alessandro et al., 2012). Two other transcripts from the CA family were targets of bta-mir-22: carbonic anhydrase VIII (CA8) and carbonic anhydrase XII (CA12). Results found in the RIF and PIF co-expression analysis identified two miRNAs as potential regulators and also confirmed that some transcripts like MB, ENO3, and CA3, previously described by other studies, are good markers of tenderness.

Proteolysis impacts the structure of the sarcomere/myofibril (Ertbjerg and Puolanne, 2017). The thin filaments are made of the protein actin, while the thick filaments are composed by the protein myosin (Ertbjerg and Puolanne, 2017). PIF analysis identified myosin-binding protein C, fast type (MYBPC2) as a candidate biomarker for tenderness. Myosin is an important structural protein of the thick filament of the sarcomere, responsible for contraction velocity and power (Zhao et al., 2012b). Both myosin and actin interact for contraction to occur, forming the complex called actomyosin that is ATP dependent, which means that the actin and myosin can dissociate in the presence of ATP (Goll et al., 1984). Actomyosin bonds become irreversible in postmortem muscles when the supply of ATP is depleted. Myosin heavy chains (MyH3) are present in the tender meat due to the proteolytic activity of the muscle that continues postmortem (Zhao et al., 2012b). The amount that a myofibril shortens affects tenderness, and meat with very short sarcomeres might be tougher (Locker and Hagyard, 1963). Myosin regulatory light chain (MLC2) and myosin binding protein H (MyBP-H) expression are also associated with meat tenderness (Carvalho et al., 2014). In addition, histone methylation and modifications are related to beef tenderness in Angus based on a study in fast and slow skeletal muscle fibers, which identified differential histone modification patterns in myosin heavy chain genes (Zhao et al., 2015).

Conclusion

Differential gene expression analysis identified transcripts that could be responsible for variation in beef tenderness in Nellore cattle, which differ from those previously identified in Bos taurus breeds. Systems biology analyses using PCIT, PIF and RIF helped to elucidate differential gene networks that may be associated with beef tenderness, many of which were related to proteolysis, calcium signaling, and apoptosis. Furthermore, two microRNAs, bta-mir-133a-2 and bta-mir-22, were identified in the RIF analysis as potential regulators of tenderness. These miRNA were associated with proteolysis and apoptosis as well as a number other genes identified in this study to impact tenderness. While most of the genes identified in this study have not been previously implicated in beef tenderness, we conclude that the identified pathways comprised by these genes were consistent and interconnected with previously identified genes that regulate beef tenderness. These findings provide a better understanding of the molecular mechanisms involved in the regulation of beef tenderness in Nellore cattle.

Data Availability Statement

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

Author Contributions

TG performed the bioinformatics and data analysis, and wrote the manuscript. LdAR co-designed the experiment and co-directed the work. JK, AC, and SdSA contributed to the bioinformatics and data analyses. GBM co-designed the experiment. GG contributed to the sequencing analysis. GCM revised all the drafts of the manuscript. EF-W contributed to the data analysis. JR supervised the data analysis. LC supervised the data analysis, co-designed the experiment, and co-directed the work. All the authors have read and approved the final manuscript.

Funding

This study was conducted with funding from EMBRAPA (Macroprogram 1, 01/2005) and FAPESP (Processes No. 2012/23934-6 and 2013/26121-9). This study was also financed in part by the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior – Brazil (CAPES) – Finance Code 001.

Conflict of Interest Statement

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

Acknowledgments

We are thankful to Marcelo M. Brandão and Horácio Montenegro from the Laboratório Molecular de Plantas, ESALQ – USP for the access to the Thunder Server, to Rymer Ramiz Tullio and Renata Tieko Nassu for phenotype collection. LdAR, GBM, and LC are recipient of productivity fellowship from CNPq. This manuscript is part of the first author’s Master’s thesis cited as Gonçalves (2015).

Supplementary Material

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

Abbreviations

DE, differentially expressed; DH, differentially hubbed; EBVSF, shear force estimated breeding values; H, highest values of EBVSF; L, lowest values of EBVSF; LD, Longissimus dorsi; PCIT, partial correlation with information theory; PIF, phenotypic impact factor; RIF, regulatory impact factor; SF, shear force.

References

Abo-Ismail, M., Kelly, M., Squires, E., Swanson, K., Bauck, S., and Miller, S. (2013). Identification of single nucleotide polymorphisms in genes involved in digestive and metabolic processes associated with feed efficiency and performance traits in beef cattle. J. Anim. Sci. 91, 2512–2529. doi: 10.2527/jas.2012-5756

PubMed Abstract | CrossRef Full Text | Google Scholar

Aganna, E., Burrin, J. M., Hitman, G. A., and Turner, M. D. (2006). Involvement of calpain and synaptotagmin Ca²⁺ sensors in hormone secretion from excitable endocrine cells. J. Endocrinol. 190, R1–R7.

Google Scholar

Alarcon-Rojo, A., and Dransfield, E. (1995). Alteration of post-mortem ageing in beef by the addition of enzyme inhibitors and activators. Meat Sci. 41, 163–178. doi: 10.1016/0309-1740(95)99781-V

PubMed Abstract | CrossRef Full Text | Google Scholar

Altschul, S. F., Madden, T. L., Schäffer, A. A., Zhang, J., Zhang, Z., Miller, W., et al. (1997). Gapped BLAST and PSI-BLAST: a new generation of protein database search programs. Nucleic Acids Res. 25, 3389–3402. doi: 10.1093/nar/25.17.3389

PubMed Abstract | CrossRef Full Text | Google Scholar

Anders, S., and Huber, W. (2010). Differential expression analysis for sequence count data. Genome Biol. 11:R106. doi: 10.1186/gb-2010-11-10-r106

PubMed Abstract | CrossRef Full Text | Google Scholar

Anderson, D. M., George, R., Noyes, M. B., Rowton, M., Liu, W., Jiang, R., et al. (2012). Characterization of the DNA-binding properties of the Mohawk homeobox transcription factor. J. Biol. Chem. 287, 35351–35359. doi: 10.1074/jbc.M112.399386

PubMed Abstract | CrossRef Full Text | Google Scholar

Andrews, S. (2010). FastQC: a Quality Control Tool for High Throughput Sequence Data. Cambridge: Babraham Institute.

Google Scholar

Armenta-Medina, A., Lepe-Soltero, D., Xiang, D., Datla, R., Abreu-Goodger, C., and Gillmor, C. S. (2017). Arabidopsis thaliana miRNAs promote embryo pattern formation beginning in the zygote. Dev. Biol. 431, 145–151. doi: 10.1016/j.ydbio.2017.09.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Bai, J., Lai, L., Yeo, H. C., Goh, B. C., and Tan, T. M. (2004). Multidrug resistance protein 4 (MRP4/ABCC4) mediates efflux of bimane-glutathione. Int. J. Biochem. Cell Biol. 36, 247–257. doi: 10.1016/S1357-2725(03)00236-X

CrossRef Full Text | Google Scholar

Barabasi, A. L., and Oltvai, Z. N. (2004). Network biology: understanding the cell’s functional organization. Nat. Rev. Genet. 5, 101–113. doi: 10.1038/nrg1272

PubMed Abstract | CrossRef Full Text | Google Scholar

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

Google Scholar

Bersani, F., Lingua, M. F., Morena, D., Foglizzo, V., Miretti, S., Lanzetti, L., et al. (2016). Deep sequencing reveals a novel miR-22 regulatory network with therapeutic potential in rhabdomyosarcoma. Cancer Res. 76, 6095–6106. doi: 10.1158/0008-5472.CAN-16-0709

PubMed Abstract | CrossRef Full Text | Google Scholar

Betel, D., Wilson, M., Gabow, A., Marks, D. S., and Sander, C. (2008). The microRNA. org resource: targets and expression. Nucleic Acids Res. 36, D149–D153.

Google Scholar

Bongiorni, S., Gruber, C., Bueno, S., Chillemi, G., Ferre, F., Failla, S., et al. (2016). Transcriptomic investigation of meat tenderness in two Italian cattle breeds. Anim. Genet. 47, 273–287. doi: 10.1111/age.12418

PubMed Abstract | CrossRef Full Text | Google Scholar

Boudida, Y., Gagaoua, M., Becila, S., Picard, B., Boudjellal, A., Herrera-Mendez, C. H., et al. (2016). Serine protease inhibitors as good predictors of meat tenderness: which are they and what are their functions? Crit. Rev. Food Sci. Nutr. 56, 957–972. doi: 10.1080/10408398.2012.741630

PubMed Abstract | CrossRef Full Text | Google Scholar

Bullard, J. H., Purdom, E., Hansen, K. D., and Dudoit, S. (2010). Evaluation of statistical methods for normalization and differential expression in mRNA-Seq experiments. BMC Bioinformatics 11:94. doi: 10.1186/1471-2105-11-94

PubMed Abstract | CrossRef Full Text | Google Scholar

Burden, C. J., Qureshi, S. E., and Wilson, S. R. (2014). Error estimates for the analysis of differential expression from RNA-seq count data. PeerJ 2:e576. doi: 10.7717/peerj.576

PubMed Abstract | CrossRef Full Text | Google Scholar

Calkins, C. R., and Sullivan, G. (2007). Adding Enzymes to Improve Beef Tenderness. Beef Facts Product Enhancement, National Cattleman’s Beef Association. Centennial: Cattlemen’s Beef Board.

Google Scholar

Caputo, A., Caci, E., Ferrera, L., Pedemonte, N., Barsanti, C., Sondo, E., et al. (2008). TMEM16A, a membrane protein associated with calcium-dependent chloride channel activity. Science 322, 590–594. doi: 10.1126/science.1163518

PubMed Abstract | CrossRef Full Text | Google Scholar

Carvalho, M. E., Gasparin, G., Poleti, M. D., Rosa, A. F., Balieiro, J. C. C., Labate, C. A., et al. (2014). Heat shock and structural proteins associated with meat tenderness in Nellore beef cattle, a Bos indicus breed. Meat Sci. 96, 1318–1324. doi: 10.1016/j.meatsci.2013.11.014

PubMed Abstract | CrossRef Full Text | Google Scholar

Casas, E., Duan, Q., Schneider, M., Shackelford, S., Wheeler, T. L., Cundiff, L. V., et al. (2014). Polymorphisms in calpastatin and mu-calpain genes are associated with beef iron content. Anim. Genet. 45, 283–284. doi: 10.1111/age.12108

PubMed Abstract | CrossRef Full Text | Google Scholar

Cassar-Malek, I., Boby, C., Picard, B., Reverter, A., and Hudson, N. J. (2017). Molecular regulation of high muscle mass in developing Blonde d’Aquitaine cattle fetuses. Biol. Open 6, 1483–1492. doi: 10.1242/bio.024950

PubMed Abstract | CrossRef Full Text | Google Scholar

Cassar-Malek, I., and Picard, B. (2016). Expression marker-based strategy to improve beef quality. Sci. World J. 2016:2185323. doi: 10.1155/2016/2185323

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Choi, Y., and Kim, B. (2009). Muscle fiber characteristics, myofibrillar protein isoforms, and meat quality. Livest. Sci. 122, 105–118. doi: 10.1016/j.livsci.2008.08.015

CrossRef Full Text | Google Scholar

Choi, Y., Lee, S., Choe, J., Rhee, M., Lee, S., Joo, S., et al. (2010). Protein solubility is related to myosin isoforms, muscle fiber types, meat quality traits, and postmortem protein changes in porcine longissimus dorsi muscle. Livest. Sci. 127, 183–191. doi: 10.1016/j.livsci.2009.09.009

CrossRef Full Text | Google Scholar

Chuang, H. N., Hsiao, K. M., Chang, H. Y., Wu, C. C., and Pan, H. (2014). The homeobox transcription factor Irxl1 negatively regulates MyoD expression and myoblast differentiation. FEBS J. 281, 2990–3003. doi: 10.1111/febs.12837

PubMed Abstract | CrossRef Full Text | Google Scholar

Cunningham, F., Amode, M. R., Barrell, D., Beal, K., Billis, K., Brent, S., et al. (2014). Ensembl 2015. Nucleic Acids Res. 43, D662–D669. doi: 10.1093/nar/gku1010

PubMed Abstract | CrossRef Full Text | Google Scholar

D’Alessandro, A., Marrocco, C., Rinalducci, S., Mirasole, C., Failla, S., and Zolla, L. (2012). Chianina beef tenderness investigated through integrated Omics. J. Proteomics 75, 4381–4398. doi: 10.1016/j.jprot.2012.03.052

PubMed Abstract | CrossRef Full Text | Google Scholar

De Bruin, E., Meersma, D., De Wilde, J., Den Otter, I., Schipper, E., Medema, J., et al. (2003). A serine protease is involved in the initiation of DNA damage-induced apoptosis. Cell Death Differ. 10, 1204–1212. doi: 10.1038/sj.cdd.4401296

PubMed Abstract | CrossRef Full Text | Google Scholar

de la Fuente, A. (2010). From ‘differential expression’to ‘differential networking’–identification of dysfunctional regulatory networks in diseases. Trends Genet. 26, 326–333. doi: 10.1016/j.tig.2010.05.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Di Cera, E. (2009). Serine proteases. IUBMB Life 61, 510–515. doi: 10.1002/iub.186

PubMed Abstract | CrossRef Full Text | Google Scholar

Drag, M., Skinkyté-Juskiené, R., Do, D. N., Kogelman, L. J., and Kadarmideen, H. N. (2017). Differential expression and co-expression gene networks reveal candidate biomarkers of boar taint in non-castrated pigs. Sci. Rep. 7:12205. doi: 10.1038/s41598-017-11928-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Enright, A. J., John, B., Gaul, U., Tuschl, T., Sander, C., and Marks, D. S. (2003). MicroRNA targets in Drosophila. Genome Biol. 5:R1.

Google Scholar

Enright, A. J., and Ouzounis, C. A. (2001). BioLayout—an automatic graph layout algorithm for similarity visualization. Bioinformatics 17, 853–854. doi: 10.1093/bioinformatics/17.9.853

CrossRef Full Text | Google Scholar

Ertbjerg, P., and Puolanne, E. (2017). Muscle structure, sarcomere length and influences on meat quality: a review. Meat Sci. 132, 139–152. doi: 10.1016/j.meatsci.2017.04.261

PubMed Abstract | CrossRef Full Text | Google Scholar

Fonseca, L. F. S., Gimenez, D. F. J., dos Santos Silva, D. B., Barthelson, R., Baldi, F., Ferro, J. A., et al. (2017). Differences in global gene expression in muscle tissue of Nellore cattle with divergent meat tenderness. BMC Genomics 18:945. doi: 10.1186/s12864-017-4323-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Fukuda, M., and Itoh, T. (2004). Slac2-a/melanophilin contains multiple PEST-like sequences that are highly sensitive to proteolysis. J. Biol. Chem. 279, 22314–22321. doi: 10.1074/jbc.M401791200

PubMed Abstract | CrossRef Full Text | Google Scholar

Glavan, G., Schliebs, R., and Živin, M. (2009). Synaptotagmins in neurodegeneration. Anat. Rec. 292, 1849–1862. doi: 10.1002/ar.21026

PubMed Abstract | CrossRef Full Text | Google Scholar

Glickman, M. H., and Ciechanover, A. (2002). The ubiquitin-proteasome proteolytic pathway: destruction for the sake of construction. Physiol. Rev. 82, 373–428. doi: 10.1152/physrev.00027.2001

PubMed Abstract | CrossRef Full Text | Google Scholar

Goll, D. E., Robson, R. M., and Stromer, M. H. (1984). “Skeletal muscle, nervous system, temperature regulation, and special senses,” in Dukes’ Physiology of Domestic Animals, 10th Edn, ed. M. J. Swenson (Ithaca, NY: Cornell University Press), 548–580.

Google Scholar

Gonçalves, T. M. (2015). Differential Expression of Genes Related with Meat Tenderness in Nellore Cattle. Master thesis, University of São Paulo, Piracicaba.

Google Scholar

Graifer, D., and Karpova, G. (2015). Roles of ribosomal proteins in the functioning of translational machinery of eukaryotes. Biochimie 109, 1–17. doi: 10.1016/j.biochi.2014.11.016

PubMed Abstract | CrossRef Full Text | Google Scholar

Gubser, C., Bergamaschi, D., Hollinshead, M., Lu, X., van Kuppeveld, F. J., and Smith, G. L. (2007). A new inhibitor of apoptosis from vaccinia virus and eukaryotes. PLoS Pathog. 3:e17. doi: 10.1371/journal.ppat.0030017

PubMed Abstract | CrossRef Full Text | Google Scholar

Guillemin, N., Jurie, C., Cassar-Malek, I., Hocquette, J. F., Renand, G., and Picard, B. (2011). Variations in the abundance of 24 protein biomarkers of beef tenderness according to muscle and animal type. Animal 5, 885–894. doi: 10.1017/S1751731110002612

PubMed Abstract | CrossRef Full Text | Google Scholar

Gurha, P., Abreu-Goodger, C., Wang, T., Ramirez, M. O., Drummond, A. L., Van Dongen, S., et al. (2012). Targeted deletion of microRNA-22 promotes stress induced cardiac dilation and contractile dysfunction. Circulation 25, 2751–2761. doi: 10.1161/CIRCULATIONAHA.111.044354

PubMed Abstract | CrossRef Full Text | Google Scholar

Hammond, C. L., Marchan, R., Krance, S. M., and Ballatori, N. (2007). Glutathione export during apoptosis requires functional multidrug resistance-associated proteins. J. Biol. Chem. 282, 14337–14347. doi: 10.1074/jbc.M611019200

PubMed Abstract | CrossRef Full Text | Google Scholar

Han, C., Jiao, Y., Zhao, Q., and Lu, B. (2014). Mex3c mutation reduces adiposity partially through increasing physical activity. J. Endocrinol. 221, 457–468. doi: 10.1530/JOE-14-0071

PubMed Abstract | CrossRef Full Text | Google Scholar

Harris, S., Huff-Lonergan, E., Lonergan, S. M., Jones, W., and Rankins, D. (2001). Antioxidant status affects color stability and tenderness of calcium chloride-injected beef. J. Anim. Sci. 79, 666–677. doi: 10.2527/2001.793666x

PubMed Abstract | CrossRef Full Text | Google Scholar

Henderson, C. R. (1973). Sire evaluation and genetic trends. J. Anim. Sci. 1973, 10–41. doi: 10.1093/ansci/1973.Symposium.10

CrossRef Full Text | Google Scholar

Hill, W. G., Goddard, M. E., and Visscher, P. M. (2008). Data and theory point to mainly additive genetic variance for complex traits. PLoS Genet. 4:e1000008. doi: 10.1371/journal.pgen.1000008

PubMed Abstract | CrossRef Full Text | Google Scholar

Holt, L. J., Turner, N., Mokbel, N., Trefely, S., Kanzleiter, T., Kaplan, W., et al. (2012). Grb10 regulates the development of fiber number in skeletal muscle. FASEB J. 26, 3658–3669. doi: 10.1096/fj.11-199349

PubMed Abstract | CrossRef Full Text | Google Scholar

Horak, M., Novak, J., and Bienertova-Vasku, J. (2016). Muscle-specific microRNAs in skeletal muscle development. Dev. Biol. 410, 1–13. doi: 10.1016/j.ydbio.2015.12.013

PubMed Abstract | CrossRef Full Text | Google Scholar

Huang, D. W., Sherman, B. T., and Lempicki, R. A. (2008). Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat. Protoc. 4, 44–57. doi: 10.1038/nprot.2008.211

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Hui, Y., Yang, Y., Liu, R., Zhang, Y., Xiang, C., Liu, Z., et al. (2013). Significant association of APOA5 and APOC3 gene polymorphisms with meat quality traits in Kele pigs. Genet. Mol. Res. 12, 3643–3650. doi: 10.4238/2013.September.13.8

PubMed Abstract | CrossRef Full Text | Google Scholar

Issop, L., Fan, J., Lee, S., Rone, M. B., Basu, K., Mui, J., et al. (2015). Mitochondria-associated membrane formation in hormone-stimulated Leydig cell steroidogenesis: role of ATAD3. Endocrinology 156, 334–345. doi: 10.1210/en.2014-1503

PubMed Abstract | CrossRef Full Text | Google Scholar

Jacobson, A. D., MacFadden, A., Wu, Z., Peng, J., and Liu, C. W. (2014). Autoregulation of the 26S proteasome by in situ ubiquitination. Mol. Biol. Cell 25, 1824–1835. doi: 10.1091/mbc.E13-10-0585

PubMed Abstract | CrossRef Full Text | Google Scholar

Jiao, Y., Bishop, C. E., and Lu, B. (2012a). Mex3c regulates insulin-like growth factor 1 (IGF1) expression and promotes postnatal growth. Mol. Biol. Cell 23, 1404–1413. doi: 10.1091/mbc.E11-11-0960

PubMed Abstract | CrossRef Full Text | Google Scholar

Jiao, Y., George, S. K., Zhao, Q., Hulver, M. W., Hutson, S. M., Bishop, C. E., et al. (2012b). Mex3c mutation reduces adiposity and increases energy expenditure. Mol. Cell. Biol. 32, 4350–4362. doi: 10.1128/MCB.00452-12

PubMed Abstract | CrossRef Full Text | Google Scholar

Johnson, S. L., Franz, C., Kuhn, S., Furness, D. N., Rüttiger, L., Münkner, S., et al. (2010). Synaptotagmin IV determines the linear Ca²⁺ dependence of vesicle fusion at auditory ribbon synapses. Nat. Neurosci. 13, 45–52. doi: 10.1038/nn.2456

PubMed Abstract | CrossRef Full Text | Google Scholar

Kadarmideen, H. N., von Rohr, P., and Janss, L. L. (2006). From genetical genomics to systems genetics: potential applications in quantitative genomics and animal breeding. Mamm. Genome 17, 548–564. doi: 10.1007/s00335-005-0169-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Kataoka, T., Holler, N., Micheau, O., Martinon, F., Tinel, A., Hofmann, K., et al. (2001). Bcl-rambo, a novel Bcl-2 homologue that induces apoptosis via its unique C-terminal extension. J. Biol. Chem. 276, 19548–19554. doi: 10.1074/jbc.M010520200

PubMed Abstract | CrossRef Full Text | Google Scholar

Keogh, K., Waters, S. M., Cormican, P., Kelly, A. K., O’Shea, E., and Kenny, D. A. (2017). Effect of dietary restriction and subsequent re-alimentation on the transcriptional profile of bovine ruminal epithelium. PLoS One 12:e0177852. doi: 10.1371/journal.pone.0177852

PubMed Abstract | CrossRef Full Text | Google Scholar

Koch, R., Cundiff, L., and Gregory, K. (1982). Heritabilities and genetic, environmental and phenotypic correlations of carcass traits in a population of diverse biological types and their implications in selection programs 1, 2. J. Anim. Sci. 55, 1319–1329. doi: 10.2527/jas1982.5561319x

CrossRef Full Text | Google Scholar

Koesterke, L., Milfeld, K., Vaughn, M. W., Stanzione, D., Koltes, J. E., Weeks, N. T., et al. (2013). Optimizing the PCIT algorithm on stampede’s Xeon and Xeon Phi processors for faster discovery of biological networks. Paper Presented at the Proceedings of the Conference on Extreme Science and Engineering Discovery Environment: Gateway to Discovery, San Diego. doi: 10.1145/2484762.2484794

CrossRef Full Text | Google Scholar

Kogelman, L. J., Cirera, S., Zhernakova, D. V., Fredholm, M., Franke, L., and Kadarmideen, H. N. (2014). Identification of co-expression gene networks, regulatory genes and pathways for obesity based on adipose tissue RNA Sequencing in a porcine model. BMC Med. Genomics 7:57. doi: 10.1186/1755-8794-7-57

PubMed Abstract | CrossRef Full Text | Google Scholar

Koohmaraie, M. (1994). Muscle proteinases and meat aging. Meat Sci. 36, 93–104. doi: 10.1016/0309-1740(94)90036-1

CrossRef Full Text | Google Scholar

Kuninger, D., Wright, A., and Rotwein, P. (2006). Muscle cell survival mediated by the transcriptional coactivators p300 and PCAF displays different requirements for acetyltransferase activity. Am. J. Physiol. Cell Physiol. 291, C699–C709. doi: 10.1152/ajpcell.00056.2006

PubMed Abstract | CrossRef Full Text | Google Scholar

Kuniyoshi, K., Takeuchi, O., Pandey, S., Satoh, T., Iwasaki, H., Akira, S., et al. (2014). Pivotal role of RNA-binding E3 ubiquitin ligase MEX3C in RIG-I–mediated antiviral innate immunity. Proc. Natl. Acad. Sci. U.S.A. 111, 5646–5651. doi: 10.1073/pnas.1401674111

PubMed Abstract | CrossRef Full Text | Google Scholar

Laville, E., Sayd, T., Morzel, M., Blinet, S., Chambon, C., Lepetit, J., et al. (2009). Proteome changes during meat aging in tough and tender beef suggest the importance of apoptosis and protein solubility for beef aging and tenderization. J. Agric. Food Chem. 57, 10755–10764. doi: 10.1021/jf901949r

PubMed Abstract | CrossRef Full Text | Google Scholar

Lecker, S. H., Goldberg, A. L., and Mitch, W. E. (2006). Protein degradation by the ubiquitin–proteasome pathway in normal and disease states. J. Am. Soc. Nephrol. 17, 1807–1819. doi: 10.1681/ASN.2006010083

PubMed Abstract | CrossRef Full Text | Google Scholar

Lek, A., Evesson, F. J., Lemckert, F. A., Redpath, G. M., Lueders, A.-K., Turnbull, L., et al. (2013). Calpains, cleaved mini-dysferlinC72, and L-type channels underpin calcium-dependent muscle membrane repair. J. Neurosci. 33, 5085–5094. doi: 10.1523/JNEUROSCI.3560-12.2013

PubMed Abstract | CrossRef Full Text | Google Scholar

Lepe-Soltero, D., Armenta-Medina, A., Xiang, D., Datla, R., Gillmor, C. S., and Abreu-Goodger, C. (2017). Annotating and quantifying pri-miRNA transcripts using RNA-Seq data of wild type and serrate-1 globular stage embryos of Arabidopsis thaliana. Data Brief 15, 642–647. doi: 10.1016/j.dib.2017.10.019

PubMed Abstract | CrossRef Full Text | Google Scholar

Lian, T., Wang, L., and Liu, Y. (2013). A new insight into the role of calpains in post-mortem meat tenderization in domestic animals: a review. Asian-Australas. J. Anim. Sci. 26, 443–454. doi: 10.5713/ajas.2012.12365

PubMed Abstract | CrossRef Full Text | Google Scholar

Lisak, D. A., Schacht, T., Enders, V., Habicht, J., Kiviluoto, S., Schneider, J., et al. (2015). The transmembrane Bax inhibitor motif (TMBIM) containing protein family: tissue expression, intracellular localization and effects on the ER Ca²⁺ -filling state. Biochim. Biophys. Acta 1853, 2104–2114. doi: 10.1016/j.bbamcr.2015.03.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Locker, R. H., and Hagyard, C. J. (1963). A cold shortening effect in beef muscles. J. Sci. Food Agric. 14, 787–793. doi: 10.1002/jsfa.2740141103

CrossRef Full Text | Google Scholar

Lonergan, E. H., Zhang, W., and Lonergan, S. M. (2010). Biochemistry of postmortem muscle—lessons on mechanisms of meat tenderization. Meat Sci. 86, 184–195. doi: 10.1016/j.meatsci.2010.05.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Lund, S. P., Nettleton, D., McCarthy, D. J., and Smyth, G. K. (2012). “Detecting differential expression in RNA-sequence data using quasi-likelihood with shrunken dispersion estimates,” in Statistical Applications in Genetics and Molecular Biology, ed. M. P. H. Stumpf (Berlin: de Gruyter).

Google Scholar

Mateescu, R., Reecy, J., and Garrick, D. (2017). Network analysis reveals putative genes affecting meat quality in Angus cattle. Front. Genet. 8:171. doi: 10.3389/fgene.2017.00171

PubMed Abstract | CrossRef Full Text | Google Scholar

McCarthy, F. M., Gresham, C. R., Buza, T. J., Chouvarine, P., Pillai, L. R., Kumar, R., et al. (2010). AgBase: supporting functional modeling in agricultural organisms. Nucleic Acids Res. 39, D497–D506. doi: 10.1093/nar/gkq1115

PubMed Abstract | CrossRef Full Text | Google Scholar

Mellgren, R. L., Miyake, K., Kramerova, I., Spencer, M. J., Bourg, N., Bartoli, M., et al. (2009). Calcium-dependent plasma membrane repair requires m-or μ-calpain, but not calpain-3, the proteasome, or caspases. Biochim. Biophys. Acta 1793, 1886–1893. doi: 10.1016/j.bbamcr.2009.09.013

PubMed Abstract | CrossRef Full Text | Google Scholar

Metzker, M. L. (2010). Sequencing technologies—the next generation. Nat. Rev. Genet. 11, 31–46. doi: 10.1038/nrg2626

PubMed Abstract | CrossRef Full Text | Google Scholar

Miao, S., Wu, K., Zhang, B., Weng, Z., Zhu, M., Lu, Y., et al. (2014). Synuclein γ compromises spindle assembly checkpoint and renders resistance to antimicrotubule drugs. Mol. Cancer Ther. 13, 699–713. doi: 10.1158/1535-7163.MCT-13-0671

PubMed Abstract | CrossRef Full Text | Google Scholar

Mihalyi, E. (1953). Trypsin digestion of muscle proteins II. The kinetics of the digestion. J. Biol. Chem. 201, 197–209.

Google Scholar

Mrode, R. A. (2014). Linear Models for the Prediction of Animal Breeding Values, 2nd Edn, Wallingford: CAB International. doi: 10.1079/9781780643915.0000

CrossRef Full Text | Google Scholar

Nonneman, D., Shackelford, S., King, D., Wheeler, T., Wiedmann, R., Snelling, W., et al. (2013). Genome-wide association of meat quality traits and tenderness in swine. J. Anim. Sci. 91, 4043–4050. doi: 10.2527/jas.2013-6255

PubMed Abstract | CrossRef Full Text | Google Scholar

O’Brien, P. J., Alborn, W. E., Sloan, J. H., Ulmer, M., Boodhoo, A., Knierman, M. D., et al. (2005). The novel apolipoprotein A5 is present in human serum, is associated with VLDL, HDL, and chylomicrons, and circulates at very low concentrations compared with other apolipoproteins. Clin. Chem. 51, 351–359. doi: 10.1373/clinchem.2004.040824

PubMed Abstract | CrossRef Full Text | Google Scholar

Ordway, G. A., and Garry, D. J. (2004). Myoglobin: an essential hemoprotein in striated muscle. J. Exp. Biol. 207, 3441–3446. doi: 10.1242/jeb.01172

PubMed Abstract | CrossRef Full Text | Google Scholar

Ouali, A., Herrera-Mendez, C. H., Coulis, G., Becila, S., Boudjellal, A., Aubry, L., et al. (2006). Revisiting the conversion of muscle into meat and the underlying mechanisms. Meat Sci. 74, 44–58. doi: 10.1016/j.meatsci.2006.05.010

PubMed Abstract | CrossRef Full Text | Google Scholar

Patterson, H. D., and Thompson, R. (1971). Recovery of inter-block information when block sizes are unequal. Biometrika 58, 545–554. doi: 10.1093/biomet/58.3.545

CrossRef Full Text | Google Scholar

Pearce, K. L., Rosenvold, K., Andersen, H. J., and Hopkins, D. L. (2011). Water distribution and mobility in meat during the conversion of muscle to meat and ageing and the impacts on fresh meat quality attributes—A review. Meat Sci. 89, 111–124. doi: 10.1016/j.meatsci.2011.04.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Picard, M., Hepple, R. T., and Burelle, Y. (2011). Mitochondrial functional specialization in glycolytic and oxidative muscle fibers: tailoring the organelle for optimal function. Am. J. Physiol. Cell Physiol. 302, C629–C641. doi: 10.1152/ajpcell.00368.2011

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Renand, G., Picard, B., Touraille, C., Berge, P., and Lepetit, J. (2001). Relationships between muscle characteristics and meat quality traits of young Charolais bulls. Meat Sci. 59, 49–60. doi: 10.1016/S0309-1740(01)00051-1

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Reverter, A., Hudson, N. J., Wang, Y., Tan, S.-H., Barris, W., Byrne, K. A., et al. (2006). A gene coexpression network for bovine skeletal muscle inferred from microarray data. Physiol. Genomics 28, 76–83. doi: 10.1152/physiolgenomics.00105.2006

PubMed Abstract | CrossRef Full Text | Google Scholar

Roderburg, C., Luedde, M., Cardenas, D. V., Vucur, M., Mollnow, T., Zimmermann, H. W., et al. (2013). miR-133a mediates TGF-β-dependent derepression of collagen synthesis in hepatic stellate cells during liver fibrosis. J. Hepatol. 58, 736–742. doi: 10.1016/j.jhep.2012.11.022

PubMed Abstract | CrossRef Full Text | Google Scholar

Saraiva, N., Prole, D. L., Carrara, G., Johnson, B. F., Taylor, C. W., Parsons, M., et al. (2013). hGAAP promotes cell adhesion and migration via the stimulation of store-operated Ca²⁺ entry and calpain 2. J. Cell Biol. 202, 699–713. doi: 10.1083/jcb.201301016

PubMed Abstract | CrossRef Full Text | Google Scholar

Saxton, A. M., and SAS Institute (2004). Genetic Analysis of Complex Traits Using SAS. Cary, NC: SAS Institute.

Google Scholar

Sentandreu, M., Coulis, G., and Ouali, A. (2002). Role of muscle endopeptidases and their inhibitors in meat tenderness. Trends Food Sci. Technol. 13, 400–421. doi: 10.1016/S0924-2244(02)00188-7

CrossRef Full Text | Google Scholar

Shibata, M., Sato, T., Nukiwa, R., Ariga, T., and Hatakeyama, S. (2012). TRIM45 negatively regulates NF-κB-mediated transcription and suppresses cell proliferation. Biochem. Biophys. Res. Commun. 423, 104–109. doi: 10.1016/j.bbrc.2012.05.090

PubMed Abstract | CrossRef Full Text | Google Scholar

Sosnicki, A. A., and Newman, S. (2010). The support of meat value chains by genetic technologies. Meat Sci. 86, 129–137. doi: 10.1016/j.meatsci.2010.04.034

PubMed Abstract | CrossRef Full Text | Google Scholar

Steinhardt, R. A., Bi, G., and Alderton, J. M. (1994). Cell membrane resealing by a vesicular mechanism similar to neurotransmitter release. Science 263, 390–393. doi: 10.1126/science.7904084

PubMed Abstract | CrossRef Full Text | Google Scholar

Storey, J. D., and Tibshirani, R. (2003). Statistical significance for genomewide studies. Proc. Natl. Acad. Sci. U.S.A. 100, 9440–9445. doi: 10.1073/pnas.1530509100

PubMed Abstract | CrossRef Full Text | Google Scholar

Strydom, P., Naude, R., Smith, M., Scholtz, M., and Van Wyk, J. (2000). Characterisation of indigenous African cattle breeds in relation to meat quality traits. Meat Sci. 55, 79–88. doi: 10.1016/S0309-1740(99)00128-X

PubMed Abstract | CrossRef Full Text | Google Scholar

Südhof, T. C. (2002). Synaptotagmins: why so many? J. Biol. Chem. 277, 7629–7632. doi: 10.1074/jbc.R100052200

PubMed Abstract | CrossRef Full Text | Google Scholar

Tizioto, P. C., Decker, J. E., Taylor, J. F., Schnabel, R. D., Mudadu, M. A., Silva, F. L., et al. (2013a). Genome scan for meat quality traits in Nelore beef cattle. Physiol. Genomics 45, 1012–1020. doi: 10.1152/physiolgenomics.00066.2013

PubMed Abstract | CrossRef Full Text | Google Scholar

Tizioto, P. C., Gasparin, G., Souza, M. M., Mudadu, M. A., Coutinho, L. L., Mourão, G. B., et al. (2013b). Identification of KCNJ11 as a functional candidate gene for bovine meat tenderness. Physiol. Genomics 45, 1215–1221. doi: 10.1152/physiolgenomics.00137.2012

PubMed Abstract | CrossRef Full Text | Google Scholar

Trapnell, C., Pachter, L., and Salzberg, S. L. (2009). TopHat: discovering splice junctions with RNA-Seq. Bioinformatics 25, 1105–1111. doi: 10.1093/bioinformatics/btp120

PubMed Abstract | CrossRef Full Text | Google Scholar

Trapnell, C., Williams, B. A., Pertea, G., Mortazavi, A., Kwan, G., Van Baren, M. J., et al. (2010). Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation. Nat. Biotechnol. 28, 511–515. doi: 10.1038/nbt.1621

PubMed Abstract | CrossRef Full Text | Google Scholar

Tuckow, A. P., Jefferson, S. J., Kimball, S. R., and Jefferson, L. S. (2011). Simvastatin represses protein synthesis in the muscle-derived C2C12 cell line with a concomitant reduction in eukaryotic initiation factor 2B expression. Am. J. Physiol. Endocrinol. Metab. 300, E564–E570. doi: 10.1152/ajpendo.00383.2010

PubMed Abstract | CrossRef Full Text | Google Scholar

Uytterhaegen, L., Claeys, E., and Demeyer, D. (1994). Effects of exogenous protease effectors on beef tenderness development and myofibrillar degradation and solubility. J. Anim. Sci. 72, 1209–1223. doi: 10.2527/1994.7251209x

PubMed Abstract | CrossRef Full Text | Google Scholar

Uzureau, S., Coquerelle, C., Vermeiren, C., Uzureau, P., Acker, A., Pilotte, L., et al. (2016). Apolipoproteins L control cell death triggered by TLR3/TRIF signaling in dendritic cells. Eur. J. Immunol. 46, 1854–1866. doi: 10.1002/eji.201546252

PubMed Abstract | CrossRef Full Text | Google Scholar

Van Wynsberghe, P. M., Chan, S. P., Slack, F. J., and Pasquinelli, A. E. (2011). Analysis of MicroRNA Expression and Function. In Methods in Cell Biology, Vol. 106, Cambridge, MA: Academic Press, 219–252. doi: 10.1016/B978-0-12-544172-8.00008-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Vandenabeele, P., Orrenius, S., and Zhivotovsky, B. (2005). Serine Proteases and Calpains Fulfill Important Supporting Roles in the Apoptotic Tragedy of the Cellular Opera. London: Nature Publishing Group.

Google Scholar

Voet, D., Voet, J. G., and Pratt, C. W. (2013). Fundamentals of Biochemistry: Life at the Molecular Level, 4th Edn, Hoboken, NJ: Wiley.

Google Scholar

Wang, Y., Zhang, C., Fang, X., Zhao, Y., Chen, X., Sun, J., et al. (2014). Identification and profiling of microRNAs and their target genes from developing Caprine skeletal muscle. PLoS One 9:e96857. doi: 10.1371/journal.pone.0096857

PubMed Abstract | CrossRef Full Text | Google Scholar

Wheeler, T., Cundiff, L., Shackelford, S., and Koohmaraie, M. (2005). Characterization of biological types of cattle (Cycle VII): Carcass, yield, and longissimus palatability traits 1 2. J. Anim. Sci. 83, 196–207. doi: 10.2527/2005.831196x

PubMed Abstract | CrossRef Full Text | Google Scholar

Wheeler, T., Shackelford, S., and Koohmaraie, M. (1997). Standardizing collection and interpretation of Warner-Bratzler shear force and sensory tenderness data. Paper Presented at the Proceedings Of the Recip. Meat Conf. Clay Center, Washington, DC.

Google Scholar

Wing, S. S. (2013). Deubiquitinases in skeletal muscle atrophy. Int. J. Biochem. Cell Biol. 45, 2130–2135. doi: 10.1016/j.biocel.2013.05.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Wu, J., Lei, G., Mei, M., Tang, Y., and Li, H. (2010). A novel C53/LZAP-interacting protein regulates stability of C53/LZAP and DDRGK domain-containing Protein 1 (DDRGK1) and modulates NF-κB signaling. J. Biol. Chem. 285, 15126–15136. doi: 10.1074/jbc.M110.110619

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, X., Koltes, J. E., Park, C. A., Chen, D., and Reecy, J. M. (2015). Gene Co-expression network analysis provides novel insights into myostatin regulation at three different mouse developmental timepoints. PLoS One 10:e0117607. doi: 10.1371/journal.pone.0117607

PubMed Abstract | CrossRef Full Text | Google Scholar

Zeng, Y., Yi, R., and Cullen, B. R. (2005). Recognition and cleavage of primary microRNA precursors by the nuclear processing enzyme Drosha. EMBO J. 24, 138–148. doi: 10.1038/sj.emboj.7600491

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, Y., Zan, L., and Wang, H. (2011). Screening candidate genes related to tenderness trait in Qinchuan cattle by genome array. Mol. Biol. Rep. 38, 2007–2014. doi: 10.1007/s11033-010-0323-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhao, C., Carrillo, J. A., Tian, F., Zan, L., Updike, S. M., Zhao, K., et al. (2015). Genome-wide H3K4me3 analysis in Angus cattle with divergent tenderness. PLoS One 10:e0115358. doi: 10.1371/journal.pone.0115358

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhao, C., Tian, F., Yu, Y., Liu, G., Zan, L., Updike, M. S., et al. (2012a). miRNA-dysregulation associated with tenderness variation induced by acute stress in Angus cattle. J. Anim. Sci. Biotechnol. 3:12. doi: 10.1186/2049-1891-3-12

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhao, C., Tian, F., Yu, Y., Luo, J., Hu, Q., Bequette, B. J., et al. (2012b). Muscle transcriptomic analyses in Angus cattle with divergent tenderness. Mol. Biol. Rep. 39, 4185–4193. doi: 10.1007/s11033-011-1203-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhbannikov, I. Y., Hunter, S. S., and Settles, M. L. (2013). SeqyClean User Manual. New York, NY: Elsevier.

Google Scholar

Zimin, A. V., Delcher, A. L., Florea, L., Kelley, D. R., Schatz, M. C., Puiu, D., et al. (2009). A whole-genome assembly of the domestic cow. Bos taurus. Genome Biol. 10:R42. doi: 10.1186/gb-2009-10-4-r42

PubMed Abstract | CrossRef Full Text | Google Scholar