Exploration of non-coding RNAs related to intramuscular fat deposition in Xinjiang Brown cattle and Angus × Wagyu cattle

Abstract

Non-coding RNAs (ncRNAs) serve as crucial regulatory elements in the process of adipogenesis in animals; however, the specific roles and interrelationships of ncRNAs in bovine fat deposition remain poorly understood. This study aims to investigate the differentially expressed ncRNAs in the longissimus dorsi muscle of Xinjiang Brown cattle (XB) and Angus $\times$ Wagyu cattle (AW), to elucidate the regulatory mechanisms underlying lipidogenesis that may involve ncRNAs. Four Xinjiang Brown cattle and four Angus $\times$ Wagyu cattle were selected, ensuring they are subjected to identical feeding conditions, in order to evaluate the intermuscular fat (IMF) of longissimus dorsi muscles. The fat content of muscle tissue was quantified using the Soxhlet extraction method, revealing that the fat levels in the AW group were significantly elevated compared to those in the XB group. Taking muscle samples for paraffin sectioning and observing their morphology, it was found that the fat richness of the AW group was significantly higher than that of the XB group. Utilizing high-throughput RNA sequencing technology, we conducted an extensive transcriptomic analysis of longissimus dorsi muscles of XB and AW to identify significant ncRNAs implicated in fat metabolism and adipogenesis. The miRNA analysis yielded between 109,343,831 117,258,570 clean reads, whereas the lncRNA and circRNA analyses produced between 81,607,756 102,917,174 clean reads. Subsequent analysis revealed the identification of 53 differentially expressed miRNAs, 176 differentially expressed lncRNAs, and 234 differentially expressed circRNAs. KEGG enrichment analysis revealed that the target genes of differentially expressed miRNAs, lncRNAs, and circRNAs are significantly enriched in 2, 17, and 22 distinct pathways, respectively. The pathways associated with the differential enrichment of miRNA target genes involve processes such as phosphorylation and protein modification. Concurrently, the pathways linked to the varying enrichment of lncRNA target genes encompass G protein-coupled receptor signaling, regulation of cell death and apoptosis, activities related to GTPase activation, and functions governing nucleotide triphosphatases, among others. The circRNA exhibiting differential expression are significantly enriched in a variety of biological processes, including signal transduction, nucleic acid synthesis, cellular architecture, GTPase activation, and phosphatase activities, among others. The analysis of the ncRNA interaction network suggests that AGBL1, THRB, and S100A13 may play pivotal roles in the formation and adipogenic differentiation of adipocytes. In conclusion, we conducted a comprehensive analysis and discussion of the complete transcriptome of intermuscular fat tissue from the longissimus dorsi muscles in Xinjiang Brown cattle and Angus $\times$ Wagyu cattle. This study provides a theoretical foundation for enhancing our understanding of the molecular mechanisms underlying fat metabolism and deposition in beef cattle.

Background

The assessment of beef quality involves the examination of its intramuscular fat content and marbling, which are closely interrelated. Marbling serves as an indicator of the richness of intramuscular fat, while the fat content directly influences the sensory characteristics and palatability of beef. An optimal level of fat can enhance tenderness, juiciness, and contribute to a more complex flavor profile in the meat, thereby positively impacting the overall quality of beef [1]. Nevertheless, the modern beef industry continues to encounter a significant challenge in optimizing intramuscular fat deposition. The process of adipogenesis consists of two stages: (1) The differentiation of pluripotent stem cells into adipocyte precursors. (2) Preadipocytes proliferate through multiple rounds and accumulate lipids in response to specific stimuli, ultimately undergoing differentiation into mature adipocytes [2]. Hence, it is imperative to comprehensively comprehend and precisely delineate the molecular mechanisms accountable for beef fat accumulation to ensure the generation of premium-quality beef [3].

The Xinjiang Brown Cattle has been developed through crossbreeding, incorporating a small number of Alauta cattle and Costrom cattle as the male parent, and local Kazakh cattle from Xinjiang as the female parent. This breed exhibits resilience in challenging environmental conditions and with coarse feed, along with resistance to freezing temperatures. Furthermore, it is renowned for its high-quality meat and milk production, has become the dominant breed in the cattle industry of Xinjiang region. However, the intramuscular fat content of Xinjiang Brown Cattle is not high, resulting in an inferior flavor profile. Wagyu and their hybrids are renowned for its high levels of intramuscular fat, and in comparison to other beef cattle breeds, it exhibits significantly increased marbling [4, 5]. Elucidating the factors contributing to the differential intramuscular fat deposition between Wagyu cattle and Xinjiang Brown cattle may aid in a more comprehensive understanding of the fat formation process.

In recent years, transcriptome sequencing technology has been widely utilized for the identification of genes exhibiting differential expression and for the exploration of novel transcription factors, thereby facilitating more comprehensive research on ncRNAs [6,7,8,9,10]. NcRNAs mainly comprise microRNAs (miRNAs), long noncoding RNAs (lncRNAs) and circular RNAs (circRNAs). The process of adipose tissue deposition is highly intricate, and there is mounting evidence indicating that ncRNAs plays a pivotal role in the differentiation and metabolism of adipocytes [11,12,13,14,15]. Research has demonstrated that miRNA is involved in the differentiation of human adipocytes and the regulation of lipid metabolism [16, 17]; Guo et al. reported increased expression levels of miR-142-3p, miR-379, and miR-196a in the adipose tissue of Wagyu cattle relative to those in Holstein cattle [18]; Yang et al. discovered that miR-330 targets SESN3 to activate the Akt-mTOR pathway, promoting adipogenesis in bovine intramuscular preadipocytes [19]. Zhang et al. have identified linc-ADAL1, a non-conserved LncRNA, as a pivotal regulator of human adipogenesis. It facilitates this process through its interactions with hnRNPU and IGF2BP2 in distinct subcellular compartments [20] lncRNA ADINR has been established as a pivotal regulator of the CCAAT/enhancer-binding protein alpha (C/EBP$\alpha$) gene, which is a critical gene for the differentiation of preadipocytes. Through specific interaction with PA1, a component of the histone methylation complex MLL3/4, ADINR stimulates the activation of the MLL3/4 methyltransferase complex, leading to an augmentation of the trimethylation of H3K4me3 and ultimately facilitating the transcriptional activation of the C/EBP$\alpha$ gene [21] Jiang et al. demonstrated that cirfut10 acts as a sponge for let-7c, thereby promoting the proliferation of bovine adipocytes while inhibiting their differentiation through the targeting of PPARGC1B [22]. The aforementioned studies collectively demonstrate that ncRNAs are crucial in regulating the proliferation and differentiation of adipose tissue.

