Expression profile and bioinformatics analysis of circRNA and its associated ceRNA networks in longissimus dorsi from Lufeng cattle and Leiqiong cattle
BMC Genomics volume 24, Article number: 499 (2023)
This paper aims to explore the role of circRNA expression profiles and circRNA-associated ceRNA networks in the regulation of myogenesis in the longissimus dorsi of cattle breeds surviving under subtropical conditions in southern China by RNA sequencing and bioinformatics analysis. It also aims to provide comprehensive understanding of the differences in muscle fibers in subtropical cattle breeds and to expand the knowledge of the molecular networks that regulate myogenesis. With regard to meat quality indicators, results showed that the longissimus dorsi of LQC had lower pH (P < 0.0001), lower redness (P < 0.01), lower shear force (P < 0.05), and higher brightness (P < 0.05) than the longissimus dorsi of LFC. With regard to muscle fiber characteristics, the longissimus dorsi of LQC had a smaller diameter (P < 0.0001) and higher density of muscle fibers (P < 0.05). The analysis results show that the function of many circRNA-targeted mRNAs was related to myogenesis and metabolic regulation. Furthermore, in the analysis of the function of circRNA source genes, we hypothesized that btacirc_00497 and btacirc_034497 may regulate the function and type of myofibrils by affecting the expression of MYH6, MYH7, and NEB through competitive linear splicing.
According to the Domestic Animal Diversity Information System (FAO-DAD-IS), there are over 70 cattle breeds widely distributed in China, with 53 of them being indigenous [1, 2]. Native breeds are divided into three categories in accordance with their geographical distribution, including northern-distributed group, central-distributed group, and southern-distributed group, which are distributed in northern China, middle and lower reaches of the Yellow River, and southern China, respectively [3, 4]. Native Chinese cattle were once bred as draft animals, which has allowed them to adapt to the climate in which they grow and to be highly resistant to adversity . Natural selection and environmental characteristics have shaped the genomes and expression profiles of cattle grown in each environment. Cattle that have been grown in southern China for a long time are adapted to hot and humid climates, and they have excellent meat quality and high disease resistance [5, 6].
Meat quality is an important element in the assessment of the meat value of livestock . Beef myofibers can affect the traits to a certain extent, such as tenderness, color, and pH drop of beef . Thus, among the factors that can affect meat quality, myofiber is the strongest [9, 10]. However, the growth and development of muscle fibers include ontogeny at different embryonic stages, as well as hypertrophy and transformation at postnatal stages [9, 11]. In addition, the process of myogenesis is a complex biological process regulated by multiple factors, such as myogenic regulatory factors , signaling pathways , noncoding RNAs (ncRNAs) , and genes. Among these factors, circular RNA (circRNA) is an ncRNA with a covalently closed continuous loop structure that can play a role in gene regulation by competing with linear splicing , and this factor plays an important regulatory role in myogenesis . In addition, the function of circRNA as a sponge of miRNA may indirectly affect the translation of mRNA . CircRNA can be found in bovine muscle, and it is involved in the regulation of myogenesis [18,19,20].
Leiqiong cattle (LQN) and Lufeng cattle (LFN) are native cattle breeds in South China, which are primarily distributed in Guangdong and Hainan Provinces. These two kinds of cattle live in a subtropical climate environment for a long time, and they are less affected by artificial introduction. In addition, the breeding scale is considerable. The genetic relationship between LQN and LFN cattle is close, but their physical appearance, such as is quite different . In this study, we compared the circRNA transcripts in longissimus dorsi from two types of southern Chinese cattle, then constructed and analyzed the circRNA-associated ceRNA network. These results could provide comprehensive understanding of the differences in muscle fiber development in cattle grown under subtropical conditions and extend our understanding of the molecular networks that regulate beef meat quality.
Characteristics in meat quality and myofber
The result related to meat quality were shown in Table 1. Compared potential of hydrogen (pH) of longissimus dorsi from two varieties of cattle, pH value of LQC was significantly lower than LFC (P < 0.0001). As for the muscle color, the muscle lightness (L*) of LFC was significantly lower than that of LQC (P < 0.05), and the muscle redness (a*) of LFC was significantly higher than that of LQC (P < 0.01). It was no significantly differences between LFC and LQC on yellowness (b*) (P > 0.05). In the assessment of tenderness, tangential stress of LQC is lower than LFC (P < 0.05), that may preliminarily indicate that beef from LQC is more tender than that from LFN.
Compared the longissimus dorsi fibers from two kinds of cattle (Fig. 1), in the same level viewing area, myofibers quantity of LQC was significantly more than that of LFC (P < 0.05). For the area and diameter of longissimus dorsi fibers, LQC has smaller myofibers area and smaller diameter compared with LFC (P < 0.0001).
miRNA and transcriptome expression analysis
miRNA-seq generated 24,041,508 row reads from 8 samples, and after filtering, 23,312,472 clean reads were obtained, accounting for 97.00 ± 0.05% of the row reads. The reference genome, Bos taurus, gene data was downloaded in miRBase and analyzed using miRDeep2 software, 92.21 ± 1.88% of clean reads could be mapped to the reference genome, and a total of 493known miRNAs as well as 538 miRNA precursors were obtained. Predictive analysis of sequences not annotated to any information using mireap yielded 34.88 ± 8.22 new miRNAs (Table S1). Among all the annotated miRNAs, 59 and 51 miRNAs were unique to LQC and LFC, respectively, and 572 miRNAs were shared (Fig. 2(a)). Further analysis of miRNAs differentially expressed (DEseq, Expression difference fold: |log2FoldChange| > 1, Significance: P-value < 0.05) between LQC and LFC by the longissimus dorsi, with LQC as the control group, revealed 7 up-regulated and 6 down-regulated (Fig. 2(b)& Table S2).
For the RNA sequencing of mRNA and circRNA, a total of 8 samples from LQC and LFC yielded an average of 108,794,942 reads per sample. After quality control and trimming using Cutadapt, there are 94,857,498 per sample clean reads were yielded. TopHat2’s upgraded HISAT2 software (http://ccb.jhu.edu/software/hisat2/index.shtml) was used to map the reads, and 97.14%±0.28% of the clean reads can be mapped to the reference genome, Bos taurus (Table S1). A total of 17,313 transcripts were found in all samples detected, of which 16,810 transcripts were found in samples from LQC and 993 were unique; 16,320 were found in samples from LFC and 503 were unique; and 15,817 transcripts were shared by LQC and LFC (Fig. 2(d)). Differential analysis of gene expression was performed using DESeq (Expression difference fold: |log2FoldChange| > 1, Significance: P-value < 0.05), with LQC as the control group, 155 up-regulated and 444 down-regulated (Fig. 2(e)&Table S2). The PCA plot of miRNA and mRNA expression profile is shown in Fig. 2(c& f). GO (Gene Ontology) analysis revealed that these differentially expressed transcripts were significantly enriched (P < 0.05), and the vast majority (77.7%) of GO terms were Biological Process (Table S3).
After alignment with the reference genome, the unmatched Reads double-end 20 bp from the HISAT2 alignment results were intercepted as Anchors sequences and re-matched to the genome using Bowtie2 to detection of circRNA. 4,660,322 of these sequences were re-matched to the genome, representing 85.80 ± 3.50% of all Anchors sequences. Using find_circ to identify and annotate circRNAs, a total of 5,715 annotated entries were obtained (Table S4). Expression profile about circRNAs in LQC and LFC was showed in Fig. 3. The sourcegenes of identified circRNA vast majority were annotated exons (annot exons), and the spliced length of most identified circRNAs were concentrated about 500 bp (Fig. 3(a& b)). In the other hand, identified circRNAs were abundant on chromosome 1, 2, 3, 10 and 11. Compared with each other, there were 1,310 unique expressed circRNAs in LQC and 1,075 in LFC, 3,330 circRNAs were shared (Fig. 3(c& d)). PCA plot (Fig. 3(e)) showed that obviously separation exist between the samples of LQC and the samples of LFC, indicating that there was a certain difference in expression of circRNAs between the two types of cattle.
Furthermore, differential expressed circRNAs (DE circRNAs) between LQC and LFC were showed in Fig. 4& Table S2 (DEseq, Expression difference fold: |log2FoldChange| > 1, Significance: P-value < 0.05). Compared to LQC, 19 circRNAs were up regulated and 10 were down regulated in LFC (Fig. 4(a& c)), and accounts for 65.5% and 34.5% of all DE circRNAs (Fig. 4(b)), respectively.
Construction of ceRNA-networks and Enrichment analysis
After screening the DE circRNAs, psRobot and miRanda databases were used to predict targeted miRNA. As the co-expression networks in Fig. 5 and Table S5 shown, 13 miRNA and 135 mRNAs were associated with 8 circRNAs. Above them, there were 4 up-regulated and 4 down-regulated circRNAs, and influenced by them, 7 miRNAs and 35 mRNAs were up-regulated, and 6 miRNAs and 100 mRNAs were down-regulated.
GO and KEGG (Kyoto Encyclopedia of Genes and Genomes) enrichment analysis were performed on mRNAs that were indirectly influenced by circRNAs. As shown in Fig. 6(a& b) and Table S6, the top20 GO terms that targeted genes enriched in were mostly associated with Cellular Component and Biological Process. Further refinement, the top20 GO terms were closely related to peptidyl-tyrosine phosphorylation, peptidyl-tyrosine modification, regulation of peptidyl-tyrosine phosphorylation, plasma membrane part, myelin sheath, intrinsic component of plasma membrane, external side of plasma membrane, integral component of plasma membrane, regulation of cytosolic calcium ion concentration, cellular calcium ion homeostasis, myelination, ensheathment of neurons, axon ensheathment, calcium ion homeostasis, regulation of biological quality, positive regulation of cytosolic calcium ion concentration, cellular divalent inorganic cation homeostasis, regulation of midbrain dopaminergic neuron differentiation, regulation of planar cell polarity pathway involved in axis elongation, negative regulation of planar cell polarity pathway involved in axis elongation and so on. And the KEGG enrichment analysis results shown by Fig. 6(c& d) and Table S7, top20 KEGG Orthologys distributed in 5 categories, included of Pyrimidine metabolism, Nicotinate and nicotinamide metabolism, Rap1 signaling pathway, Platelet activation, MAPK signaling pathway, Purine metabolism, Complement and coagulation cascades, Calcium signaling pathway, Thyroid hormone signaling pathway, Focal adhesion, Glycine, serine and threonine metabolism, Toxoplasmosis, Vitamin B6 metabolism, Ras signaling pathway, PI3K-Akt signaling pathway, Regulation of actin cytoskeleton, Vascular smooth muscle contraction, Central carbon metabolism in cancer, Human cytomegalovirus infection and so on.
Enrichment analysis of sourcegene of circRNAs
CircRNA is usually produced by exons or introns of their sourcegenes, and the cyclization of circRNA probably affects the expression of the sourcegene . We analyzed sourcegenes of 29 DE circRNAs from LQC and LFC, then KEGG and GO enrichment analysis were performed on them (Fig. 7). The result of top 25 GO terms in Fig. 7(a) illustrate that most terms were enriched in Cellular Component and 6 terms were enriched in Biological Progress. The information of all terms that were enriched have been stored in Table S8, and it is worth noting that most of the enriched genes were associated with the function of myofibers, such as sarcomere, contractile fiber part, myofibril, contractile fiber, stress fiber, myosin complex, contractile actin filament bundle, actin filament bundle, actomyosin, muscle filament sliding, actin-myosin filament sliding and so on. Figure 7(b) and Table S9 shown that top 25 KEGG Orthology of related genes, among them, there were 1, 1, 2, 3, 5 and 13 KEGG Orthology were classified to Cellular Processes, Environmental Information Processing, Genetic Information Processing, Human Diseases, Organismal Systems and Metabolism, respectively.
In order to verify the reliability of the RNA sequencing results, 3 circRNAs, 3 miRNAs and 8 mRNAs that differentially expressed in LQC and LFC were randomly selected and analyzed by RT-qPCR. The LQC served as a control group, RNA-seq and RT-qPCR results were consistent in the expression levels of each selected genes (Figure S1). Information about the primers used to amplify the selected genes has been placed in the Table S10.
Differences in meat quality between longissimus dorsi from LQC and LFC may be related to the differences in myofiber composition
LFC and LQC are native Chinese cattle breeds, which have adapted to the hot and humid climate of southern China after long evolution; therefore, it is important to enhance their meat value potential. Most of the indicators used for meat quality assessment, such as pH decline, meat color, and tenderness, are related to the structure and metabolism of the muscle fibers [22, 23]. The rate of pH decline is often related to the amount of glycogen reserves in muscle tissue before slaughter and the content of mitochondria in different muscle fibers. When the glycogen level in muscle tissue is low, the rate of glycolysis decreases; therefore, the pH declines gradually because of the slow accumulation of lactic acid . Muscles with more oxidative properties had a higher pH than those with more glycolytic properties after the animals were slaughtered , and the meat with a higher percentage of fast-twitch glycolytic fibers has a higher rate and degree of pH decline, it is commonly known in the field of red muscle. Furthermore, the rate of glycolysis and rate of pH decrease can affect changes of some proteins in muscle fibers at the post-mortem, such as tropomyosin, actin, troponin, and some enzymes on glycolysis, and these protein changes will be directly reflected in meat tenderness, flavor substances, color, and other measures of meat quality . Meat color is an important point in determining meat quality, and it is a marketability criterion of meat because consumers often use this criterion to select and buy meat. Previous studies [26, 27] have shown that meat with low glycogen reserves tends to have higher ultimate pH, which leads to low light scattering and high oxygen consumption in meat surface, resulting in low lightness of the meat. The hue and chromaticity of meat color are primarily dominated by myoglobin (Mb) because the changes in the biochemical state of Mb cause the changes in meat color, particularly the degree of oxidation and reduction of Mb. In addition, Mb in fresh meat is present in four states, namely, deoxymyoglobin, oxymyoglobin (OxyMb), carboxymyoglobin (COMb), and metmyoglobin, in which OxyMb and COMb provide the meat with a bright cherry red color that is typical of fresh meat . However, the four redox states of Mb are not constant, and many meat-endogenous factors can affect the color of meat by influencing the state of Mb, of which pH, muscle source, presence of antioxidants, lipid oxidation, and mitochondrial activity and the most prominent . Beef tenderness has received considerable attention because it affects the return decision and satisfaction of consumers. The characteristics of the muscle fibers largely affect the tenderness of the meat, although the tenderness and texture of beef are also affected by the connective tissue and intermuscular fat in the muscle . The activity of some protein hydrolytic enzymes changes with the decrease of the rate of pH, which leads to differences in glycogen metabolism by muscle fibers that indirectly affect meat tenderness [31, 32]. In addition, the number of muscle fibers affect the tenderness of meat. In general, tender meat has more muscle fibers per unit area and a smaller diameter . In the present study, the pH of the longissimus dorsi of LQC decreased faster than that of LFC, which may be one of the reasons for the high lightness of the longissimus dorsi of LQC. With regard to redness, the longissimus dorsi of LFC has a higher level of redness than that of LQC, which may be due to some factors that keep the Mb more in the state of OxyMb or COMb in LFC. By analyzing the quantitative characteristics of myofibers per unit area, we found that the number of myofibers in LQC was less than that in LFC, and the diameter of myofibers was smaller. Such quantity characteristics of myofibers in LQC also confer a smaller shear force in LQC.
ceRNA network points to the regulation of PI3K-Akt, MAPK, and calcium pathways
Considering that muscle fibers profoundly influence meat quality, the regulatory mechanism of myogenesis has been discussed comprehensively and extensively. Myogenesis is a physiological process that is influenced by various internal or external factors starting from fetal life. Among the factors, circRNA as a ncRNA has received considerable attention for its wide range of regulatory functions. circRNAs are dynamically expressed and abundant in muscle tissue of many species, including humans, cattle, goat, pig, chicken, and mice . Although the functions of circRNAs remain largely unexplored, they can serve as miRNA sponges and further contribute to mRNA stability or protein production . For example, circPTPN4 can sponge miR-499-3p, which regulates NAMPT expression, thereby promoting myoblast proliferation and differentiation and activating the fast-twitch muscle phenotype . The overexpression of circCPE counteracts the inhibitory effect of miR-138 on cell proliferation and the accelerating effect on differentiation and apoptosis, and circRNAs can reduce the inhibitory effect of miR-138 on FOSC1, which is involved in myogenesis . With the continuous improvement of sequencing technology and database, the construction of a ceRNA network has become an important tool for predicting circRNA function. Therefore, we performed Pearson correlation analysis to assess the association between DE circRNA and DE mRNA that exist between LQC and LFC, and enrichment analysis was performed for target mRNAs. Based on the results of KEGG analysis, we found that target mRNAs were primarily related to PI3K-Akt, MAPK, and calcium signaling pathways. PI3K-Akt has been widely proven to be associated with muscle hypertrophy and atrophy [37,38,39,40]. Several recent studies have also demonstrated the ability of circRNA to regulate the PI3K/AKT pathway, for example, the newly identified circRILPL1 as miR-145 sponge promotes myogenic cell growth by regulating the expression of the IGF1R gene to reduce the inhibitory effect of miR-145 on the PI3K/AKT signaling pathway . The proliferation and differentiation of bovine myoblasts can be inhibited by circMEF2D, and circMEF2D can regulate the PI3K-AKT signaling pathway by competitively binding miR-486 . Previous studies suggest that MAPK may be associated with myoblast cell cycle arrest, which is critical for initiating muscle differentiation in myogenic cells . In addition, MAPK is a signaling pathway that drives the metabolic adaptation of skeletal muscle to exercise . Numerous studies have shown that MAPK increases insulin-dependent glucose uptake and oxidative metabolism as well as mitochondrial oxidative phosphorylation in muscle during exercise [44,45,46,47]. In muscle tissue, intracellular stores of Ca2+, upon release, trigger the formation of actin crossbeams and the generation of contractile force . Myocytes can sense and respond to changes in workload and activation patterns by regulating the gene expression and cellular metabolism of calcium signaling pathways [49, 50]. Stored calcium influx has become a mechanism by which the calcium signaling pathway is activated in response to the changing demands of myocytes [48, 51]. Most of the target mRNAs identified by the ceRNA network using LQC and LFC longissimus dorsi were enriched in pathways related to muscle metabolism and myogenesis. Therefore, these mRNA-related circRNAs and miRNAs might be effective tools for regulating muscle production and development.
Cyclization of btacirc_00497 and btacirc_034497 may affect myogenesis
Apart from being a sponge for miRNA, circRNA can compete with linear splicing to play a functional role in gene regulation . In this study, we attempted to use GO and KEGG enrichment analyses to organize and predict the function of genes that were affected by and were the source of DE circRNA. Based on the results of enrichment analysis, MYH6, MYH7, and nebulin (NEB) have high frequency in multiple pathways and GO terms that may associate with myogenesis and meat quality. Myosin is a motor protein that plays an important role in the contraction of animal skeletal muscle. Myosin consists of six subunits, two of which are myosin heavy chain (MyHC) subunits. MyHC isoforms play an irreplaceable role in muscle contraction because they have ATPase activity, which provides energy for muscle contraction . In mammals, 11 known genes can encode MyHC, and these genes are highly conserved during evolution . The four MyHC isoforms, namely, I, IIa, IIx, and IIb, are classified by different coding genes, which leads to four myofiber types. MyHC I is expressed in type I fibers, and IIa, IIx, and IIb are expressed in IIA, IIX, and IIB fiber types, respectively [53, 54]. The two subunits of the MyHC I heavy chain are composed of β-cardiac and α-cardiac encoded by MYH7 and MYH6, respectively . Theoretically, the high expression level of MyHC I indicates that muscle fiber metabolism will have less efficient glycolysis, resulting in less lactate accumulation . Low levels of lactic acid in the muscle tissue of live cattle can lead to a gradual decrease in beef pH rate postmortem . Insufficient pH reduction in beef impairs meat color, tenderness, and shelf life . Many studies have shown that feeding high-energy diet before slaughter is an effective way to increase muscle glycogen content and improve muscle pH reduction after slaughter, and this strategy may not become dependent on muscle fibers that have high MyHC I expression level [58, 59]. In our study, btacirc_00497 in LFC longissimus dorsi has a high expression level, which is sourced from MYH6 and MYH7. However, our findings on the meat quality of LQC and LFC do not support our conjecture, that is, circRNA reduces β-MHC and α-MHC expression by competing for linear shearing of the transcripts of MYH6 and MYH7 and further influences the metabolic process of MyHC. We hypothesize that the expression level of MYH6 and MYH7 transcripts of LFC may be higher than that of LQC, so the probability of cyclization to btacirc_00497 is also increased. Thus, although the factors affecting the cyclization of btacirc_00497 are unclear, this gene may still serve as a potential lynchpin for the regulation of muscle fiber types. In addition, we identified btacirc_034497, which shares the same origin gene as NEB. Actin is present and functionally important in most eukaryotic cells, and the control of actin filament organization and structure is critical for many cellular functions. NEB stabilizes actin filaments in thin filament architecture, thereby regulating the filament length [60, 61]. Moreover, NEB regulates skeletal muscle contraction in skeletal muscle, and muscles from NEB knockout mice produce significantly less force than their wild-type siblings [62, 63]. Therefore, the factors affecting btacirc_034497 cyclization have great application potential as important tools for NEB regulation.
Materials and methods
All experimental animals were sourced from the breeding farm (Meixian Country, Meizhou City, Guangdong province). Four Leiqiong cattle (LQC) and four Lufeng cattle (LFC) each, four-month-old, were humanely slaughtered, and the cattle of the same breed are similar in weight and body condition and raised in the same farm environment. Then longissimus dorsi tissues from 8 cattle were collected and snap-frozen in liquid nitrogen immediately to extract total RNA. Another portion of the longissimus dorsi tissue was trimmed into 0.5×0.5×1.0 cm pieces and immediately fixed in 4% paraformaldehyde for observation in tissue sections. In addition, any anesthesia or euthanizing agent was not used in our study.
Analysis of muscle quality properties
After 24 hours of slaughter, pH was determined using Meat pH Meter HI99163 (HANNA, Italy) on the longissimus dorsi. The colorimetric parameters of the muscles were calculated using the L* a * b * system with a colorimeter OPTO-STAR Meat Color Tester (Matthaus, Germany) from the average of three random readings of each sample. And the evaluation of steaming losses was performed according to the project of Honikel . During the shear force measurement, three samples of the longissimus dorsi after 72 hours of aging of each cattle (1cm×1cm×3cm) were measured three times perpendicular to the fiber direction. The final shear force was the average of the three readings, and the measurements were expressed in Pascal(Pa) .
Muscle fiber characteristics
The tissues were after fixing by 4% paraformaldehyde for 24h, and then immersed in xylene alcohol (1:1, v/v), infiltrated and embedded in paraffin. Cross sections of 3µm thickness were prepared, stained with hematoxylin and eosin, observed under a microscope, and photographed (400×magnification). The number of muscle fibers and total cross-sectional area were subsequently evaluated with Image-Pro Plus V6.0 (Media Cybernetics Inc., Rockville, MD, USA).
Preparation and sequencing for RNA-seq
TRIzol (Thermo Fisher, Shanghai, China) used to extract total RNA from the tissue samples. RNA quantity and purity were determined using Agilent 2100 Bioanalyzer and the RNA 6000 Nano Labchip Kit (Agilent, Santa Clara, USA). Epicentre Ribo-Zero™rRNA Removal Kit (Epicentre, America) were used to remove rRNA before the RNA-seq library construction. Then, the rRNA-depleted RNAs were then fragmented and reverse-transcribed to obtain cDNA libraries using the TruseqTM RNA sample prep Kit (Illumina, America). Small RNA-seq library construction was performed according to the instructions of NEBNext Multiplex Small RNA library Prep (Illumina, America), which was followed by TRIzol extraction. Then, the library was amplified and enriched using PCR technique, proprietary indexed adapters were then ligated to 5ʹ- and 3ʹ-termini, and later electrophoresed using 15% concentration of agarose gel to obtain the target fragments. Finally, all sequenced libraries were sent to the sequencing company (Personalbio, Shanghai, China) and sequenced by Illumina Hiseq4000 platform for paired-end sequencing after quality control.
Bioinformatics Primary Analysis
The reference genome used in this experiment was Bos taurus ftp:(ftp.ensembl.org/pub/release-101/fasta/bos_taurus/dna). For miRNA, Hiseq Single-End mode sequencing data were mapped between reference genome sequences after quality control using miRDeep2  software. The mature miRNA and precursor sequences of the species miRNAs were downloaded from miRbase (https://www.mirbase.org/). Then the de-duplicated sequences were mapped and annotate. Sequences that were not annotated with any information were analyzed using mireap (http://sourceforge.net/projects/mireap). DESeq V1.18.0  was performed to analyze miRNAs for differential expression and filtered by |log2FoldChange|>1 and P-value < 0.05 to filter out differentially conserved miRNAs.
For the circRNA and mRNA, downstream data were used for initial evaluation of row data using FastQC v0.11.9 software after RNA-seq was completed. Cutadapt v2.6 was used to filter lower quality data as well as adepters. The obtained clean data were used for the subsequent analysis of circRNA and mRNA. During the analysis of mRNA, the upgraded HISAT2 (http://ccb.jhu.edu/software/hisat2/index.shtml) software of TopHat2 was used to map reads with reference genes. The differentially expressed transcripts were analyzed using DEseq, and the differentially expressed genes were screened for |log2FoldChange| > 1, P-value < 0.05. During the circRNA prediction analysis, the 20 bp of two ends of the reads that were not matched on the HISAT2 were realignment and used as Anchors. Ffind_circ  was used for circRNA identification. CircRNA expression was analyzed for differential expression, and the screened standard was |log2FoldChange| > 1 and P-value < 0.05. GO (Gene Ontology) and KEGG (Kyoto Encyclopedia of Genes and Genomes) enrichment analysis [69,70,71] were performed using DAVID (https://david.ncifcrf.gov/). In addition, miRNA targeting relationships to circRNA and mRNA were predicted using miRanda and psRobot software.
After obtained sequencing results, RNA was extracted from the longissimus dorsi tissue using TRIzol (Thermo Fisher, Shanghai, China). RNA was reverse-transcribed to cDNA using the method provided in the introductions of PrimeScript RT Reagent Kit-Prefect Real Time (Takara, Beijing, China) for validation of mRNA and circRNA. miRNA 1st Strand cDNA Synthesis Kit (by stem-loop) (Vazyme, Nanjing, China) was used for miRNA reverse transcription. 2×Ultra SYBR Green qPCR Mix (CISTRO, Shanghai, China) was used for RT-qPCR in an Applied Biosystems QuantStudio 5 (Thermo Fisher, Shanghai, China). The reaction mixture containing, reaction conditions and primer sequences for the RT-qPCR validation procedure are stored in Table S10. The results were statistically analyzed using the 2-ΔΔCT method . All data are expressed as the average of 4 independent experiments.
The comparative analysis of two groups was performed using unpaired independent t-test. SPSS 25.0 (SPSS Inc., Chicago, IL) and GraphPad prism 9 was used for statistical analyses. Results were considered statistically significant differences at P < 0.05 and were expressed as the mean ± Standard Deviation unless otherwise stated.
In this study, we found that the pH drop and brightness of beef from LQC were significantly higher than those from LFC. In addition, the small cross-sectional area and diameter of muscle fibers produced better tenderness for beef from LQC. However, LFC was found to have superior redness compared with LQC. Furthermore, we identified several circRNAs related to muscle production and metabolism by ceRNA network building and analysis. By analyzing the source genes of DE circRNAs, we found that btacirc_00497 and btacirc_034497 may regulate the expression of MYH6, MYH7, and NEB by competing for linear shear, thereby altering the muscle fiber structure. Therefore, the ceRNA network-related genes and circRNA-derived genes identified in this study could be used as tools for regulating muscle production and meat quality of cattle grown in a subtropical environment.
The raw sequence data reported in this paper have been deposited in the Genome Sequence Archive (Genomics, Proteomics & Bioinformatics 2021) in National Genomics Data Center (Nucleic Acids Res 2022), China National Center for Bioinformation / Beijing Institute of Genomics, Chinese Academy of Sciences (CRA007241 and CRA007242) that are publicly accessible at https://ngdc.cncb.ac.cn/gsa.
Zhang W, Gao X, Zhang Y, Zhao Y, Zhang J, Jia Y, et al. Genome-wide assessment of genetic diversity and population structure insights into admixture and introgression in chinese indigenous cattle. BMC Genet. 2018;19(1):114. https://doi.org/10.1186/s12863-018-0705-9.
Huai Q, Zhiyong J, Zhijie C. A survey of cattle production in China. World Anim Rev. 1993, 76.
Zeng L, Chen N, Ning Q, Yao Y, Chen H, Dang R, et al. PRLH and SOD1 gene variations associated with heat tolerance in chinese cattle. Anim Genet. 2018;49(5):447–51. https://doi.org/10.1111/age.12702.
Hong C, Huai Q. Studies on sex chromosome polymorphism of four local cattle (Bos taurus) breeds in China. Hereditas. 1993.
Liu Y, Xu L, Yang L, Zhao G, Li J, Liu D, et al. Discovery of genomic characteristics and selection signatures in Southern Chinese local cattle. Front Genet. 2020;11:533052. https://doi.org/10.3389/fgene.2020.533052.
Lu X, Arbab AAI, Zhang Z, Fan Y, Han Z, Gao Q, et al. Comparative transcriptomic analysis of the Pituitary gland between cattle breeds differing in growth: yunling cattle and Leiqiong cattle. Anim (Basel). 2020;10(8). https://doi.org/10.3390/ani10081271.
Santiago GG, Siqueira F, Cardoso FF, Regitano LCA, Ventura R, Sollero BP, et al. Genomewide association study for production and meat quality traits in Canchim beef cattle. J Anim Sci. 2017;95(8):3381–90. https://doi.org/10.2527/jas.2017.1570.
Picard B, Gagaoua M. Muscle Fiber Properties in cattle and their Relationships with meat qualities: an overview. J Agric Food Chem. 2020;68(22):6021–39. https://doi.org/10.1021/acs.jafc.0c02086.
Zhuang X, Lin Z, Xie F, Luo J, Chen T, Xi Q, et al. Identification of circRNA-associated ceRNA networks using longissimus thoracis of pigs of different breeds and growth stages. BMC Genomics. 2022;23(1):294. https://doi.org/10.1186/s12864-022-08515-7.
Lee SH, Joo ST, Ryu YC. Skeletal muscle fiber type and myofibrillar proteins in relation to meat quality. Meat Sci. 2010;86(1):166–70. https://doi.org/10.1016/j.meatsci.2010.04.040.
Lefaucheur L. Myofiber typing and pig meat production. Slovenski Veterinarski Zbornik. 2001:5–29.
He S, Fu T, Yu Y, Liang Q, Li L, Liu J, et al. IRE1α regulates skeletal muscle regeneration through myostatin mRNA decay. J Clin Invest. 2021;e143737. https://doi.org/10.1172/JCI143737.
Coudert L, Osseni A, Gangloff YG, Schaeffer L, Leblanc P. The ESCRT-0 subcomplex component Hrs/Hgs is a master regulator of myogenesis via modulation of signaling and degradation pathways. BMC Biol. 2021;19(1):153. https://doi.org/10.1186/s12915-021-01091-4.
Archacka K, Ciemerych MA, Florkowska A, Romanczuk K. Non-coding RNAs as regulators of myogenesis and postexercise muscle regeneration. Int J Mol Sci. 2021;22(21). https://doi.org/10.3390/ijms222111568.
Ashwal-Fluss R, Meyer M, Pamudurti NR, Ivanov A, Bartok O, Hanan M, et al. circRNA biogenesis competes with pre-mRNA splicing. Mol Cell. 2014;56(1):55–66. https://doi.org/10.1016/j.molcel.2014.08.019.
Zhang P, Chao Z, Zhang R, Ding R, Wang Y, Wu W, et al. Circular RNA regulation of myogenesis. Cells. 2019;8(8). https://doi.org/10.3390/cells8080885.
Yue B, Wang J, Song C, Wu J, Cao X, Huang Y, et al. Biogenesis and ceRNA role of circular RNAs in skeletal muscle myogenesis. Int J Biochem Cell Biol. 2019;117:105621. https://doi.org/10.1016/j.biocel.2019.105621.
Liu R, Liu X, Bai X, Xiao C, Dong Y. Identification and characterization of circRNA in Longissimus Dorsi of different breeds of cattle. Front Genet. 2020;11:565085. https://doi.org/10.3389/fgene.2020.565085.
Huang K, Chen M, Zhong D, Luo X, Feng T, Song M, et al. Circular RNA profiling reveals an abundant circEch1 that promotes myogenesis and differentiation of bovine skeletal muscle. J Agric Food Chem. 2021;69(1):592–601. https://doi.org/10.1021/acs.jafc.0c06400.
Zhang R-M, Pan Y, Zou C-X, An Q, Cheng J-R, Li P-J, et al. CircUBE2Q2 promotes differentiation of cattle muscle stem cells and is a potential regulatory molecule of skeletal muscle development. BMC Genomics. 2022;23(1):267. https://doi.org/10.1186/s12864-022-08518-4.
Zhang Y, Yang Y, Ju H, He X, Sun P, Tian Y, et al. Comprehensive profile of circRNAs in formaldehyde induced heart development. Food Chem Toxicol. 2022;162:112899. https://doi.org/10.1016/j.fct.2022.112899.
Lefaucheur L. A second look into fibre typing–relation to meat quality. Meat Sci. 2010;84(2):257–70. https://doi.org/10.1016/jmeatsci200905004.
Klont RE, Brocks L, Eikelenboom G. Muscle fibre type and meat quality. Meat Sci. 1998;49S1:219–S229.
Wicks J, Beline M, Gomez JFM, Luzardo S, Silva SL, Gerrard D. Muscle energy metabolism, growth, and meat quality in beef cattle. Agriculture-Basel. 2019;9(9). https://doi.org/10.3390/agriculture9090195.
Whipple G, Koohmaraie M, Dikeman ME, Crouse JD. Predicting beef-longissimus tenderness from various biochemical and histological muscle traits. J Anim Sci. 1990;68(12):4193–9. https://doi.org/10.2527/1990.68124193x.
Hughes JM, Clarke FM, Purslow PP, Warner RD. Meat color is determined not only by chromatic heme pigments but also by the physical structure and achromatic light scattering properties of the muscle. Compr Rev Food Sci Food Saf. 2020;19(1):44–63. https://doi.org/10.1111/1541-4337.12509.
Purslow PP, Warner RD, Clarke FM, Hughes JM. Variations in meat colour due to factors other than myoglobin chemistry; a synthesis of recent findings (invited review). Meat Sci. 2020;159:107941. https://doi.org/10.1016/j.meatsci.2019.107941.
Suman SP, Joseph P. Myoglobin chemistry and meat color. Annu Rev Food Sci Technol. 2013;4:79–99. https://doi.org/10.1146/annurev-food-030212-182623.
Mancini RA, Hunt MC. Current research in meat color. Meat Sci. 2005;71(1):100–21. https://doi.org/10.1016/jmeatsci200503003.
Listrat A, Gagaoua M, Normand J, Gruffat D, Andueza D, Mairesse G, et al. Contribution of connective tissue components, muscle fibres and marbling to beef tenderness variability in longissimus thoracis, rectus abdominis, semimembranosus and semitendinosus muscles. J Sci Food Agric. 2020;100(6):2502–11. https://doi.org/10.1002/jsfa.10275.
White A, O’Sullivan A, Troy DJ, O’Neill EE. Manipulation of the pre-rigor glycolytic behaviour of bovine M. longissimus dorsi in order to identify causes of inconsistencies in tenderness. Meat Sci. 2006;73(1):151–6. https://doi.org/10.1016/j.meatsci.2005.11.021.
Gagaoua M, Picard B, Monteils V. Assessment of cattle inter-individual cluster variability: the potential of continuum data from the farm-to-fork for ultimate beef tenderness management. J Sci Food Agric. 2019;99(8):4129–41. https://doi.org/10.1002/jsfa.9643.
Albrecht E, Teuscher F, Ender K, Wegner J. Growth- and breed-related changes of muscle bundle structure in cattle. J Anim Sci. 2006;84(11):2959–64. https://doi.org/10.2527/jas.2006-345.
Li X, Yang L, Chen L-L. The Biogenesis, Functions, and Challenges of Circular RNAs. Mol Cell. 2018;71(3):428–42. https://doi.org/10.1016/j.molcel.2018.06.034.
Cai B, Ma M, Zhou Z, Kong S, Zhang J, Zhang X, et al. circPTPN4 regulates myogenesis via the miR-499-3p/NAMPT axis. J Anim Sci Biotechnol. 2022;13(1):2. https://doi.org/10.1186/s40104-021-00664-1.
Ru W, Qi A, Shen X, Yue B, Zhang X, Wang J, et al. The circular RNA circCPE regulates myoblast development by sponging miR-138. J Anim Sci Biotechnol. 2021;12(1):102. https://doi.org/10.1186/s40104-021-00618-7.
Glass DJ. Signalling pathways that mediate skeletal muscle hypertrophy and atrophy. Nat Cell Biol. 2003;5(2):87–90. https://doi.org/10.1038/ncb0203-87.
Yu M, Wang H, Xu Y, Yu D, Li D, Liu X, et al. Insulin-like growth factor-1 (IGF-1) promotes myoblast proliferation and skeletal muscle growth of embryonic chickens via the PI3K/Akt signalling pathway. Cell Biol Int. 2015;39(8):910–22. https://doi.org/10.1002/cbin.10466.
Stitt TN, Drujan D, Clarke BA, Panaro F, Timofeyva Y, Kline WO, et al. The IGF-1/PI3K/Akt pathway prevents expression of muscle atrophy-induced ubiquitin ligases by inhibiting FOXO transcription factors. Mol Cell. 2004;14(3):395–403. https://doi.org/10.1016/s1097-2765(04)00211-4.
Kumar A, Xie L, Ta CM, Hinton AO, Gunasekar SK, Minerath RA, et al. SWELL1 regulates skeletal muscle cell size, intracellular signaling, adiposity and glucose metabolism. Elife. 2020;9:e58941. https://doi.org/10.7554/eLife.58941.
Shen X, Tang J, Jiang R, Wang X, Yang Z, Huang Y, et al. CircRILPL1 promotes muscle proliferation and differentiation via binding miR-145 to activate IGF1R/PI3K/AKT pathway. Cell Death Dis. 2021;12(2):142. https://doi.org/10.1038/s41419-021-03419-y.
Zhang XY, Yang SL, Kang ZH, Ru WX, Shen XM, Li M et al. circMEF2D negatively regulated by HNRNPA1 inhibits proliferation and differentiation of Myoblasts via miR-486-PI3K/AKT Axis. J Agric Food Chem8145–63https://doi.org/10.1021/acs.jafc.2c01888.
Lee J, Hong F, Kwon S, Kim SS, Kim DO, Kang HS, et al. Activation of p38 MAPK induces cell cycle arrest via inhibition of Raf/ERK pathway during muscle differentiation. Biochem Biophys Res Commun. 2002;298(5):765–71. https://doi.org/10.1016/s0006-291x(02)02562-7.
Bengal E, Aviram S, Hayek T. p38 MAPK in glucose metabolism of skeletal muscle: beneficial or harmful? Int J Mol Sci. 2020;21(18). https://doi.org/10.3390/ijms21186480.
Somwar R, Koterski S, Sweeney G, Sciotti R, Djuric S, Berg C, et al. A dominant-negative p38 MAPK mutant and novel selective inhibitors of p38 MAPK reduce insulin-stimulated glucose uptake in 3T3-L1 adipocytes without affecting GLUT4 translocation. J Biol Chem. 2002;277(52):50386–95. https://doi.org/10.1074/jbc.M205277200.
Akimoto T, Pohnert SC, Li P, Zhang M, Gumbs C, Rosenberg PB, et al. Exercise stimulates Pgc-1alpha transcription in skeletal muscle through activation of the p38 MAPK pathway. J Biol Chem. 2005;280(20):19587–93. https://doi.org/10.1074/jbc.m408862200.
Pogozelski AR, Geng T, Li P, Yin X, Lira VA, Zhang M, et al. p38gamma mitogen-activated protein kinase is a key regulator in skeletal muscle metabolic adaptation in mice. PLoS ONE. 2009;4(11):e7934. https://doi.org/10.1371/journal.pone.0007934.
Kuo IY, Ehrlich BE. Signaling in muscle contraction. Cold Spring Harb Perspect Biol. 2015;7(2):a006023. https://doi.org/10.1101/cshperspect.a006023.
Tarasova NV, Vishnyakova PA, Logashina YA, Elchaninov AV. Mitochondrial calcium Uniporter structure and function in different types of muscle tissues in Health and Disease. Int J Mol Sci. 2019;20(19):4823. https://doi.org/10.3390/ijms20194823.
Stiber JA, Rosenberg PB. The role of store-operated calcium influx in skeletal muscle signaling. Cell Calcium. 2011;49(5):341–9. https://doi.org/10.1016/j.ceca.2010.11.012.
Berridge MJ. The Inositol Trisphosphate/Calcium Signaling Pathway in Health and Disease. Physiol Rev. 2016;96(4):1261–96. https://doi.org/10.1152/physrev.00006.2016.
Schiaffino S, Reggiani C. Fiber types in mammalian skeletal muscles. Physiol Rev. 2011;91(4):1447–531. https://doi.org/10.1152/physrev.00031.2010.
Bottinelli R, Reggiani C. Human skeletal muscle fibres: molecular and functional diversity. Prog Biophys Mol Biol. 2000;73(2–4):195–262. https://doi.org/10.1016/s0079-6107(00)00006-7.
Ebarb SM, Phelps KJ, Drouillard JS, Maddock-Carlin KR, Vaughn MA, Burnett DD, et al. Effects of anabolic implants and ractopamine-HCl on muscle fiber morphometrics, collagen solubility, and tenderness of beef longissimus lumborum steaks. J Anim Sci. 2017;95(3):1219–31. https://doi.org/10.2527/jas.2016.1263.
Scheffler TL, Leitner MB, Wright SA. Technical note: protocol for electrophoretic separation of bovine myosin heavy chain isoforms and comparison to immunohistochemistry analysis. J Anim Sci. 2018;96(10):4306–12. https://doi.org/10.1093/jas/sky283.
Bao G, Liu X, Wang J, Hu J, Shi B, Li S, et al. Effects of Slaughter Age on myosin heavy chain isoforms, muscle fibers, fatty acids, and Meat Quality in muscle of Tibetan Sheep. Front Vet Sci. 2021;8689589. https://doi.org/10.3389/fvets.2021.689589.
Immonen K, Ruusunen M, Hissa K, Puolanne E. Bovine muscle glycogen concentration in relation to finishing diet, slaughter and ultimate pH. Meat Sci. 2000;55(1):25–31. https://doi.org/10.1016/s0309-1740(99)00121-7.
Della Rosa MM, Pouzo LB, Pavan E. Meat and fat quality traits of grazing steers supplemented with corn grain and increasing amounts of flaxseed. Livest Sci. 2018;208:51–4. https://doi.org/10.1016/j.livsci.2017.12.004.
McGilchrist P, Alston CL, Gardner GE, Thomson KL, Pethick DW. Beef carcasses with larger eye muscle areas, lower ossification scores and improved nutrition have a lower incidence of dark cutting. Meat Sci. 2012;92(4):474–80. https://doi.org/10.1016/j.meatsci.2012.05.014.
Witt CC, Burkart C, Labeit D, McNabb M, Wu Y, Granzier H, et al. Nebulin regulates thin filament length, contractility, and Z-disk structure in vivo. EMBO J. 2006;25(16):3843–55. https://doi.org/10.1038/sj.emboj.7601242.
Conover GM, Gregorio CC. The desmin coil 1B mutation K190A impairs nebulin Z-disc assembly and destabilizes actin thin filaments. J Cell Sci. 2011;124(Pt 20):3464–76. https://doi.org/10.1242/jcs.087080.
Gokhin DS, Bang M-L, Zhang J, Chen J, Lieber RL. Reduced thin filament length in nebulin-knockout skeletal muscle alters isometric contractile properties. Am J Physiol Cell Physiol. 2009;296(5):C1123. https://doi.org/10.1152/ajpcell.00503.2008.
Bang M-L, Li X, Littlefield R, Bremner S, Thor A, Knowlton KU, et al. Nebulin-deficient mice exhibit shorter thin filament lengths and reduced contractile function in skeletal muscle. J Cell Biol. 2006;173(6):905–16. https://doi.org/10.1083/jcb.200603119.
Honikel KO. Reference methods for the assessment of physical characteristics of meat. Meat Sci. 1998;49(4):447–57. https://doi.org/10.1016/S0309-1740(98)00034-5.
Wang Y, Wang Z, Hu R, Peng Q, Xue B, Wang L. Comparison of carcass characteristics and meat quality between Simmental crossbred cattle, cattle-yaks and Xuanhan yellow cattle. J Sci Food Agric. 2021;101(9):3927–32. https://doi.org/10.1002/jsfa.11032.
Friedländer MR, Mackowiak SD, Li N, Chen W, Rajewsky N. miRDeep2 accurately identifies known and hundreds of novel microRNA genes in seven animal clades. Nucleic Acids Res. 2012;40(1):37–52. https://doi.org/10.1093/nar/gkr688.
Anders S, Huber W. Differential expression analysis for sequence count data. Genome Biol. 2010;11(10):R. https://doi.org/10.1186/gb-2010-11-10-r106.
Gao Y, Wang J, Zhao F. CIRI: an efficient and unbiased algorithm for de novo circular RNA identification. Genome Biol. 2015;16:4. https://doi.org/10.1186/s13059-014-0571-3.
Kanehisa M, Furumichi M, Sato Y, Kawashima M, Ishiguro-Watanabe M. KEGG for taxonomy-based analysis of pathways and genomes. Nucleic Acids Res. 2023;51(D1):D587. https://doi.org/10.1093/nar/gkac963.
Kanehisa M. Toward understanding the origin and evolution of cellular organisms. Protein Sci. 2019, 28(11):1947–1951.https://doi.org10.1002/pro.3715.
Kanehisa M, Goto S. KEGG: kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 2000;28(1):27–30. https://doi.org/10.1093/nar/28.1.27.
Bustin SA, Benes V, Garson JA, Hellemans J, Huggett J, Kubista M, et al. The MIQE guidelines: minimum information for publication of quantitative real-time PCR experiments. Clin Chem. 2009;55(4):611–22. https://doi.org/10.1373/clinchem.2008.112797.
We gratefully acknowledge the invaluable help from lab technicians of Gene Bank of Guangdong Livestock and Poultry in the sample collection process and providing computational resources in this work.
This research has received funding from the Key Realm R&D Program of Guangdong Province (2022B0202110002), Guangdong Basic and Applied Basic Research Foundation of China (2019B1515210020).
Our study was conducted in accordance with ARRIVE guidelines and no anesthetic or euthanasia agents were used in our study. All experimental animal procedures were in accordance with the relevant guidelines and regulations for the management and welfare of experimental animals approved by the Ethics Committee of South China Agricultural University. This study was approved by the Ethics Committee of South China Agricultural University in Guangdong Province.
Consent for publication
The authors declare no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Electronic supplementary material
Below is the link to the electronic supplementary material.
About this article
Cite this article
Yang, C., Wu, L., Guo, Y. et al. Expression profile and bioinformatics analysis of circRNA and its associated ceRNA networks in longissimus dorsi from Lufeng cattle and Leiqiong cattle. BMC Genomics 24, 499 (2023). https://doi.org/10.1186/s12864-023-09566-0