Abstract Background Intramuscular fat (IMF) content is a crucial determinant of beef quality and a key indicator in cattle breeding and production. However, the molecular regulatory mechanisms governing IMF deposition remain poorly understood. Results This study preliminarily explored the molecular mechanisms underlying IMF deposition by integrating weighted gene co-expression network analysis (WGCNA) and competitive endogenous RNA (ceRNA) network analysis. Sequencing of longissimus dorsi muscle samples from crossbred Wagyu cattle with varying IMF deposition levels revealed 172 differentially expressed circular RNAs (circRNAs), which were subsequently annotated and used to construct regulatory networks. Protein-protein interaction (PPI) network analysis predicted possible several lipid metabolism-related genes, including EZH2, AKT3, APP and SMARCA5. By combining the miRNA and mRNA data from our previous studies, we constructed circRNA-mRNA coexpression networks and circRNA-miRNA-mRNA regulatory networks. Functional enrichment analysis revealed that the identified circRNAs are involved primarily in lipid metabolism-related pathways, including phosphatidylinositol metabolism and the cGMP-PKG signaling pathway. Additionally, several circRNAs were predicted to function as molecular sponges based on coexpression patterns. Conclusion This study provides novel insights into the molecular mechanisms underlying IMF deposition in hybrid cattle and provides candidate regulatory mechanisms for further validation in selective breeding. Supplementary Information The online version contains supplementary material available at 10.1186/s12864-025-12097-5. Keywords: CircRNA, WGCNA, CeRNA network, Intramuscular fat, Cattle Introduction IMF is a specialized type of adipose tissue predominantly found within skeletal muscle and is composed of white adipocytes. These cells are arranged in a distinctive pattern, creating a marbled texture commonly known as marbling fat. The content of IMF plays a pivotal role in determining critical meat quality traits, such as juiciness, color, tenderness, and flavor [[46]1, [47]2]. Beef with a high IMF content exhibits a characteristic marbling pattern that is highly desirable in premium meat products. Among cattle breeds, Qinchuan cattle are renowned for their robust growth performance and adaptability, whereas Wagyu cattle are distinguished by their exceptional meat quality and superior IMF deposition [[48]3]. Notably, Wagyu cattle possess an extraordinary ability to accumulate substantial amounts of IMF compared with other breeds, leading to a significantly greater capacity for marbling fat production [[49]4]. Crossbreeding Qinchuan cattle with Wagyu cattle has been proposed as a strategy to combine the desirable traits of both breeds, aiming to increase beef quality, particularly by increasing the IMF content. This crossbreeding approach aims to increase IMF levels in Qinchuan beef while preserving the breed’s advantageous growth performance. Unlike most linear messenger RNAs (mRNAs) and long noncoding RNAs (lncRNAs), which have N7-methylguanosine (m7G) caps and 3’ polyadenylic acid tails, circRNAs are covalently closed, single-stranded RNA molecules [[50]5]. Recent studies investigating the role of circRNAs in IMF deposition have focused primarily on economically important species such as yaks and buffalo. Among these studies, peroxisome proliferator-activated receptor γ (PPARG) has been identified as a key regulator of IMF deposition. PPARG is positively correlated with IMF deposition and can inhibit the transcription of hormone-sensitive lipase (HSL), a lipolytic enzyme, through the generation of exon-derived circPPARG1. This circRNA interacts with PPARG proteins to increase lipogenic differentiation and lipid accumulation [[51]6]. In yak IMF adipocytes (YIMAs), circCWC22 is produced via reverse splicing of exons 5 to 10 from its source gene and has been shown to promote differentiation while inhibiting the proliferation of YIMAs through the circCWC22/miR-3059-x/HMGCL axis [[52]7]. Similarly, transcriptomic analyses of bovine adipocytes from calves and adult cattle revealed that circFUT10 promotes cell proliferation but inhibits adipogenic differentiation by sequestering let-7c, a microRNA that targets PPARG cofactor 1b (PPARGC1B). This regulation impacts the expression of key adipogenic markers, including PPARG and CCAAT/enhancer-binding protein β (C/EBPβ) [[53]8].In another study, circINSR was found to reduce the inhibitory effects of miR-15/16 on FOXO1 and EPT1, thereby suppressing adipogenic gene expression and lipid accumulation in preadipocytes [[54]9]. Additionally, overexpression of circFLT1 significantly decreased miR-93 levels in preadipocytes, potentially promoting adipocyte differentiation. Conversely, miR-93 overexpression downregulated the expression of critical adipogenic markers, such as PPARG and CEBPA, suggesting that circFLT1 enhances adipocyte differentiation by modulating miR-93 activity [[55]10]. WGCNA further identified circMARK3 as a promoter of adipogenic differentiation in buffalo adipocytes and 3T3-L1 cells, acting through the upregulation of the adipogenic marker genes PPARG, C/EBPα, and FABP4 [[56]11]. Collectively, these findings indicate that circRNAs play essential roles in fat metabolism and IMF deposition in bovine species, such as yaks and buffalo. These insights into the regulatory functions of circRNAs in adipogenesis provide a foundation for studying their potential roles in improving the IMF content in economically important crossbred cattle, such as Wagyu. While numerous studies have identified genes involved in the regulatory pathways of IMF deposition [[57]12], the molecular mechanisms underlying IMF deposition in crossbred Wagyu cattle remain poorly understood, and studies investigating the associated regulatory networks are limited. To address this gap, we collected longissimus dorsi muscle samples from crossbred Wagyu bovines, analyzed the expression of differential circRNAs, and performed bioinformatics analysis of the genes encoding these differentially expressed circRNAs (DECs) to elucidate their potential regulatory roles. In addition, on the basis of previous miRNA and mRNA sequencing data from the same samples in our group, we constructed a WGCNA and a ceRNA regulatory network to uncover key circRNA-miRNA-mRNA interactions associated with hybridization and bovine IMF deposition. These findings provide new insights into the molecular mechanisms of IMF regulation and provide a valuable basis for further research on the formation of IMF deposition. Materials and methods Sample collection and IMF content determination A total of 30 healthy crossbred Wagyu cattle (♀ Qinchuan cattle × ♂ Wagyu cattle), aged 28 ± 2 months and weighing 528 ± 77.39 kg, were selected from a farm in Beijing (Have obtained informed consent from the owner to use the animals in this study). All the cattle were reared under uniform nutritional and management conditions to minimize environmental variability. Before slaughter, the animals were fasted for 24 h with ad libitum access to water. Following standard slaughter procedures (The bovine euthanasia procedure employed a two-stage electrical stunning protocol in strict accordance with OIE Terrestrial Animal Health Code Chap. 7.5 of animal slaughter. This protocol ensures immediate unconsciousness, sustained unconsciousness until brain death, and minimization of the stress response), the longissimus dorsi muscle was collected, with a portion snap-frozen in liquid nitrogen and stored at −80 °C for RNA extraction. Approximately 500 g of the longissimus dorsi muscle was also collected for IMF content determination via the Soxhlet extraction method [[58]13]. On the basis of the IMF content measurements, three individuals with high IMF (H group; 32.96 ± 3.26%) and three individuals with low IMF (L group; 10.91 ± 1.85%) were selected for further analysis. The differences in the IMF content between the two groups were statistically significant (P < 0.05). The phenotypic characteristics and comparisons between the high- and low-IMF groups have been documented in detail in our previous study [[59]14]. RNA extraction, library construction, and sequencing Total RNA was extracted from longissimus dorsi muscle samples via TRIzol Reagent (Invitrogen, Carlsbad, CA, USA) following the manufacturer’s protocol. The quality and integrity of the RNA were evaluated via an Agilent 2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA, USA) and RNase-free agarose gel electrophoresis. The RNA samples meeting the quality standards were further processed. Linear RNA was selectively degraded by RNase R treatment, and the remaining RNA was purified with the RNeasy MinElute Cleanup Kit (Qiagen, Venlo, The Netherlands). Strand-specific RNA libraries for Illumina sequencing were constructed following the manufacturer’s instructions via the VAHTS Total RNA-Seq (H/M/R) Library Preparation Kit (Vazyme, Nanjing, China). Ribosomal RNA (rRNA) was first removed to enrich circRNAs, which were subsequently fragmented into short segments via fragmentation buffer. Complementary DNA (cDNA) was synthesized through reverse transcription with random primers. cDNA second-strand synthesis was performed using second-strand synthesis reagents (Vazyme). The resulting cDNA fragments were purified via VAHTSTM DNA Clean Beads (Vazyme, Nanjing, China) for end-repair, A-tailing, and ligation of Illumina sequencing adapters. To ensure strand specificity, the second strand of cDNA was digested with uracil-N-glycosylase (UNG), and the digested product was purified via VAHTSTM DNA Clean Beads. Finally, the prepared Libraries were subjected to PCR amplification and sequenced on the Illumina NovaSeq 6000 platform. Quality control and circrna identification To ensure data quality for subsequent analyses, raw reads generated via Illumina sequencing were processed to remove adapters and low-quality sequences. Specifically, raw reads containing adapter sequences, low-quality bases, or ambiguous nucleotides were filtered via fastp (v0.23.4) [[60]15]. The filtering criteria were as follows: (i) removal of reads containing adapter sequences; (ii) exclusion of reads with more than 10% ambiguous bases (N); and (iii) removal of reads where more than 50% of the bases had a Phred quality score ≤ 20. This step ensured the generation of high-quality clean reads for downstream analyses. Given that the efficiency of rRNA removal may vary due to differences in sample quality and species-specific variations, rRNA contamination was assessed and removed via Bowtie2 (v2.2.8) by mapping the reads to an rRNA database [[61]16]. Reads mapped to the rRNA database were discarded, and the remaining nonrRNA reads were used for subsequent analyses. The filtered reads were then mapped to the reference genome via TopHat2 (v2.1.1) with default parameters [[62]17]. Reads that aligned to the reference genome were discarded, while the unmapped reads were subjected to further analysis for circRNA identification. To identify candidate circRNAs, the terminal 20 nucleotides from both ends of each unmapped read were extracted and aligned to the reference genome to identify unique anchor positions spanning back-splicing junctions. Reads showing reverse (head-to-tail) orientation at the splice sites were identified as potential circRNA candidates. CircRNA identification was performed via find_circ (v1.2), which extends the anchored alignment to ensure that the entire read aligns to the genome while verifying the presence of canonical splice sites (GU/AG motifs) at the junctions. CircRNAs supported by at least two uniquely aligned back-spliced reads in at least one sample were classified as candidate circRNAs and retained for downstream analyses. This systematic approach ensured the identification of high-confidence circRNA candidates, providing a robust dataset for subsequent characterization and functional studies. Quantification of circrnas and differential expression analysis To quantify circRNA abundance, the number of back-spliced junction reads was normalized to the number of reads per million mapped reads (RPM), ensuring comparability across samples with varying sequencing depths. The RPM for each circRNA was calculated via the following formula: graphic file with name d33e550.gif In this equation, C represents the number of uniquely mapped back-spliced junction reads corresponding to circRNAs, whiereas N denotes the total number of back-spliced junction reads. The RPM normalization method effectively mitigates the impact of varying sequencing depths on circRNA expression calculations, allowing for direct comparisons of differential expression across samples. To identify DECs between experimental groups, the edgeR software package (v 3.12.1) was employed. CircRNAs exhibiting a |log[2]FoldChange| ≥ 1 and a P-value < 0.05 were classified as significantly differentially expressed. This rigorous statistical approach ensures the reliability of the identified circRNAs, facilitating further biological interpretation and functional analysis. Construction of the PPI network and identification of key genes PPI network analysis was conducted to explore the functional associations between proteins encoded by DEC-associated genes. Interaction data were obtained from the STRING database, a comprehensive resource for functional protein association networks. Confidence interaction scores (≥ 0.4) were utilized to ensure the reliability of the network. The PPI network was visualized and constructed via Cytoscape software (v3.10.3), a widely used tool for network biology. To identify key genes within the network, the CytoNCA plug-in (v2.1.6) was employed to perform topological analysis on the basis of centrality measures. Specifically, betweenness centrality (BC) was used as the primary metric for evaluating the importance of each node. Genes with BC values exceeding the median threshold were prioritized as hub genes, reflecting their potential significance in maintaining network structure and functionality. The top-ranked hub genes were further analyzed to investigate their biological relevance in the context of the study. CircRNA‒mRNA coexpression network analysis To explore the regulatory relationships between circRNAs and mRNAs, WGCNA was performed by integrating circRNA and mRNA expression data derived from the same samples. The analysis was based on total RNA sequencing data, which included circRNAs and coding gene mRNAs identified in all six samples. Differentially expressed genes were selected on the basis of the following thresholds: false discovery rate (FDR) < 0.05 and |log[2]FoldChange| ≥1. Coexpression networks were constructed via the one-step modular algorithm of the WGCNA package (v1.73) in R software (v4.4.0). The network parameters were optimized as follows: R2 = 0.9, networkType = Signed, corFnc = bicor, minModuleSize = 30, and mergeCutHeight = 0.25. These settings ensured the accurate identification of modules with distinct coexpression patterns. Modules significantly correlated with IMF were identified, and the mRNAs within these modules were further analyzed for functional enrichment. Gene Ontology (GO) molecular function analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis were performed on the identified mRNAs significantly associated with IMF via the DAVID database. The results were visualized to highlight enriched biological processes, molecular functions, and pathways, providing insights into the potential roles of circRNAs in regulating gene expression and their contributions to cellular functions. Construction of a circrna‒mirna‒mrna competitive endogenous RNA network To explore the regulatory role of circRNAs in IMF deposition, a circRNA-miRNA‒mRNA ceRNA network was constructed on the basis of circRNAs identified in the WGCNA modules significantly associated with IMF. First, normalized expression data for miRNAs were obtained from small RNA sequencing of the same six samples. A total of 564 miRNAs with detectable expression in these samples were included in the analysis. Second, the TargetScan database was used to identify miRNA target sites in circRNA sequences, and only the predicted circRNA‒miRNA interactions were retained for further analysis. Third, miRNA‒mRNA interactions were inferred via both the TargetScan and the miRDB databases, with the analysis limited to mRNAs present in the coexpression modules identified via WGCNA. Common miRNAs capable of binding both the circRNAs and these module-specific mRNAs were identified, enabling the establishment of connections between the circRNAs and their downstream mRNAs via shared miRNAs. Fourth, Spearman rank correlation analysis was performed on the normalized expression data of the circRNAs, miRNAs, and mRNAs. Interactions were retained only if they demonstrated significant negative correlations for circRNA‒miRNA and miRNA‒mRNA pairs, as well as positive correlations for circRNA‒mRNA pairs (Spearman’s P < 0.05). This correlation analysis ensured the reliability of the predicted ceRNA relationships. Finally, the maximum clique computation (MCC) algorithm within the CytoHubba plug-in of Cytoscape software was used to identify the top 60 nodes within the ceRNA network. After independent nodes were excluded, the final ceRNA regulatory network, which represented key interactions among circRNAs, miRNAs, and mRNAs, was constructed. This ceRNA network highlights the potential roles of circRNAs as regulators of miRNA-mediated gene expression in IMF deposition. These findings provide a robust foundation for further experimental validation and a deeper understanding of the molecular mechanisms underlying IMF accumulation. Real-time quantitative PCR validation To validate the accuracy of the sequencing results, five circRNAs were randomly selected for quantification via real-time quantitative PCR (RT-qPCR). Total RNA was extracted from the samples and assessed for integrity via agarose gel electrophoresis. The extracted RNA was then reverse-transcribed into cDNA via the PrimeScript RT Kit (Takara, Dalian, China) according to the manufacturer’s protocol. Primers specific to the selected circRNAs were designed via Primer5.0 software (v5.5.0) with GAPDH used as the internal reference gene. The primers were synthesized by Sangon Bioengineering Co., Ltd. (Shanghai, China). The detailed sequences of the primers are provided in Supplementary Table S1. The qPCRs were performed using TB Green Premix Ex Taq II (Takara, Dalian, China) as the fluorescent dye reagent. Amplification was carried out on a CFX96 Real-Time PCR Detection System (Bio-Rad, Hercules, CA, USA) using the following program: initial denaturation at 95 °C for 30 s, followed by denaturation at 95 °C for 5 s and annealing at 60 °C for 30 s, for a total of 40 cycles. Melting curve analysis was conducted to confirm the specificity of the amplification, with conditions set at 60 °C for 1 min and 95 °C for 15 s. For detailed information on the qPCR system, refer to Supplementary Table S2. Each sample were replicateed three times, and the relative expression levels of the selected circRNAs were calculated via the 2^−ΔΔCt method. The expression levels were normalized and further analyzed via log[2]FoldChange to ensure the accuracy of the sequencing data. Results Summary of RNA-Seq data High-throughput RNA sequencing generated 306,863,386 and 275,544,550 clean reads from the high- and low-IMF groups, respectively. Approximately 93.94% of the clean reads were successfully aligned to the Bos taurus reference genome(ARS-UCD2.0). The sequencing quality was high, with Q30 values ranging from 92.69 to 93.56%, and the GC content of the reads ranged from 53.84 to 56.73%. Among the reads successfully aligned to the genome, an average of 88.96% were uniquely mapped and subsequently used for transcript reconstruction and downstream analyses. Detailed statistics of the sequencing data are provided in Supplementary Table S3. Identification and characterization of circrnas To identify circRNAs, the positional distribution of reads mapped to the Bos taurus reference genome was analyzed on the basis of the total mapped reads. As shown in Fig. [63]1A, the majority of reads were aligned to exonic regions, followed by intronic and intergenic regions. A total of 10,382 circRNAs were identified, and their chromosomal distribution was mapped (Fig. [64]1B). Among these, chromosome 2 presented the greatest number of circRNAs (937), whereas chromosome 25 presented the lowest number(190). The abundance of circRNAs was further analyzed on the basis of the number of backspliced reads (Fig. [65]1C). CircRNAs with 1–2 backspliced reads were the most prevalent, comprising approximately 20% of the total identified circRNAs. The overall abundance of circRNAs was found to be low, with more than half of the identified circRNAs having fewer than 10 backspliced reads. CircRNA annotation revealed that the majority were derived from exonic regions. The length distribution of the identified circRNAs was also examined (Fig. [66]1D), revealing that circRNAs with lengths between 300 and 400 bp were the most abundant (1,451 circRNAs), accounting for approximately 14% of the total. Fig. 1. [67]Fig. 1 [68]Open in a new tab Identification and characterization of circRNAs. A: Distribution of sequenced reads across genomic regions. B: Chromosomal distribution of identified circRNAs. C: Abundance analysis of circRNA reads. D: Length distribution of identified circRNAs Differential expression of circrnas To investigate the potential role of circRNAs in regulating IMF deposition, the expression profiles of circRNAs in the dorsal longissimus dorsi muscle of crossbred Wagyu cows with high and low IMF contents were analyzed. A total of 172 DECs were identified via the filtering criteria of a |log[2]FoldChange| ≥ 1 and a P-value < 0.05. The results are presented as a volcano plot and a clustered heatmap (Fig. [69]2A and B). Among these, 89 circRNAs were upregulated, whereas 83 were downregulated. Notably, approximately 79% of the DECs presented a |log[2]FoldChange | >2. Among the upregulated DECs, 68 circRNAs were derived from annot_exons, whereas only 2 originated from intronic regions. Similarly, in the downregulated DECs, 58 circRNAs were derived from annot_ exons, while both intronic and antisense regions contributed a minimum of 2 circRNAs. The detailed genomic distribution of these DECs is provided in Supplementary Table S4. Annotation of these 172 DECs revealed that they were derived from a total of 158 host genes (Supplementary Table S5). Additionally, differential expression analyses of mRNAs and miRNAs from the same samples revealed 359 differentially expressed mRNAs, of which 297 were upregulated and 62 were downregulated (Fig. [70]2C). Furthermore, 38 differentially expressed miRNAs were detected, including 18 upregulated and 20 downregulated miRNAs (Fig. [71]2D). Fig. 2. [72]Fig. 2 [73]Open in a new tab Screening and characterization of DECs. A: Volcano plot of DECs. B: Cluster heatmap of DECs. C: Number of differential mRNAs expressed up- and downregulated. D: Number of differential miRNAs expressed up- and downregulated Functional enrichment analysis of the source genes of DECs To investigate the potential biological roles of circRNAs associated with varying IMF contents, GO and KEGG enrichment analyses were conducted for the host genes of DECs. GO analysis revealed that these host genes were significantly enriched in 24 biological processes, predominantly involving cellular processes, regulatory pathways, and metabolic processes; 15 cellular components, including cells, cell parts, and organelles; and 7 molecular functions, such as binding, catalytic activity, and regulation of molecular function (Fig. [74]3A). KEGG pathway enrichment analysis revealed the top 20 pathways ranked by significance, with several adipose-related pathways including phosphatidylinositol metabolism and the cGMP-PKG signaling pathway (Fig. [75]3B). Fig. 3. [76]Fig. 3 [77]Open in a new tab Functional enrichment analysis of DEC-derived genes. A: GO enrichment analysis of DEC-derived genes. B: KEGG enrichment analysis of DEC-derived genes PPI network construction and identification of key genes The protein interactions of 158 genes derived from DECs were analyzed, with an interaction score threshold set at 0.4. A total of 154 genes were identified to participate in 126 interactions, forming a preliminary network consisting of 154 nodes and 126 edges. To refine this network, the CytoNCA plugin was utilized to evaluate the topological properties of the network, specifically focusing on the BC scores. On the basis of the high BC scores, a more condensed PPI network comprising 72 nodes and 97 edges was constructed (Fig. [78]4). Several genes with the highest BC values (genes with BC scores exceeding 900; Supplementary Table S6) were identified as potentially critical for bovine IMF deposition. These genes included EZH2, AKT3, MYSM1, GNAQ, BRAF, APP, and SMARCA5. Notably, EZH2, APP, and SMARCA5 have been previously reported to be involved in fat metabolism and the regulation of lipids, as well as other key biological processes [[79]18–[80]20]. Fig. 4. [81]Fig. 4 [82]Open in a new tab PPI network construction of DEC-derived genes. Each node represents a protein encoded by a single gene, with the node size corresponding to the BC value. Larger nodes indicate higher BC values. Pink nodes represent BC values in the top 20, and cyan nodes represent BC values in the bottom 50 Coexpression network of circrnas and mRNAs To further investigate the relationships between circRNAs and mRNAs, we constructed a weighted gene co-expression network via two RNA-seq datasets from the same transcripts. This analysis aimed to identify key mRNAs potentially associated with IMF traits. On the basis of the hypothesis that coexpressed genes are Likely involved in the same or related biological processes, we inferred the potential functions of circRNAs through their coexpressed mRNAs. A soft threshold of 13 was determined by a scale-free topology fit index of 0.9 (Fig. [83]5A). Using this threshold, 10 distinct gene expression modules were identified, and a clustered dendrogram of gene expression modules was constructed (Fig. [84]5B). Each module was summarized by its module eigengene (ME), which represents the first principal component of circRNA and mRNA expression in that module. To further explore the relationships between modules, we performed a correlation analysis and visualized the results via a heatmap (Fig. [85]5C). Finally, module‒trait associations were analyzed to identify modules significantly correlated with IMF traits. Among the identified modules, two gene modules, MEblue and MEblack, showed significant associations with IMF (Fig. [86]5D), highlighting their potential relevance in the regulation of IMF deposition. Fig. 5. [87]Fig. 5 [88]Open in a new tab Analysis of the WGCNA coexpression network. A: Demonstration of the degree of scale-free topology and module connectivity for different soft thresholds via power function exponentials. B: Dendrogram of gene module division when a soft threshold of 13 was determined. C: Correlations between modules. D: Correlations between modules and IMF deposition traits To investigate the functional implications of mRNAs significantly correlated with sample groups in each module (FDR < 0.05), we performed GO molecular function and KEGG pathway enrichment analyses. In the MEblue module, genes were filtered using stringent criteria of |GS| >0.9 (gene significance) and |MM| >0.9 (module membership) (Supplementary Fig. [89]S1A). Owing to the smaller number of genes in the MEblack module, the selection thresholds were adjusted to |GS| >0.85 and |MM| >0.85 (Supplementary Fig. [90]S1B). As a result, 158 mRNAs were identified in the MEblue module(Supplementary Table S7), whereas 58 mRNAs were identified in the MEblack module (Supplementary Table S8). To highlight the biological significance of these modules, the top 10 GO terms and KEGG pathways based on P-values were selected for visualization (Fig. [91]6). For the MEblue module, GO enrichment analysis revealed significant enrichment in biological processes such as actin cytoskeleton organization, lipid binding, and actin filament assembly (Fig. [92]6A). KEGG pathway analysis revealed that the MEblue module was associated with pathways including beta-alanine metabolism, aldosterone synthesis and secretion, the cAMP signaling pathway, and general metabolic pathways (Fig. [93]6B). In the MEblack module, GO enrichment analysis revealed processes such as fibroblast migration and actin binding (Fig. [94]6C). KEGG pathway analysis further revealed enrichment of pathways related to insulin resistance, cytoskeletal organization in muscle cells, and other relevant metabolic pathways (Fig. [95]6D). Fig. 6. [96]Fig. 6 [97]Open in a new tab GO and KEGG pathway enrichment analysis for the MEblue and MEblack modules. A, B: GO terms and KEGG pathways associated with genes in the MEblue module; C, D: GO terms and KEGG pathways associated with genes in the MEblack modules CircRNA‒miRNA‒mRNA regulatory network construction CircRNAs have emerged as crucial regulators of gene expression, often acting as miRNA sponges that repress miRNA activity and consequently regulate their target mRNAs. To identify circRNAs with potential miRNA sponging activity in the context of high and low IMF contents, we integrated circRNA and mRNA expression profiles with miRNA sequencing data from the same samples. This allowed for the construction of a circRNA-miRNA‒mRNA regulatory network. The relationships between circRNAs, miRNAs, and mRNAs were determined on the basis of significant negative correlations between circRNA-miRNA and miRNA‒mRNA interactions, along with significant positive correlations between circRNA‒mRNA interactions (Spearman’s P-value < 0.05) (Fig. [98]7A). In the coexpression network modules significantly correlated with IMF content, we identified positively correlated circRNAs and mRNAs that may compete with shared miRNAs (Supplementary Table S9). The constructed circRNA‒miRNA‒mRNA network consisted of 207 differentially expressed genes (DEGs), 70 DECs, 57 differentially expressed miRNAs (DEMs), and 1,435 interaction edges (Fig. [99]7B). This analysis identified 70 circRNAs with potential miRNA sponging functions, explaining the observed coexpression patterns between circRNAs and mRNAs associated with the IMF content. For example, novel_circ_009207 and its target gene NFIC are both regulated by bta-miR-150. Similarly, novel_circ_001732, novel_circ_003151, novel_circ_007171, and novel_circ_000946—along with their shared target genes NFIC and CKRL—are all regulated by bta-miR-204. Moreover, NFIC is targeted by both bta-miR-150 and bta-miR-204. Previous studies have implicated NFIC, bta-miR-150, and bta-miR-204 in adipocyte differentiation and other regulatory processes [[100]21–[101]23]. These findings suggest that NFIC is a candidate regulatory gene potentially involved in IMF deposition, with circRNAs potentially mediating its regulation through their miRNA sponge activity. To further investigate this regulatory relationship, we constructed a refined ceRNA network by visualizing the top 60 genes via the CytoHubba algorithm (Fig. [102]7C). This subnetwork comprised 16 DEGs, 7 DECs, 29 DEMs, and 113 interaction edges. The analysis revealed that these seven circRNAs are prioritized candidates for functional studies on adipogenesis regulation. These findings highlight the importance of circRNA‒miRNA‒mRNA networks in IMF deposition and adipogenesis, providing novel insights into the molecular mechanisms underlying these processes. Fig. 7. [103]Fig. 7 [104]Open in a new tab CircRNA-miRNA‒mRNA regulatory network construction and analysis. A: Schematic representation of the circRNA-miRNA‒mRNA regulatory network. B: Comprehensive circRNA-miRNA‒mRNA network. C: The circRNA‒miRNA‒mRNA network of the top 60 genes. Node color is directly proportional to the maximum community centrality (MCC) score. The darker the node color, the higher its MCC value RT‒qPCR to verify the expression of differential genes Five DECs were randomly selected, and their expression levels were detected via RT‒qPCR. The expression levels obtained from the experiment and sequencing were log[2]FoldChange normalized (averaged in three replicates). The RT‒qPCR results were consistent with the RNA-Seq data (Fig. [105]8A). Furthermore, the correlation coefficient between the RNA-Seq and RT‒qPCR data was 0.99 (Fig. [106]8B), which was highly correlated. The above results confirm that the circRNA analysis results on the basis of RNA-Seq data are relatively reliable. Fig. 8. [107]Fig. 8 [108]Open in a new tab RT-qPCR to verify differentially expressed genes. A: The expression levels of DECs were validated by RT-qPCR. B: Comparison between RNA-Seq and RT-qPCR by correlational analyses Discussion The accumulation of IMF is a critical factor affecting beef quality. Beef cattle with varying IMF deposition levels exhibit distinct gene expression patterns associated with this trait. IMF formation is a complex process influenced by multiple factors, and understanding the underlying genetic mechanisms of the IMF is crucial. Recent studies have highlighted the regulatory roles of noncoding RNAs, particularly circRNAs, in this process. Previous analyses of muscle and adipose tissues from cattle of different ages have constructed regulatory networks of mRNAs and non-coding RNAs associated with IMF deposition [[109]24]. These studies have identified key regulators and provided insights into their roles. For example, circSSBP2, a novel circRNA, has been shown to influence the proliferation of bovine intramuscular preadipocytes. Knockdown of circSSBP2 promotes preadipocyte proliferation, whereas its overexpression inhibits it [[110]25]. Similarly, circPPARA was found to bind to miR-429 and miR-200b, facilitating intramuscular adipogenesis in pigs [[111]6]. These findings underscore the importance of circRNAs in the regulation of IMF deposition. In this study, a PPI network was constructed using genes derived from DECs, leading to the identification of key candidate genes that may be associated with adipose tissue, including EZH2, AKT3, APP, and SMARCA5. CircEZH2, a circRNA transcribed from the EZH2 gene, selectively inhibits EZH2 expression, resulting in a reduction in lipid accumulation in Huh-7 cells [[112]26]. Furthermore, circEZH2 regulates milk lipid metabolism via its sponge activity on miR-378b, and its overexpression promotes the expression of fatty acid transport- and metabolism-related proteins, such as CD36, LPL, FADS1, and SCD1, underscoring its regulatory role in lipid metabolism [[113]27]. In addition, AKT3 was identified as a key regulator of adipogenesis in mouse models, where it was found to inhibit alternative pathways of adipogenesis by modulating WNK lysine-deficient protein kinase-1 (WNK1) and serum/glucocorticoid-regulated kinase 1 (SGK1), thus playing a pivotal role in the regulation of lipid synthesis and metabolism [[114]28]. Another critical gene, APP, which is primarily recognized for its role in the pathogenesis of Alzheimer’s disease, was observed to be highly expressed in adipose tissue and upregulated in the adipocytes of obese individuals, with its expression levels positively correlated with fat deposition [[115]29, [116]30]. Consistent with these findings, our results revealed that APP expression was significantly greater in the high-IMF group than in the low-IMF group, further suggesting a potential positive association. Finally, SMARCA5, a chromatin remodeling factor, was implicated in adipogenesis through its interaction with Werner syndrome protein (WRN), which regulates adipogenesis in Werner syndrome [[117]20, [118]30]. While SMARCA5 was previously linked to body size traits in Hulunbeier sheep [[119]31], our results suggest a distinct role in adipogenesis, possibly due to species-specific regulatory mechanisms. These discrepancies highlight the need for further investigation to clarify the precise role of SMARCA5 in adipogenesis. Enrichment analysis of GO and KEGG pathways for DEC-derived genes revealed significant enrichment in pathways such as phosphatidylinositol metabolism and the cGMP-PKG signaling pathway. Phosphatidylinositol metabolism plays a pivotal role in adipocyte differentiation and lipid storage by regulating cell signaling pathways involved in lipid metabolism. Notably, inositol phosphates, such as IP3, promote the release of intracellular calcium ions in adipocytes, thereby influencing fatty acid synthesis and catabolism [[120]32]. Additionally, IP6K1, a key metabolite in the inositol phosphate metabolic pathway, and inositol pyrophosphate biosynthesis exert critical effects on lipid metabolism, as evidenced by studies showing that IP6K1 deficiency leads to increased basal lipolysis [[121]33]. The cGMP‒PKG signaling pathway, which is mediated by guanylate cyclase (GC) and catalyzed by cyclic guanosine monophosphate (cGMP), can be activated by nitric oxide (NO), subsequently activating PKG, which promotes fatty acid oxidation and inhibits lipid synthesis [[122]34]. Moreover, this pathway has been implicated in the adipogenic differentiation of adipose-derived stem cells and may interact with other signaling pathways to further modulate adipogenesis and lipid catabolism [[123]35, [124]36]. Pathways identified in enrichment analyses of WGCNA-related modules included cAMP signaling, aldosterone synthesis and secretion, lipid binding, insulin resistance, and metabolic pathways. The cAMP signaling pathway functions through cyclic adenosine monophosphate (cAMP) as a second messenger, which activates cAMP-dependent protein kinase (PKA). PKA phosphorylates downstream target proteins and activates CREB, which in turn induces the expression of C/EBPβ to promote adipogenesis [[125]37]. Additionally, exchange proteins directly activated by cAMP (EPAC) synergize with PKA to promote adipocyte differentiation and lipid accumulation [[126]38]. Aldosterone has been shown to influence lipogenesis and lipolysis by regulating insulin sensitivity and promoting adipocyte differentiation. Insulin, the primary hormone responsible for lowering blood glucose levels, also facilitates the synthesis of energy storage molecules such as fat, glycogen, and proteins, while aldosterone is hypothesized to modulate energy metabolism through its effects on insulin sensitivity [[127]39–[128]41]. However, canonical pathways frequently reported to regulate bovine IMF deposition, including the Wnt signaling pathway [[129]42], MAPK signaling pathway [[130]43], and AMPK signaling pathway [[131]44], were not significantly enriched in this study. The reasons underlying these discrepancies remain unclear and warrant further investigation. CeRNAs represent a class of RNA molecules, including mRNAs, lncRNAs, and circRNAs, that competitively bind to microRNAs (miRNAs), thereby indirectly regulating mRNA expression through modulation of miRNA activity [[132]45]. In the ceRNA regulatory network constructed in this study, 70 circRNAs with potential roles as molecular sponges were identified, although their specific regulatory mechanisms require further validation. Among them, PPARA, NFIC, bta-miR-320a, bta-miR-204, and bta-miR-150 have previously been reported to be involved in lipid storage and adipocyte differentiation [[133]21–[134]23, [135]46, [136]47]. Nuclear factor IC (NFIC) has been shown to regulate adipocyte differentiation through modulation of the classical Wnt signaling pathway, thereby balancing the lipogenic and osteogenic differentiation of progenitor cells [[137]21]. Bta-miR-204 has been identified as a key regulator of bovine adipocyte differentiation and fatty acid metabolism [[138]22, [139]48], while overexpression of bta-miR-150 has been demonstrated to upregulate the mRNA expression of marker genes and increase the expression of proteins such as CDK1, CDK2, and PCNA [[140]23]. This overexpression promotes adipocyte proliferation, inhibits adipocyte differentiation, and reduces lipid droplet formation. In the ceRNA network constructed in this study, novel_circ_009207 was predicted to interact with bta-miR-150, potentially modulating NFIC expression, whereas novel_circ_001732, novel_circ_003151, novel_circ_007171, and novel_circ_000946 were predicted to regulate NFIC through bta-miR-204. However, the specific mechanisms by which novel_circ_009207, novel_circ_001732, novel_circ_003151, novel_circ_007171, and novel_circ_000946 interact with bta-miR-204 and bta-miR-150 remain unclear, and the literature on the regulation of NFIC by bta-miR-204 and bta-miR-150 is Limited. Furthermore, within the ceRNA network constructed for the top 60 genes, seven circRNAs were identified as potentially associated with IMF deposition. These findings underscore the need for further investigation into these circRNAs to elucidate their regulatory roles and mechanisms in IMF deposition. Despite certain limitations, such as the small sample size in the construction of WGCNA networks, which may affect the reliability of the findings, the overall study design remains sound and methodologically rigorous. Through integrative analysis of miRNAs and circRNAs, this study identified their critical roles within the ceRNA regulatory network, offering a comprehensive perspective on the molecular mechanisms underlying fat metabolism and IMF deposition. To increase the statistical reliability and credibility of future analyses, increasing the sample size is essential. Furthermore, functional validation experiments, including gene knockout, gene overexpression, and dual-luciferase reporter assays, should be conducted for key circRNAs, miRNAs, and target genes to verify the predicted regulatory interactions and elucidate their biological significance. These efforts will provide deeper insights into the ceRNA network and its contribution to fat metabolism and IMF deposition. Conclusion This study preliminarily deciphers the circRNA-mediated regulatory landscape governing IMF deposition in crossbred Wagyu cattle through integrative transcriptomic analysis. By synergizing WGCNA and ceRNA network approaches, we identified critical circRNA‒miRNA‒mRNA regulatory axes, including novel_circ_009207‒bta-miR-150‒NFIC and novel_circ_001732‒bta‒miR-204-NFIC, which orchestrate adipocyte differentiation and lipid metabolism. Functional enrichment and PPI network analyses further revealed that lipid-associated pathways, such as phosphatidylinositol metabolism and cGMP-PKG signaling, are central to IMF regulation, with hub genes (EZH2, AKT3, APP, and SMARCA5) serving as pivotal nodes in adipogenesis. Crucially, this work establishes circRNA as dynamic regulators of bovine IMF deposition, bridging transcriptional and posttranscriptional control mechanisms. These findings not only offer new molecular insights into beef marbling but also provide a framework for precision breeding strategies targeting lipid metabolism. Future studies should prioritize the functional validation of candidate circRNAs and their interactions across diverse cattle populations to translate these insights into actionable genetic markers. Supplementary Information [141]Supplementary material 1.^ (415.1KB, zip) Acknowledgements