Whole transcriptome sequencing technology provided the possibility to elucidate the reasons for the differences in intramuscular fat deposition between Xinjiang Brown cattle and Angus $\times$ Wagyu cattle. This study utilized high-throughput sequencing to compare the expression levels of ncRNAs in the intermuscular fat tissue of Xinjiang Brown Cattle and Angus $\times$ Wagyu Cattle, and explored candidate ncRNAs and regulatory networks involved in adipogenesis. These findings provide novel insights into the potential mechanisms of ncRNAs mediated fat deposition in beef.

Result

Comparison of IMF content and morphological observation of the longissimus dorsi muscle in two cattle breeds

This study revealed a significant difference in intramuscular fat (IMF) content between the AW and XB groups within the longissimus dorsi muscle ($p < 0.0001$). The IMF content in the AW group was approximately 31.48%, markedly exceeding the 15.01% observed in the XB group (Fig. 1A). Muscle samples from both AW and XB were sectioned, as illustrated in the figure below, which presents cross-sectional views of their respective longissimus dorsi muscles. Notably, there was a pronounced disparity in IMF deposition, with a greater accumulation evident in the AW group compared to that of the XB group (Fig. 1B).

Fig. 1

IMF content and morphological observation of the longissimus dorsi muscle in XB and AW. A IMF content in the longissimus dorsi muscles of XB and AW. Results were expressed as mean ± standard deviation, n = 4. ****$p < 0.001$. B Histological Analysis of Longissimus Dorsi Muscle in XB and AW. Hematoxylin eosin staining was used, the nucleus was dyed to blue, muscle fibers were dyed to red, adipose tissue is white. Magnification: 100$\times$

Sequencing and profiling of miRNAs, lncRNAs and circRNAs

We performed miRNA sequencing on 8 muscle cDNA libraries, yielding a total of 109,343,831-117,258,577 high-quality reads. Through alignment with the reference genome sequence, we identified 673 known mature miRNAs and discovered an additional 121 novel mature miRNAs.

A total of 81,607,756-10,291,7174 clean reads were obtained from the cDNA libraries of 8 muscle samples for lncRNA and circRNA analysis. Subsequently, 2,199 known lncRNAs, 481 novel lncRNAs, and 8,648 newly identified circRNAs were characterized. The majority of circular RNAs exhibit a length predominantly within the range of 300 to 450 nucleotides (Fig. 2).

Fig. 2

The distribution of sequence lengths observed in the XB and AW libraries derived from circRNA sequencing results. Different colors represent the length distribution of circRNA in the longissimus dorsi muscle of various individuals

We compared the expression levels of three types of ncRNAs between the AW and XB groups. A Venn diagram illustrates that 361 miRNAs, 1,672 lncRNAs, and 7,480 circRNAs are expressed in the AW group, whereas 363 miRNAs, 1557 lncRNAs, and 5894 circRNAs are present in the XB group (Fig. 3A, B, C). Through the analysis of gene expression distribution within the samples, we observed that the expression levels of various ncRNAs across both groups were relatively uniform (Fig. 3D, E, F).

Fig. 3

The expression levels of three ncRNAs in XB and AW were analyzed. A miRNAs that are co-expressed in both XB and AW at the transcript level. B LncRNAs that are co-expressed in both XB and AW at the transcript level. C circRNAs that are co-expressed in both XB and AW at the transcript level. D The expression level distribution of miRNAs in XB and AW. E The expression level distribution of lncRNAs in XB and AW. F The expression level distribution of circRNAs in XB and AW

The consistency in gene expression among samples is a critical factor for assessing experimental stability and the appropriateness of sample selection. In this study, the R$^2$ value for miRNA surpassed 0.91, demonstrating a robust correlation among the samples (Fig. 4A). Concurrently, the R$^2$ values for lncRNA and circRNA exceeded 0.77, also indicating a significant degree of correlation between these samples (Fig. 4B, C).

Fig. 4

A correlation heatmap illustrating the relationships among various samples. A Heatmap depicting miRNA correlation across samples. B Heatmap depicting lncRNA correlation across samples. C Heatmap depicting circRNA correlation across samples

A comprehensive analysis of differential expression profiles of miRNAs, lncRNAs and circRNAs across various cattle breeds

In comparison to XB, a total of 53 miRNAs demonstrate differential expression in AW, with 38 exhibiting upregulation and 15 downregulation (Fig. 5A). Additionally, there are 176 differentially expressed lncRNAs in AW relative to XB, consisting of 80 that are upregulated and 96 that are downregulated (Fig. 6A). Furthermore, compared to XB, AW reveals the differential expression of 234 circRNAs, including 146 that are upregulated and 88 that are downregulated (Fig. 7A). The heatmaps illustrating the clustering of three ncRNAs demonstrate a pronounced clustering effect for genes that exhibit differential expression between the two groups (Figs. 5B, 6B and 7B).

Fig. 5

