RNA sequencing and transcriptomes analysis
A total of 682,547,111 raw paired-end reads (204.77 Gb) were generated from RNA sequencing of beef cut samples. Specifically, an average of 22.75, 23.61, 23.26, 22.11, and 22.02 million reads were obtained from the tenderloin, longissimus lumborum, rump, neck, and chuck cuts, respectively (Additional file 1: Table S1). In total, 96.97% of the raw reads passed the quality control, and an average of 95.27% (ranging from 94.23 to 97.02%) clean reads were mapped to the bovine reference genome ARS-UCD1.2. To quantify the gene expression level, the FPKM values of genes were calculated based on the length of the gene and the read counts mapped to the gene, and 15,701 genes were identified using StringTie software [12]. In this study, we found the total number of expressed genes was significantly lower in longissimus lumborum (9352 ± 600) than those in the chuck (10,071 ± 8), neck (10,060 ± 24), rump (10,074 ± 4), and tenderloin cuts (10,071 ± 11) (P < 0.05 was considered as the significant level) (Fig. 1a, Additional file 1: Table S1). In addition, 12,086 genes were commonly expressed in five beef cuts, and 209, 168, 159, 203, and 163 genes were uniquely expressed in tenderloin, longissimus lumborum, rump, neck, and chuck cuts, respectively (Fig. 1b). To eliminate the influence of confounding factors at the experimental level, we retained the genes with FPKM values greater than one in the six biological replicate samples. A total of 4511 genes were obtained for the downstream analyses.
Region-specific expression patterns analysis
Based on the similar detection methods as described by a previous study [13], we identified a total of 80 RSGs in 4511 genes (ranging from 3045 to 4126) from five types of beef cuts (Fig. 2a, Additional file 2: Table S2). Among beef cuts, we detected the largest (n = 30) RSGs in the longissimus lumborum cut, followed by 19 RSGs in the neck cut (Fig. 2a). The functional annotation and pathway enrichment of RSGs indicated the meat quality characteristics of the beef cut. For instance, the chuck-related RSGs were significantly enriched in 2-oxocarboxylic acid metabolism (P value = 1.00e-02) and biosynthesis of amino acids (P value = 3.70e-02), the longissimus lumborum-related RSGs for ribosome (P value = 3.09e-03) and proteasome complex (P value = 4.43e-02), the tenderloin-related RSGs for calcium ion binding (P value = 3.32e-03) and actin cytoskeleton (P value = 1.02e-02) (Fig. 2b, Additional file 3: Table S3).
To identify master regulators of RSGs involved in biological processes, we performed TF analysis on five types of beef cuts. In the chuck cut, a total of four TFs were identified including ESRRA (NES = 8.317), RXRA (NES = 4.130), RCOR1 (NES = 3.422), and ESR1 (NES = 3.153) (Fig. 2c, Additional file 4: Table S4). In the tenderloin and neck cut, six TFs were identified, respectively, namely TEAD4, POLR2A, TAF1, JUN, FOSL, and MAFK, as well as SIN3A, BCL3, JUND, CTCF, TEAD4, and EP300 (Fig. 2d, e). The TFs identified in other beef cuts were shown in Additional file 5: Fig. S1. Meanwhile, we investigated the PPI network using our list of RSGs in different beef cuts, the sub-networks revealed that the major clusters consist of a large number of RSGs in beef cuts. We found a strikingly consistent pattern among RSGs. Severn RSGs associated with ribosome function showed a high connection in sub-network, while four out of nine muscle fiber structure genes displayed a similar pattern (Additional file 5: Fig. S2). In addition, we performed an RSG expression pattern profile and observed that a cluster of RSGs was highly expressed in the tenderloin cut. Meanwhile, the hierarchical cluster analysis found the tenderloin cut was separated from other beef cuts, and the RSGs in cluster 1 were highly expressed in the tenderloin cut, indicating the differences among them (Fig. 3). We also observed that RSGs in cluster 2 have high expression levels than other RSGs (Fig. 3).
Region-specific modules analysis
To explore the associations between gene expression and meat quality, we performed co-expression analysis on the expression levels of genes in five types of beef cuts. Using a WGCNA approach, we divided the filtered 4511 genes into 13 co-expressed gene modules. To obtain region-specific modules, we assessed the association between 13 modules and five types of beef cuts. Under the criteria of the correlation coefficient (r > 0.60) and P-value (P < 1.0e-2), we identified seven region-specific modules in our analysis (Fig. 4a, Additional file 6: Table S5). For instance, the MEpink module (r = 0.77, P-value = 2.00e-6) was positively correlated with longissimus lumborum, the MEblue module (r = − 0.98, P-value = 2.00e-19) was negatively correlated with longissimus lumborum. (Fig. 4a). The MEmagenta module was significantly correlated with chuck (r = 0.65, P-value = 2.00e-04) and neck (r = 0.57, P-value = 2.00e-3), respectively. The MEgreen module (r = 0.97, P-value = 3.00e-18) and MElightyellow module (r = − 0.63, P-value = 3.00e-4) were closely related to tenderloin (Fig. 4a).
Pathway enrichment analysis showed that genes included in the MEpink module were involved in fatty acid metabolism (bta01212) and oxidative phosphorylation (bta00190). The genes included in the MEblue module were mainly involved in valine, leucine, and isoleucine degradation (bta00280) and fatty acid elongation (bta00062) (Fig. 4b). The functional annotations of genes in region-specific modules indicated that the MEpink module was mainly involved in proteasome core complex (GO:0005839), ion channel binding (GO:0044325) and oxidation-reduction process (GO:0055114). The MEblue module was mainly involved in mRNA processing (GO:0006397), protein phosphorylation (GO:0006468), and fatty acid beta-oxidation (GO:0006635). The genes within the MEmagenta module were involved in positive regulation of protein catabolic process (GO:0045732), 4 iron, 4 sulfur cluster binding (GO:0051539) and thioredoxin peroxidase activity (GO:0008379) (Fig. 5). In addition, the MEgreen module was mainly involved in the transition between fast and slow fiber (GO:0014883), actin filament binding (GO:0051015) and lipid biosynthetic process (GO:0008610). The MEroyalblue module was related to the neck cut and mainly involved in protein autophosphorylation (GO:0046777), protein transport (GO:0015031) and fatty acid biosynthetic process (GO:0006633) (Fig. 5).
Candidate genes for meat quality among beef cuts
To detect the genes that affect the differences in meat quality among different beef cuts, we identified 91 candidate genes related to meat quality traits from RSGs and co-expression genes (Additional file 7: Table S6). These candidate genes were mainly involved in muscle fiber structure, fatty acids, amino acids, protein processing, energy production and conversion biological processes. For instance, 16 and 10 candidate genes were identified in muscle fiber structure and fatty acids, respectively. The gene expression patterns indicated that the expression levels of candidate genes were diverse among different beef cuts. Four genes (HSPB7, MYL12A, TNNT1, MYLK) were up-regulated in tenderloin cut compared with other beef cuts in terms of muscle fiber structure (Fig. 6, Additional file 7: Table S6). In fatty acids, we also observed that the ACADM gene was uniquely expressed in tenderloin cut. Three (ACADS, ALDH2, IVD) and six genes (CRIP2, CNBP, LIMS2, PHPT1, SNTA1, SUMO1) were involved in amino acids and ion channel binding, respectively (Additional file 5: Fig. S3, Additional file 7: Table S6). A total of 17 genes were involved in protein processing, of which seven genes belong to RSGs, (Additional file 5: Fig. S4, Additional file 7: Table S6). In addition, we found 39 genes (e.g., ATP8, COX8B, NDUFB6) were involved in energy production and conversion. These genes show significantly different expressions among beef cuts (Additional file 5: Fig. S5, Additional file 7: Table S6).
Validation of the candidate genes using RT-qPCR
To validate the expression accuracy of the candidate genes derived from RNA-seq, we performed RT-qPCR for five types of beef cuts from the same individual samples. Four candidate genes including ALDH2, CANX, IVD, PHPT1 were selected for the subsequent analysis. Consistent with gene expression changes from RNA-seq analysis, we found significant differences among the gene expression values among five types of beef cuts (Additional file 5: Fig. S6). There were significant differences in the expression of candidate genes between tenderloin and chuck, for instance, ALDH2 (P = 0.02996), CANX (P = 0.0375). IVD gene showed different expressions between tenderloin and rump. Meanwhile, significant differences were observed between longissimus lumborum and chuck for ALDH2 (P = 0.00516), PHPT1 (P = 5.48e-06). Our RT-qPCR analyses confirmed that candidate genes for meat quality were differentially expressed across beef cuts and further supported RNA-seq results.