Differential expression miRNAs in the longissimus dorsi muscles between XB and AW. A Volcano plot illustrating the differential expression of miRNAs in XB and AW. The red dots indicate miRNAs that are upregulated ($p < 0.05$), while the green dots denote downregulated miRNAs ($p < 0.05$). The blue dots represent miRNAs exhibiting no significant expression differences. The dashed line illustrates the threshold for fold change in miRNAs expression. B Differentially expressed miRNAs between XB and AW

Fig. 6

Differential expression lncRNAs in the longissimus dorsi between XB and AW. A Volcano plot illustrating the differential expression of lncRNAs in XB and AW. The red dots indicate lncRNAs that are upregulated ($p < 0.05$), while the green dots denote downregulated lncRNAs ($p < 0.05$). The blue dots represent lncRNAs exhibiting no significant expression differences. The dashed line illustrates the threshold for fold change in gene expression. B Differentially expressed lncRNAs between XB and AW

Fig. 7

Differential expression circRNAs in the longissimus dorsi between XB and AW. A Volcano plot illustrating the differential expression of miRNAs in XB and AW. The red dots indicate circRNAs that are upregulated ($p < 0.05$), while the green dots denote downregulated circRNAs ($p < 0.05$). The blue dots represent circRNAs exhibiting no significant expression differences. The dashed line illustrates the threshold for fold change in circRNAs expression. B Differentially expressed circRNAs between XB and AW

Enrichment analysis of differentially expressed RNA and their target genes

In the diverse longissimus dorsi muscles of bovines, GO enrichment analysis was conducted on the target genes of distinct miRNAs to elucidate their biological functions across different species. The findings indicated that a substantial proportion of genes demonstrated significant enrichment ($p < 0.05$) in phosphorylation and protein phosphorylation (Fig. 8A).

Fig. 8

GO enrichment analysis of ncRNAs target genes. A GO enrichment of miRNA target genes. B GO enrichment of lncRNA target genes. C GO enrichment of circRNA target genes

After performing target gene prediction for lncRNAs, we employed expression correlation to predict their target genes. Subsequently, GO enrichment analysis was conducted on the target genes of differentially expressed lncRNAs. when padj < 0.05, a total of 17 GO terms were found to be significantly enriched and revealing that a majority of the genes were enriched in the G protein-coupled receptor signaling, regulation of cell death and apoptosis, transmembrane signaling receptor activity, molecular function regulator, molecular transducer activity, GTPase regulator activity, GTPase active activity and control of nucleotide triphosphatase activity (Fig. 8B).

For circRNAs, although no significant enrichment is observed at a padj < 0.05 threshold, GO terms such as G-protein coupled receptor signaling pathway, cytoskeleton organization, and intracellular organelle organization remain noteworthy (Fig. 8C).

The study conducted GO enrichment analysis on the prospective target genes and source genes of differentially expressed ncRNAs using three distinct methods to predict the relationships between ncRNAs target genes and their sources. The GO database annotated a total of 6283 differential miRNA target genes, 18148 differential lncRNA target genes, and 1736 circRNA source genes. In the AW vs XB comparison, ncRNAs were predominantly enriched in cellular processes, signaling pathways, cellular components, organelle processes, and catalytic activities. These findings suggest that the majority of differential ncRNAs may be associated with cellular tissue components and biological activities.

In the comparison of AW and XB, a significant increase in the number of target genes for differential miRNA, lncRNA, and circRNA as well as source genes was observed in 4, 79, and 1 KEGG pathways (padj < 0.05), respectively (Fig. 9). The mechanisms underlying the differential enrichment of miRNA target genes involve endocrine resistance, systemic lupus erythematosus, longevity regulating pathway - multiple species, viral carcinogenesis (Fig. 9A). The top 20 pathways associated with enrichment of lncRNA target genes include rheumatoid arthritis, NF-kappa B signaling pathway, intestinal immune network for lgA production, chemokine signaling pathway, Fc epsilon RI signaling pathway, Fc gamma R-mediated phagocytosis and several others (Fig. 9B). The processes contributing to the distinct enrichment of circRNA source genes encompass salivary secretion (Fig. 9C).

Fig. 9

KEGG enrichment analysis of ncRNAs target genes. A KEGG enrichment of miRNA target genes. B KEGG enrichment of lncRNA target genes C KEGG enrichment of circRNA target genes

Analysis of interactions among miRNAs, lncRNAs and mRNAs

We integrated the differentially expressed lncRNAs, miRNAs, and mRNAs along with their targeted regulatory interactions and the interaction network among ncRNAs to construct a comprehensive RNA regulatory network. The regulatory network comprising miRNAs, lncRNAs, and mRNAs encompasses a total of 212 nodes and 443 edges. This network includes 43 types of miRNAs, 24 types of lncRNAs, and 146 types of mRNAs (Fig. 10A). In the investigation of the regulatory interactions among lncRNA, miRNA, and mRNA, as well as the identification of potential regulatory factors, it was observed that upregulation of miRNA corresponded with downregulation of key nodes such as, AGBL carboxypeptidase 1 (AGBL1), and Thyroid hormone receptor beta (THRB) (Fig. 10B). Upon downregulation of miRNA expression, several genes exhibiting significant variations were identified, including S100 calcium-binding protein A13 (S100A13), CD74 molecule (CD74), lactose-binding glycoprotein-9 (LGALS9), and tensin 1 (TNS1) (Fig. 10C).

Fig. 10

The establishment and visualization of the LncRNA-miRNA-mRNA interaction network. A A network diagram illustrating the interactions among the target lncRNAs, miRNAs, and differentially expressed target genes across various combinations. B Assess the expression levels of downregulated lncRNAs, upregulated miRNAs, and their corresponding target genes. C Assess the expression levels of upregulated lncRNAs, downregulated miRNAs, and their corresponding target genes

Revalidation of RNA-Seq data

For qRT-PCR validation, we selected three of the most significantly upregulated and three of the most significantly downregulated miRNAs, LncRNAs, and circRNAs, respectively. The results obtained were in concordance with the sequencing analysis, providing additional evidence for the reliability of our sequencing data (Fig. 11).

Fig. 11

RT-qPCR verification of RNA-Seq results. We employed “log2(fold change)” as the ordinate, which signifies the fold change of the experimental group (AW) compare to the control group (XB), transformed into a logarithmic value with a base of 2

Discussion

Maximize the accumulation of intermuscular fat is crucial for improving beef quality. A comprehensive understanding of the mechanisms and regulatory pathways governing IMF storage is essential for advancing breeding practices and production techniques in beef cattle. As previously noted, ncRNAs are crucial regulators of adipogenesis and lipid accumulation at the post-transcriptional level; however, the precise mechanisms through which they exert their influence on these processes remain inadequately elucidated. Gao et al. reported that the IMF content of XB ranged from 5.8% to 7.6% between the ages of 28 and 34 months [23]. In contrast, the IMF content in Wagyu cattle consistently remained at or above 30% [24], which aligns with our research findings. In this study, we assessed the IMF levels in both AW and XB using the Soxhlet extraction method. The findings revealed that the IMF content in the AW was approximately 31.48%, significantly surpassing the 15.01% observed in the XB. Subsequently, we conducted slicing experiments, which demonstrated a markedly greater intermuscular fat in the AW compared to that of the XB.

Concurrently, we performed high-throughput sequencing of the complete transcriptome to identify differentially expressed ncRNAs in the IMF of Xinjiang Brown cattle and Angus $\times$ Wagyu cattle. Compared to the XB group, the longissimus dorsi muscle in the AW group we identified 80 upregulated miRNA and 15 downregulated miRNAs. Furthermore, a total of 80 lncRNAs exhibited increased expression levels, while 96 showed decreased expression. Additionally, among the analyzed circRNAs, 146 were found to be upregulated and 88 demonstrated downregulation. The findings reveal that the transcriptomes of Xinjiang Brown cattle exhibit significant differences compared to those of Angus $\times$ Wagyu cattle, implying that regulatory factors may play a role in the accumulation of IMF.

Examine the molecular pathways through which ncRNAs may modulate adipose tissue accumulation in cattle. Notably, bta-miR-142-3p ranks among the top 10 most abundant miRNAs and exhibits differential expression levels between Xinjiang Brown cattle and Angus $\times$ Wagyu cattle. The study demonstrated that the stem-loop structure of miR-142-3p directly modulates the translation of beta-catenin by interacting with the 3’ untranslated region of Ctnnb1 mRNA, thereby inhibiting the Wnt/beta-catenin signaling pathway [25]. The classical Wnt/$\beta$-catenin signaling pathway serves as an intrinsic regulator that modulates the fate of mesenchymal cells, facilitating osteoblastogenesis while concurrently suppressing adipogenic differentiation [26]. Research indicates that elevated levels of miR-142-3p can result in a reduction of MAPK mRNA expression. Additionally, JNK, which is modulated by MAPK, phosphorylates serine residues on insulin receptor substrate 1 (IRS1). This phosphorylation disrupts insulin signaling and impedes the adipogenic differentiation of bone marrow mesenchymal stem cells (BMSCs), thereby acting as a negative regulator in the maturation of adipocytes [27]. Our findings indicate that miR-142-3p is may a crucial regulator of lipid synthesis and subsequent adipose deposition in Xinjiang Brown cattle as well as Angus $\times$ Wagyu cattle.

Additionally, we observed that the expression levels of bta-miR-147 and bta-miR-217 in Angus $\times$ Wagyu cattle were significantly elevated compared to those in Xinjiang Brown cattle, whereas the expression levels of bta-miR-432, bta-miR-1306, and bta-miR-211 were markedly reduced. Furthermore, there is a paucity of literature regarding their involvement in adipogenesis and fat deposition. The precise role of these miRNAs in adipose tissue accumulation remains ambiguous; thus, we will concentrate on this aspect and conduct further investigations.

Research has demonstrated that lncRNAs play significant roles in the metabolism and regulation of lipids. LncRNAs have been shown to enhance cellular proliferation [28], while numerous target genes associated with lncRNAs are implicated in the process of lipid deposition in mammal [29, 30]. Recent studies have demonstrated that the long ncRNAs, designated as IMFlnc1, is implicated in the deposition of IMF. It functions as a competing endogenous RNA (ceRNA) by interacting with miR-199a-5p, which modulates CAV-1 expression and subsequently influences fat accumulation within porcine muscle tissue [31]. In a related investigation, Cai et al. identified a novel long ncRNAs designated as BADLNCR1, which impedes the differentiation of bovine adipocytes by downregulating GLRX5 expression, a gene that plays a crucial role in lipogenesis and lipid synthesis [32]. Although it has been demonstrated that lncRNA is involved in the regulation of adipose tissue accumulation in mammals, the comprehensive regulatory framework governing IMF deposition remains inadequately understood. In our research, compared to the XB group, the unnamed lncRNAs in the AW group we observed significantly higher expression levels such as TCONS_00060413, TCONS_00047393, and TCONS_00034659 in AW compared to XB (p < 0.05), whereas the expression levels of TCON_00047449, TCON_00024142, and TCON_00044387 were relatively lower. To the best of our knowledge, this is the first report demonstrating the involvement of these lncRNAs in the regulation of fat deposition. Their roles in adipogenesis remain unclear and warrant further investigation in future studies, This study also provides a direction for our future research.

CircRNA plays a crucial role in the regulation of lipid metabolism. CircRNA functions is a non-linear RNA, serving as a protein sponge and miRNA sponge, while regulating the activity of target genes through post-transcriptional suppression via 3’-UTR pairing. mmu_circ_0001874 modulates lipid metabolism by targeting the miR-24-3p/Igf2/PI3K-AKTmTOR signaling axis in mice [33]. Prior research has demonstrated that the CIRFUT10/let-7c/PPARGC1B signaling pathway plays a pivotal role in the proliferation and differentiation of bovine adipocytes [22]. In this study, a considerable number of novel circRNAs exhibited differential expression between AW and XB. It is of particular significance to observe that the source gene LVRN of novel_circ_0000394 has been suggested to exert a potential influence on the proliferation and differentiation of adipocytes [34]; Meanwhile, Nishizuka et al. have demonstrated that the source gene KCNMA1 of novel_circ_0008017 is involved in the regulation of insulin signaling in mature adipocytes of mice. The functionality of LVRN and KCNMA1 may be attributed to the transcription of novel_circ_0000394 and novel_circ_0008017, which subsequently interact with relevant target miRNAs. this may constitute one of the reasons for the observed differences in novel_circ_0000394 and novel_circ_0008017 expression. This will be one of the key focuses of our future research endeavors. Based on the latest findings in the field of research, further investigations are essential to elucidate the precise mechanisms through which circRNAs influence adipocyte growth and differentiation.

The analysis of GO annotations indicates that a substantial proportion of ncRNAs target genes are predominantly engaged in cellular components and molecular functions. It is suggested that that the majority of differentially expressed ncRNAs may be intricately linked to the components of cellular tissues and their associated biological roles, specifically, regarding the relevant functions of cellular membrane receptors and transport channels. We have observed that numerous target genes of lncRNAs converge on terms related to enzyme activators and modulators, particularly those associated with GTPases. Rho GTPases constitute a crucial class of cellular signaling proteins that play pivotal regulatory roles within cells. Rho family GTPase 3 (RND3), a constituent within the protein superfamily, functions to inhibit Rho-associated kinase (ROCK) signaling and has been implicated in diverse pathological conditions, such as apoptotic cardiomyopathy, heart failure, cancer, and type 2 diabetes mellitus. This involvement is partially mediated through its modulation of cytoskeletal dynamics and insulin-facilitated glucose uptake processes. The study by Dankel et al. unveiled differential expression of RND3 in obesity and insulin resistance. By knocking down RND3, they observed decreased protein levels of ATGL and Ser650-phosphorylated HSL, which attenuated the stimulatory effects of isoproterenol and cAMP signaling on lipolysis; Baba et al. have demonstrated that Rac GTPase influences adipogenesis by participating in the regulation of the expression of receptor activator of nuclear factor-$\kappa$B ligand (RANKL). These studies collectively affirm the involvement of GTPases in the process of adipose deposition, leading us to hypothesize that enzyme activation, exemplified by GTPases, plays a significant role in this process, and the target genes of the differentially expressed ncRNAs identified in this study are implicated in these processes, suggesting a potential avenue for future research into the mechanisms of adipose deposition.

Concurrently, we have also identified the convergence of numerous cell membrane-related terms, such as membrane composition and signaling receptor pathway activity, highlighting the significant roles played by signaling receptors on the cell membrane in organisms. The TMEM proteins are a class of transmembrane protein families that contribute to the composition of cellular membranes and are widely acknowledged to play pivotal roles in numerous critical biological processes. TMEM182 is highly expressed in muscle and adipose tissues, Chen et al. found that overexpression of TMEM182 enhanced the expression of genes related to adipogenesis and promoted the differentiation of preadipocytes into adipocytes [35]; Wang et al. demonstrated that TMEM68 functions as an acyltransferase and influences the expression of adipogenesis-related genes, glycolipid metabolism, and triacylglycerol storage [36]. The aforementioned studies indicate that ncRNAs may interact with their target genes, thereby participating in and influencing the composition of cell membrane components, ultimately affecting fat deposition. Notably, G protein-coupled receptors represent the most abundant class of cell surface receptors, capable of binding to chemical ligands in the cellular microenvironment and initiating a cascade of intracellular signaling pathways, thereby exerting significant biological significance. Numerous research data indicate that G-protein-coupled receptors play a regulatory role in glucose metabolism and lipid metabolism. Previous studies have demonstrated that the gastric inhibitory polypeptide receptor (GIPR), characterized by seven transmembrane domains, activates G proteins to elevate intracellular levels of cAMP and $Ca^{2+}$. Furthermore, it initiates signaling cascades involving PI3K, PAK, and PKB to modulate the expression of downstream genes [37]. In addition, GIPR facilitates insulin secretion, resulting in a reduction of blood glucose levels. Additionally, it is involved in the metabolism and storage of adipose tissue and influences the production and release of glucose within hepatic tissues [16]. These studies collectively suggest that signaling receptor activity may play a crucial role in ncRNA-mediated biological processes during adipogenesis. Our research underscores the profound significance of these pathways and functions, which will constitute a vital research direction for us in the future. We will continue to focus on elucidating the specific mechanisms involved.

Research has demonstrated that protein phosphorylation is crucial for the reorganization of the lipid metabolism network [38]. The previously discussed evidence indicates that protein phosphorylation is implicated in the regulation of lipid metabolism. In this study, the target genes exhibit significant enrichment in KEGG pathways associated with phosphorylation, indicating a potential link between phosphorylation modifications and their influence on enzyme activity, as well as the subsequent alterations in lipid metabolism process they regulate. however, the precise mechanisms underlying this regulatory process necessitate further investigation. NF-KappaB prevents cellular apoptosis by downregulating adenine nucleotide translocase 2 (ANT2) [39]; Studies have indicated that the interaction between the chemokine system and adipocytes is involved in the regulation of systemic energy balance and glucose homeostasis [40]; Eguchi et al. have demonstrated that an adhesion molecule in adipocytes is specifically upregulated in mature adipocytes and adipose tissue of obese mice and humans, and they speculate that it may play a role in the process of adipocyte maturation [41]; Rodríguez et al. have indicated that adipocytokine signaling plays a certain role in energy metabolism [42]; In our study, a significant number of target genes were enriched in these pathways, which may provide us with a more comprehensive understanding for future exploration of the process of adipocyte deposition and will constitute a novel domain of inquiry for our research. Unfortunately, we have yet to uncover any relationship between pathways not mentioned, such as the Fc epsilon RI signaling pathway, and fat deposition. This may represent a potential research direction for future studies.

Furthermore, the interactome network encompassing miRNAs, lncRNAs, and mRNAs reveals that numerous mRNAs and lncRNAs are situated at critical junctions. Notably, when miRNA expression is upregulated, both AGBL1 and THRB exhibit downregulation. AGBL1 is correlated with circulating adiponectin levels, which may influence insulin sensitivity [43]. THRB is a nuclear receptor encoded by the THRB gene, exhibiting a high affinity for triiodothyronine (T3). Its mechanism of action primarily involves the regulation of target gene transcription through T3 binding [44]. Moreover, the expression and activity of THRB are intricately linked to lipid metabolism in the liver. In mice deficient in THRB, a marked alteration in the circadian metabolic rhythm of hepatic function is observed, resulting in elevated levels of cholesterol, triglycerides, diacylglycerols, and fatty acids. This dysregulation may contribute significantly to the pathogenesis of fatty liver disease [45]. Concurrent research indicates that the expression level of THRB is markedly diminished in obese individuals relative to those of normal weight, implying that THRB may serve as a crucial regulatory factor [46].

When the levels of miRNAs are downregulated, certain RNAs situated within the central nodes of the interaction network may contribute to the generation and accumulation of adipocytes. Our findings reveal that S100A13 is significantly upregulated within the interaction network. The S100A13 gene encodes the S100A13 protein, which independently binds to calcium and copper ions; this interaction is essential for the expression of IL1A and FGF1 in response to copper-dependent stress [47]. Interleukin 1 alpha (IL1A) is one of the three isoforms of interleukin-1 and functions as a pro-inflammatory cytokine produced by macrophages. In adipocytes from diet-induced and genetically obese mice, the expression level of IL1A is significantly elevated [48], a finding that has also been observed in the subcutaneous adipose tissue of obese patients [49, 50]. Research has demonstrated that the non-classical IL-1 - MyD88 - IRAK2 - PHB/OPA1 signaling pathway plays a crucial role in regulating mitochondrial metabolism within adipocytes, potentially exacerbating obesity [51]. The aforementioned research findings align with our own. Consequently, the differentially expressed ncRNAs identified in this study, along with the established regulatory network involving miRNAs, lncRNAs, and mRNAs, are poised to significantly influence adipocyte growth and differentiation. This may subsequently impact IMF accumulation. Despite certain limitations inherent to this research, it is imperative to further investigate these potential regulatory pathways.

Conclusion

In conclusion, we employed high-throughput whole transcriptome sequencing to investigate the levels of ncRNAs in the longest muscle of Xinjiang Brown cattle compared to Angus $\times$ Wagyu cattle in this study. The aim was to elucidate potential ncRNAs implicated in fat deposition. Our findings identified several ncRNAs, such as miR-142-3p, novel_circ_0000394 and nove_circ_0008017 that may modulate lipogenesis and fat accumulation through their target genes, source genes, and associated pathways. Furthermore, we established a regulatory network comprising miRNAs, lncRNAs, and mRNAs. We identified that proteins associated with insulin sensitivity, hepatic lipid metabolism, and interactions with calcium and copper ions constitute pivotal components of the network. These results deepen our understanding of the differential expression of ncRNAs within intermuscular fat across two distinct cattle breeds, thereby laying a foundation for further investigation into the molecular mechanisms governing beef fat deposition. This research offers significant references for identifying key regulatory ncRNAs influencing adipogenesis in cattle.

Material and methods

Experimental animals and sample preparation

In this experiment, we selected four 26-month-old XB and AW separately that had been fattened. ensuring that all the four cattle within the same group share an identical genetic background and were raised under identical environmental conditions. The feed formulation adheres to the guidelines established by the National Research Council (NRC) of the United States, incorporating whole-plant corn silage, corn stover, corn grains, wheat stover, alfalfa, cottonseed meal, soybean meal, and wheat bran. This research received approval from the Animal Ethics Committee at Xinjiang Agricultural University.

All experimental animals were subjected to humane slaughter at the designated facility of Shengyuan Cattle and Sheep in Changji, Xinjiang. Following the slaughtering process, samples from between the 12th and 13th ribs of the longissimus dorsi muscle were promptly collected using surgical scissors, chopped into smaller pieces. All muscle specimens were divided into two groups: one group was processed into paraffin sections with 4% formalin for morphological analysis, while the other group was preserved in liquid nitrogen for transcriptomic assessment and differential ncRNAs expression analysis.

Determination of IMF content

The IMF content was quantified using the Soxhlet extraction method. A 3-gram sample of longissimus dorsi muscle was homogenized and placed in a dry paper bag for weighing. The sample underwent continuous drying at 65°C in an oven for no less than 15 hours, after which it was transferred to a desiccator to cool for 10 minutes before reweighing. Subsequently, the bags were inserted into a Soxhlet extractor and immersed in ether overnight. Following this step, they were refluxed in a water bath maintained at 75°C for a duration of 10 hours, then allowed to evaporate naturally in a well-ventilated area until all ether had completely dissipated. Finally, the samples were dried in an oven set at 105°C until reaching constant weight, followed by another weighing to determine the IMF content. IMF% = (Mass of the Ether-Soluble Fraction / 3) $\times$ 100%.

Histology

In this study, we selected four samples of the longest dorsal muscle from both the AW and XB groups, which were subsequently sectioned into tissue blocks measuring 1 cm $\times$ 1 cm $\times$ 0.5 cm in accordance with the orientation of the muscle fibers. These blocks were immersed in a 4% formaldehyde solution (Wuhan Kanos Technology Co., Ltd., Wuhan, China) for a duration of 24 hours to ensure thorough fixation. Following fixation, dehydration was conducted by immersing them in a series of graded alcohol solutions (Xinjiang Aidil Biotechnology Co., Ltd, Xinjiang, China). The samples were then transferred to xylene (Heng Xing Reagent, Tianjing, China) for transparency treatment. After this process, they were embedded in paraffin wax (Shanghai Yiyang Instrument Co., Ltd., Shanghai, China), allowed to cool naturally until solidified before being sliced into sections. The paraffin sections underwent dewaxing with xylene Heng Xing Reagent, Tianjing, China) on microscope slides and were stained using hematoxylin (Zhuhai Beso Biotechnology Co., Ltd., Zhuhai, China), followed by gentle rinsing with water. Differentiation was performed using a hydrochloric acid solution at 0.1%, prepared as a mixture of concentrated hydrochloric acid and anhydrous ethanol diluted in ethanol prior to another wash with water. Eosin staining was subsequently applied (Zhuhai Beso Biotechnology Co., Ltd., Zhuhai, China), followed by washing with alcohol and completing dehydration steps. Finally, re-transparency was achieved through xylene before sealing the sections with resin (Shanghai Yiyang Instrument Co., Ltd., Shanghai, China) for preservation during analysis. Upon completion of preparation procedures, morphological characteristics of the tissue sections were examined utilizing a FEI Quanta 250 field emission scanning electron microscope.

Total RNA library construction and sequencing

Total RNAs were extracted from the muscular tissue using TRIzol (Invitrogen, Carlsbad, CA). For miRNA, Total RNA was used as input material for the RNA sample preparations. Sequencing libraries were generated using NEB Next® Multiplex Small RNA Library Prep Set for Illumina® (New England Biolabs, Ipswich, UK). Briefly, 3’ and 5’ adaptors were ligated to 3’ and 5’ end of small RNA, respectively. Subsequently, the first strand cDNA was synthesized following hybridization with a reverse transcription primer. The double-stranded cDNA library was then generated through PCR enrichment. After purification and size selection, libraries containing insertions ranging from 18 to 40 bp were prepared for sequencing using the Illumina SE50 platform. Subsequently, library quality was assessed on the Agilent 5400 system (Agilent, Santa Clara, CA) and quantified by QPCR (1.5 nM). The qualified libraries were combined and sequenced using the SE50 strategy on Illumina platforms at Novogene Bioinformatics Technology Co., Ltd (Beijing, China), in accordance with the necessary library concentration and data requirements.

For LncRNA analysis, 1$\mu$g of total RNA was utilized as input material for the RNA sample preparations. Sequencing libraries were constructed using the NEBNext® UltraTM RNA Library Prep Kit for Illumina® (New England Biolabs, Ipswich, UK) according to the manufacturer’s recommendations, and index codes were incorporated to assign sequences to each sample. In brief, Regulatory ncRNA and mRNA were isolated from total RNA using specific probes to deplete rRNA. The fragmentation process was conducted with divalent cations at an elevated temperature in the First Strand Synthesis Reaction Buffer (5X). First strand cDNA synthesis was performed using a random hexamer primer and M-MuLV Reverse Transcriptase (RNaseH). The second strand cDNA synthesis was subsequently conducted using DNA Polymerase I and RNase H. The remaining overhangs were then converted into blunt ends through the activities of exonuclease and polymerase. After adenylation of 3’ ends of DNA fragments, NEBNext Adaptor with a hairpin loop structure was ligated to facilitate hybridization. To selectively isolate cDNA fragments ranging from 370 to 420 bp in length, the library fragments were purified using the AMPure XP system (Beckman Coulter, Beverly, USA). The size-selected, adaptor-ligated cDNA was treated with 3 $\mu$L of USER Enzyme at 37°C for 15 min, followed by a denaturation step at 95°C for 5 min prior to PCR amplification. The PCR reaction was performed using Phusion High-Fidelity DNA polymerase, Universal PCR primers, and Index Primer. The PCR products were finally purified to library quality, which was assessed using the Agilent 5400 system and quantified by QPCR (at a concentration of 1.5 nM). Subsequently, the qualified libraries were pooled together and sequenced on Illumina platforms using the PE150 strategy at Novogene Bioinformatics Technology Co., Ltd (Beijing, China), based on the required effective library concentration and data amount.

The input material for RNA sample preparations in circRNA analysis consisted of a total amount of 2 $\mu$g RNA per sample. Following the manufacturer’s recommendations, sequencing libraries were generated using the NEB RNA Library Prep Kit for Illumina (USA), and index codes were incorporated to assign sequences to each individual sample. Firstly, ribosomal RNA was eliminated and the rRNA-free residue was purified through ethanol precipitation. Subsequently, the linear RNA was enzymatically degraded using 3 U of RNase R per 1 $\mu$g of RNA. Fragmentation was performed by exposing the RNA to divalent cations at an elevated temperature in First Strand Synthesis Reaction Buffer (5X). The first strand cDNA was synthesized utilizing a random hexamer primer and M-MuLV Reverse Transcriptase (RNaseH-). Subsequently, the second strand cDNA synthesis was performed employing DNA Polymerase I and RNase H. The remaining overhangs were converted into blunt ends through the activities of exonucleases and polymerases. Following adenylation of the DNA fragments’ 3’ ends, NEBNext Adaptor with a hairpin loop structure was ligated for hybridization preparation. To selectively isolate cDNA fragments ranging from 370 to 420 bp in length, the library fragments were purified using the AMPure XP system (Beckman Coulter, Beverly, USA). The size-selected, adaptor-ligated cDNA was incubated with 3 $\mu$L of USER Enzyme (NEB, USA) at 37°C for 15 minutes, followed by a denaturation step at 95°C for 5 minutes prior to PCR amplification. The subsequent PCR reaction utilized Phusion High-Fidelity DNA polymerase, Universal PCR primers, and Index (X) Primer. Finally, the products were purified and the library quality was evaluated using the Agilent 5400 system, followed by quantification through qPCR (1.5 nM). The qualified libraries were then pooled and subjected to sequencing on Illumina platforms using a PE150 strategy at Novogene Bioinformatics Technology Co., Ltd (Beijing, China), based on the required effective library concentration and data volume.

Real-time quantitative PCR (RT-qPCR)

RNA was extracted from 8 muscle samples for RT-qPCR analysis, with three biological replicates per sample. The Trizol method was utilized for total RNA extraction. miRNA was converted into cDNA using the miRNA 1st Strand cDNA Synthesis Kit (by stem-loop) (Nanjing Novogen Biotech Co., Ltd., Nanjing, China), while lncRNA and circRNA were synthesized into cDNA employing the PrimeScript^TM RT Reagent Kit (Perfect Real Time) (Dalian Boh Bioengineering Co., Ltd., Dalian, China). The analysis was performed utilizing the Bio-Rad CFX96 instrument, specifically engineered for fluorescence quantitative PCR. In this study, we randomly selected 18 ncRNAs from the RNA extracted from muscle tissues of Xinjiang Brown cattle and Angus $\times$ Wagyu cattle. The extraction process employed the miRNA Universal SYBR qPCR Master Mix Fluorescence Quantitative PCR Reagent (Nanjing Novogen Biotech Co., Ltd., Nanjing, China), while analyses for lncRNA and circRNA were conducted using the TB Green® Premix Ex Taq^TM (Tli RNaseH Plus) Quantitative PCR Kit (Dao Bioengineering Co., Ltd., Dalian, China). The sequences of the primers were meticulously designed Primer 5 software, with U6 and GAPDH serving as internal control genes. 95°C for 10 seconds and 60°C for 30 seconds (40 cycles), while the RT-qPCR amplification protocol for lncRNA and circRNA was set as 95°C for 30 seconds, 95°C for 5 seconds and 60°C for 30 seconds (40 cycles). Primer sequences were designed using Primer 5 (Table 1) and the internal reference gene was GAPDH and U6. Result analysis was conducted using the 2$^{-\Delta \Delta CT}$ method.

Table 1 The primer of differentially expressed ncRNAs for RT-qPCR

Data processing

Quality control

The raw data (raw reads) in fastq format were initially processed using custom perl scripts. During this step, clean data (clean reads) were obtained by removing reads containing adapters, ploy-N sequences, or low-quality scores from the raw data. Simultaneously, the Q20, Q30, and GC content of the clean data were calculated. All subsequent analyses were performed exclusively on the high-quality clean data.

Bioinformatics analysis of small RNA sequences

The reference genome and gene model annotation files were obtained from the genome website (https://ftp.ensembl.org/pub/release-76/fasta/bos_taurus/dna/). The reference genome was indexed using Hisat2 v2.0.5, and the paired-end clean reads were aligned to it using Hisat2 v2.0.5. For miRNA analysis, the small RNA tags were mapped to the reference sequence without allowing any mismatches using Bowtie to assess their expression levels and distribution.

The raw counts were initially normalized using the TPM method proposed by Zhou et al. Subsequently, differential expression analysis was conducted for miRNA and circRNA between two conditions/groups, with each condition having two biological replicates per group. This analysis was performed utilizing the DESeq2 R package (version 1.20.0).

The expression levels of LncRNA and circRNA were quantified using StringTie (version 1.3.3b) by calculating FPKM values. Differential expression analysis between two samples was conducted using the edgeR R package (version 3.22.5). Prior to differential gene expression analysis, read counts for each sequenced library were normalized by the edgeR program package through a scaling normalization factor, both with and without biological replicates. All resulting P-values for RNAs were adjusted using Benjamini and Hochberg’s method to control the false discovery rate. The differential expression of mRNA, circRNA, and lncRNA was determined based on the following criteria: |log2(Fold Change)$|$ $\ge$1 and adjusted P-value <0.05.

GO and KEGG enrichment analyses

The differential expression analysis of genes was performed using the clusterProfiler (3.8.1) R package, which incorporated gene length bias correction for accurate Gene Ontology (GO) enrichment analysis. GO terms with a corrected P-value below 0.05 were considered significantly enriched by the differentially expressed genes. The enrichment analysis is based on the principle of hypergeometric distribution, where the differential gene set refers to the set of genes obtained from significant difference analysis and annotated to either the GO or KEGG database. The background gene set includes all genes analyzed for significant difference and annotated to either the GO or KEGG database. The software packages ClusterProfiler (v3.10.1) and miRanda (3.3a) were employed for conducting Gene Ontology (GO) and Kyoko Encyclopedia of Genes and Genomes (KEGG) pathway analyses on the differentially expressed circRNA, lncRNA target genes, and miRNA target genes. The GO enrichment analysis encompasses biological processes (BP), cellular components (CC), and molecular functions (MF). Significantly enriched GO terms or KEGG pathways were identified based on a corrected p-value threshold of $q<0.05$.

Statistical analysis

Data were organized using Microsoft Office Excel (version 2020, Microsoft Corporation, Washington), with IMF and RT-qPCR results expressed as means ± SD. Independent t-tests were conducted via SPSS (version 20.0, IBM, New York), considering $P<0.05$ as statistically significant for differential ncRNAs. Correlation heatmaps (ggplot2), gene expression box plots (ggplot2 pack-age), clustering heatmaps (ggplot2 package), volcano plots (ggplot2 package), and Venn Diagrams of co-expression (Venn Diagram package) were generated using R (Version 4.3.2). Additionally, IMF content visualizations and differential ncRNA quantification plots were produced utilizing GraphPad Prism (version 9.3, GraphPad Software, San Diego, CA).