Multi-Omics Integration Identifies Key Pathways and Regulatory Genes Driving Marbling Formation and Meat Quality in Yunling Cattle

TL;DR

This study identifies the molecular mechanisms behind marbling in Yunling cattle, revealing key genes and pathways that contribute to high-quality beef.

Contribution

The study provides a multi-omics framework to uncover novel regulatory genes and metabolic pathways specific to Yunling cattle marbling.

Findings

Yunling cattle show distinct lipid and amino acid profiles compared to Angus and Simmental breeds.

Key genes like NDUFB1, COX7C, and ITGB3 are linked to marbling formation through co-expression network analysis.

Energy metabolism and cell communication pathways are uniquely active in Yunling cattle, supporting fat deposition.

Abstract

Marbling, the fat found inside muscle, is a key factor in determining beef quality, making the meat tender and flavorful. Yunling cattle, a breed from China, are known for their potential to produce excellent marbling, but the biological reasons behind this trait are not fully understood. This study investigated the muscle characteristics of Yunling cattle by comparing them with Angus and Simmental breeds. By analyzing gene activity, fat composition, and amino acid levels, we identified unique molecular patterns in Yunling cattle. We found that specific groups of genes work together with distinct fat molecules and amino acids to promote marbling formation. Specifically, Yunling cattle exhibit a unique energy metabolism and cell communication system that supports fat deposition. These results provide a scientific explanation for the superior meat quality of Yunling cattle. This knowledge is valuable for the beef industry as it helps breeders select the best animals, leading to genetic improvements and ensuring a consistent supply of high-quality, tasty beef for consumers.

Marbling, or intramuscular fat (IMF), is a primary determinant of high-quality beef, defining key sensory attributes and nutritional value. Yunling (YL) cattle, an indigenous breed from Yunnan, China, are renowned for their superior marbling, yet the underlying molecular mechanisms remain unclear. This study employed an integrated transcriptomic, lipidomic, and amino acid metabolomic approach to systematically compare the multi-omics profiles of the longissimus dorsi muscle among YL, Angus (AGS), and Simmental (XMTE) cattle. Transcriptome analysis identified 2053 and 2156 differentially expressed genes (DEGs) in XMTE vs. YL and AGS vs. YL, respectively. These DEGs were primarily enriched in the PI3K-Akt and MAPK signaling pathways, as well as oxidative phosphorylation. Lipidomic analysis revealed a distinct lipid profile in YL cattle, identifying 27 characteristic lipid molecules (e.g., SM(d20:0/24:1), DG(16:0/18:1(11Z)/0:0)) compared to XMTE and 17 differential lipids compared to AGS. The amino acid metabolome showed that Beta-Alanine and L-Aspartic acid levels in YL were 42.6% and 54.8% lower than in XMTE, respectively (p < 0.01), and levels of several functional amino acids were significantly reduced compared to AGS. Weighted Gene Co-expression Network Analysis (WGCNA) constructed a gene-metabolite network, identifying key modules strongly correlated with lipid and amino acid metabolism (|r| > 0.6). Within these modules, energy metabolism-related genes such as NDUFB1, COX7C, and IDH3B, along with signal transduction genes including ITGB3, PDGFRA, and FN1, were found to synergistically regulate marbling formation in YL cattle. This study systematically elucidates the molecular mechanisms underlying both marbling formation and the nutritional characteristics of meat in Yunling cattle. This provides a theoretical foundation for genetic improvement and offers potential molecular targets to enhance both marbling and overall meat quality in other indigenous cattle breeds worldwide.

Figures

Full text

1. Introduction

The global beef industry places a premium on meat quality, with intramuscular fat (IMF) content and its distribution, known as marbling, being a primary determinant of consumer preference and market value [1]. High-quality beef breeds such as Angus (AGS) and Wagyu are globally recognized for their superior marbling, which imparts exceptional tenderness, juiciness, and flavor [2,3,4]. Consequently, understanding and enhancing the molecular mechanisms of marbling formation is a central goal for the genetic improvement of beef cattle and the overall competitiveness of the beef market [5,6].

Beyond established international breeds, there is a growing interest in exploring unique local genetic resources that may harbor novel mechanisms for meat quality improvement. Yunling (YL) cattle, an indigenous breed from the Yunnan Plateau of China, represent such a resource. They are recognized for their potential to produce highly marbled beef, a trait developed through long-term adaptation to their specific ecological environment [7]. However, a significant challenge for the industrialization of YL cattle is the inconsistency of marbling expression; while some individuals exhibit superior IMF deposition, the trait is not uniformly stable across the population. This variability limits the scientific evaluation and targeted breeding efforts needed to elevate their market position. Therefore, a systematic comparison with benchmark breeds like the well-marbled Angus and the lean, fast-growing Simmental (XMTE) is crucial to decipher the unique molecular signature underlying the high-quality marbling potential in YL cattle and to identify pathways for its stabilization and enhancement.

Intramuscular fat deposition is a complex biological process involving the synergistic action of multiple factors, including preadipocyte differentiation, lipid synthesis and degradation, energy metabolism regulation, and cell signaling. Previous studies have shown that the expression levels of lipid metabolism-related genes, such as FABP4, SCD, and PPARG, are closely associated with beef marbling. Concurrently, energy metabolic pathways like oxidative phosphorylation and the TCA cycle also play crucial roles in IMF deposition. Furthermore, the amino acid metabolic network significantly regulates muscle development and fat deposition; for instance, L-Arginine influences muscle growth and fat distribution via the nitric oxide signaling pathway. However, most existing studies focus on a single omics level, lacking a multi-layered, systematic analysis. In particular, comprehensive multi-omics research on local breeds like YL is insufficient, hindering a full understanding of their marbling formation mechanisms.

The recent development of integrated multi-omics analysis has provided a new paradigm for dissecting complex biological traits [8]. Transcriptomics can reveal gene expression regulatory networks, lipidomics directly reflects the composition and metabolic state of lipid molecules, and amino acid metabolomics offers deep insights into the metabolic patterns of muscle tissue. Integrating these multi-omics datasets allows for a comprehensive elucidation of the molecular mechanisms of marbling from multiple perspectives. Notably, Weighted Gene Co-expression Network Analysis (WGCNA), a cutting-edge systems biology method, can identify gene modules associated with specific phenotypes and construct gene-metabolite correlation networks, showing great potential in studies of complex traits in agricultural animals [9]. Thus, a multi-omics integration strategy provides a technically feasible approach for an in-depth investigation of the mechanisms underlying marbling in YL.

This study employed an integrated multi-omics approach, with YL as the primary subject and Angus and Simmental cattle as controls, to systematically compare the transcriptomic, lipidomic, and amino acid metabolomic profiles of the longissimus dorsi muscle. The objective was to unravel the molecular regulatory network and mechanisms underlying marbling formation in YL. The specific aims were: (1) to identify key differentially expressed genes, characteristic lipid molecules, and differential amino acids that distinguish YL from the other two breeds; (2) to construct a gene-metabolite correlation network driving marbling formation and influencing meat quality in YL through multi-omics integration; and (3) to resolve the core signaling pathways and key regulatory nodes involved.

2. Materials and Methods

2.1. Experimental Animals and Sample Collection

All animal procedures in this study were reviewed and approved by the Life Science Ethics Committee of Yunnan Agricultural University (Ethics Approval No. 202406002). Longissimus dorsi muscle samples were collected from three cattle breeds: Yunling (YL, n = 3), Angus (AGS, n = 3), and Simmental (XMTE, n = 3). All experimental bulls were provided by Yunnan Tengchong Hengyi Dongshan Agricultural Development Co., Ltd. They were maintained under identical fattening conditions and dietary nutrient levels, and were 24 months of age at the start of the trial. At the end of the fattening period, the bulls were fasted for 24 h (with free access to water) before slaughter. Within 30 min post-slaughter, left-side samples of the Musculus longissimus dorsi were precisely collected at the 12th–13th thoracic rib interface. To ensure sample homogeneity, each sample was taken from the central part of the same anatomical location. The collected muscle tissue (approx. 5 g) was trimmed of visible fat and connective tissue, immediately snap-frozen in liquid nitrogen, and subsequently stored at −80 °C for transcriptomic, lipidomic, and amino acid metabolomic analyses (Figure 1). Two experimental groups were established: XMTE VS YL and AGS VS YL, with YL as the experimental group and XMTL and AGS as the control groups. XMTE VS AGS, was conducted with AGS as the experimental group and XMTL as the control group.

2.2. Transcriptome Sequencing and Bioinformatic Analysis

2.2.1. RNA Extraction and Quality Assessment

Total RNA was extracted using TRIzol™ Reagent (Invitrogen, Carlsbad, CA, USA). RNA concentration and purity (A260/A280 ratio) were assessed using a NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA), and RNA integrity was evaluated with an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). Only high-quality RNA samples with an RNA Integrity Number (RIN) ≥ 7.0 and an A260/A280 ratio between 1.8 and 2.2 were used for subsequent analysis.

2.2.2. Library Construction and Sequencing

Sequencing libraries were constructed from 1 μg of total RNA per sample using the Hieff NGS Ultima Dual-mode mRNA Library Prep Kit for Illumina (Yeasen Biotechnology, Shanghai, China). The procedure included mRNA enrichment using oligo(dT) magnetic beads, mRNA fragmentation, first- and second-strand cDNA synthesis, end-repair, A-tailing, adapter ligation, and PCR amplification. Library quality was assessed on an Agilent 2100 system, and sequencing was performed on an Illumina NovaSeq 6000 platform to generate 150 bp paired-end (PE150) reads.

2.2.3. Data Quality Control and Analysis

Raw data were obtained in FASTQ format and filtered to generate high-quality clean data. The filtering criteria included removing reads containing adapter sequences, reads with an N ratio > 10%, and reads where bases with a quality score Q ≤ 10 accounted for more than 50% of the entire read length. A Q30 percentage of ≥97% was required for all samples to ensure data reliability.

2.2.4. Sequence Alignment and Gene Expression Analysis

Clean reads were aligned to the bovine reference genome (ARS-UCD2.0) using Hisat2 software. The alignment efficiency for each sample ranged from 87.76% to 98.60%, indicating good data quality. Based on the alignment, transcripts were predicted and assembled using StringTie with the Reference-Based Annotation Transcript (RABT) assembly method, identifying 7915 novel genes. Gene functions were annotated using a multi-database strategy, including Nr, Pfam, KOG/COG, Swiss-Prot, KO, and Gene Ontology (GO). Gene expression levels were quantified and expressed as Fragments Per Kilobase of transcript per Million mapped reads (FPKM).

2.2.5. Differential Expression Analysis

Differential expression analysis was performed using the DESeq2 package, which is based on a negative binomial distribution model. The criteria for identifying differentially expressed genes (DEGs) were a fold change (FC) ≥ 2 and a false discovery rate (FDR) < 0.01. Functional enrichment analysis of DEGs was conducted using GO and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis, implemented with the clusterProfiler package based on the Wallenius non-central hyper-geometric distribution. A corrected p-value (FDR) < 0.05 was set as the threshold for significant enrichment.

2.3. Lipidomics Analysis

2.3.1. Lipid Extraction and Sample Preparation

Total lipids were extracted from longissimus dorsi samples using a modified Bligh-Dyer method. Briefly, 100 mg of tissue was accurately weighed into a 2 mL centrifuge tube, homogenized with a chloroform/methanol (2:1, v/v) solution, and subjected to sonication and centrifugation. The resulting organic phase was collected, dried under a nitrogen stream, and stored at −80 °C until analysis.

2.3.2. Ultra-High Performance Liquid Chromatography-Mass Spectrometry (UHPLC-MS/MS) Analysis

Lipid separation was performed on a Waters ACQUITY UPLC I-Class system equipped with an ACQUITY UPLC HSS T3 column (100 mm × 2.1 mm, 1.8 μm; Waters, Milford, MA, USA). Mobile phase A was an acetonitrile/water mixture (6:4, v/v) containing 10 mM ammonium acetate and 0.1% formic acid; mobile phase B was an isopropanol/acetonitrile mixture (9:1, v/v) containing 10 mM ammonium acetate and 0.1% formic acid. The gradient elution program was as follows: 0–2 min, 10% B; 2–15 min, 10–85% B; 15–24 min, 85–100% B; 24–30 min, 100% B; 30–31 min, 100–10% B; 31–35 min, 10% B. The column temperature was maintained at 35 °C, the flow rate was 0.3 mL/min, and the injection volume was 5 μL.

Mass spectrometry data were acquired on a SCIEX TripleTOF 5600+ mass spectrometer (SCIEX, Framingham, MA, USA) equipped with an electrospray ionization (ESI) source, operated in both positive and negative ion modes. Source parameters were: ion spray voltage, 5500 V (positive)/−4500 V (negative); temperature, 550 °C; ion source gas 1 and 2, 50 psi; curtain gas, 40 psi. Data were collected using an information-dependent acquisition (IDA) mode, with a TOF-MS scan range of m/z 100–1500 and an MS/MS scan range of m/z 50–1500.

2.3.3. Data Processing and Statistical Analysis

Raw MS data were processed using LipidSearch software (v4.1, Thermo Scientific) for peak extraction, alignment, lipid identification, and relative quantification. Lipids were identified by matching the accurate mass and MS/MS fragmentation patterns against the software’s built-in database. Principal Component Analysis (PCA) and Orthogonal Partial Least Squares Discriminant Analysis (OPLS-DA) were used to assess inter-group differences and intra-group variability. Differential lipid molecules (DLMs) were identified based on the following criteria: Variable Importance in Projection (VIP) > 1.0, Fold Change ≥ 1.5, and a Benjamini–Hochberg adjusted p-value < 0.05.

2.4. Targeted Amino Acid Metabolomics Analysis

2.4.1. Amino Acid Extraction and Sample Preparation

Precisely 25 mg of frozen muscle tissue was weighed, and 1 mL of a pre-chilled methanol:acetonitrile:water (2:2:1, v/v/v) extraction solution containing an isotope-labeled internal standard mixture was added. The mixture was homogenized (45 Hz, 4 min) and sonicated in an ice bath (5 min), with the process repeated three times. After incubation at −20 °C for 1 h to precipitate proteins, the sample was centrifuged at 12,000 rpm for 15 min at 4 °C. An aliquot of 200 μL of the supernatant was filtered through a 0.22 μm membrane and transferred to a sample vial for analysis.

2.4.2. UHPLC-MRM-MS/MS Analysis

Separation was performed on a Waters ACQUITY I-Class UHPLC system equipped with an ACQUITY UPLC BEH Amide column (150 × 2.1 mm, 1.7 μm; Waters, Milford, MA, USA). Mobile phase A was water with 0.2% formic acid, and mobile phase B was acetonitrile with 0.2% formic acid. An optimized gradient elution program was used for efficient separation. The column and autosampler temperatures were maintained at 35 °C and 10 °C, respectively, with an injection volume of 1 μL.

Mass spectrometric detection was carried out on a SCIEX QTRAP 6500+ triple quadrupole mass spectrometer (SCIEX, Framingham, MA, USA) in Multiple Reaction Monitoring (MRM) mode. Ion source parameters were: Curtain Gas = 35 psi, IonSpray Voltage = +5500 V/−4500 V, Temperature = 550 °C, Ion Source Gas 1 = 50 psi, and Ion Source Gas 2 = 55 psi. For each target metabolite, optimal ion pairs were selected for quantification by optimizing precursor ion selection and collision energy, while other pairs were used for qualitative confirmation. A total of 37 amino acid-related target metabolites were measured in this study (Table 1).

2.4.3. Method Validation and Data Analysis

The method’s linearity, sensitivity, precision, and accuracy were evaluated using a series of calibration solutions and quality control (QC) samples. Calibration curves were generated using a weighted (1/x) least-squares linear regression model. The lower limit of detection (LLOD) and lower limit of quantification (LLOQ) were defined as concentrations corresponding to signal-to-noise ratios of 3 and 10, respectively. Method precision was assessed by the relative standard deviation (RSD) of repeated QC sample injections, and accuracy was evaluated by spike-and-recovery tests.

Multivariate statistical methods, including PCA and OPLS-DA, were applied to the sample data. Differential metabolites were identified using the criteria: fold change ≥ 1, VIP value ≥ 1, and a p-value < 0.01 from a t-test. All statistical analyses were performed using R software (version 4.0.3).

2.5. Multi-Omics Data Integration Analysis

Weighted Gene Co-expression Network Analysis (WGCNA) was performed using the WGCNA R package (v1.70-3) to integrate transcriptomic and metabolomic data. A scale-free co-expression network was constructed from the gene expression matrix of all samples with the following parameters: soft-thresholding power = 8, network type = signed, minimum module size (minModuleSize) = 30, and module merging threshold (mergeCutHeight) = 0.25.

The Module Eigengene (ME) was calculated for each module, and Pearson correlation coefficients were computed between MEs and the content of differential lipids and amino acids. Key modules significantly associated with specific metabolite phenotypes (|r| > 0.6, p < 0.05) were selected for further analysis. GO and KEGG functional enrichment analyses were performed on the genes within these key modules.

A correlation network between DEGs and differential metabolites was constructed based on Pearson correlation (|r| > 0.7, p < 0.01) and visualized using Cytoscape software (v3.8.2). Unless otherwise specified, all statistical analyses were performed in R software (v4.0.3). Data are presented as mean ± standard deviation (SD).

3. Results

3.1. Transcriptome Data Analysis

3.1.1. Quality Assessment and Overall Features of Transcriptome Data

Following a stringent quality control process, a total of 177.44 Gb of high-quality clean data was generated. The percentage of Q30 bases exceeded 97.20% for all samples, providing a high-depth and high-fidelity dataset for subsequent analyses. To evaluate data reliability and consistency, we analyzed the gene expression distributions. The results indicated that all samples exhibited similar density distributions and dispersion levels of gene expression (measured by FPKM values), suggesting no significant technical bias between sequencing batches (Figure 2 and Figure 3).

Spearman’s correlation analysis was performed to assess the reliability of biological replicates. The heatmap revealed high correlations among biological replicates within the same breed, confirming the consistency of sample processing and the robustness of the experimental design. Conversely, lower correlations were observed between different breeds, indicating distinct breed-specific transcriptomic profiles (Figure 4).

3.1.2. Differential Gene Expression Analysis

Pairwise differential expression analysis was conducted among the three breeds using DESeq2, with the criteria of Fold Change ≥ 2 and FDR < 0.01. As shown in Figure 5, a substantial number of differentially expressed genes (DEGs) were identified between YL and the other two breeds. Compared to Simmental (XMTE, a non-marbled breed), 2053 DEGs were identified (1333 upregulated, 720 downregulated). Compared to Angus (AGS, a marbled breed), 2156 DEGs were identified (1230 upregulated, 926 downregulated).

Notably, only 854 DEGs (402 upregulated, 452 downregulated) were found between Angus, a marbled breed, and Simmental, a non-marbled breed. The number of DEGs between Yunling and the other two breeds far exceeded the number between Angus and Simmental. This result strongly suggests that YL possess a unique transcriptional regulatory program. The significant divergence in transcriptomic profiles highlights a distinct molecular regulatory network governing intramuscular fat deposition in YL, providing the genetic basis for its characteristic marbling phenotype.

3.2. Lipidomics Data Analysis

3.2.1. Overall Lipidomic Features and Differential Analysis

In the lipidomic analysis of the longissimus dorsi muscle from the three cattle breeds, a total of 6245 high-quality lipid molecules were annotated, spanning 8 major lipid categories and 29 subclasses from the LIPID MAPS database (Figure 6a). Chi-square analysis revealed significant differences in the distribution of lipid subclasses among the three breeds (p < 0.01). In the XMTE vs. YL comparison, sphingolipids (SPs, 13.7% vs. 18.7%) and glycerophospholipids (GPs, 19.6% vs. 22.7%) were significantly enriched in YL (Figure 6b). In the XMTE vs. AGS comparison, fatty acyls (FAs, 21.2% vs. 14.9%) and GPs (19.8% vs. 38.9%) showed marked differences (Figure 6c). In the AGS vs. YL comparison, FAs (31.6% vs. 22.1%) and glycerolipids (GLs, 22.4% vs. 13.4%) differed significantly (Figure 6d). Both Principal Component Analysis (PCA) (Figure 6e) and correlation analysis (Figure 6f) demonstrated that samples clustered within breeds and were clearly separated between breeds, confirming significant differences in their lipid compositions.

3.2.2. Analysis and Functional Enrichment of DLMs

Volcano plots visualized the distribution of differential lipid molecules (DLMs) in each comparison group (Figure 7c–e). A total of 784 DLMs (391 upregulated, 393 downregulated) were identified in the XMTE vs. YL comparison. In contrast, the number of DLMs decreased to 410 (152 upregulated, 258 downregulated) in the AGS vs. YL comparison (Figure 7a,b). Notably, the number of DLMs between XMTE vs. YL and XMTE vs. AGS was substantially greater than that between AGS vs. YL, suggesting a higher similarity in lipid composition between Angus and YL.

Lipid subclass enrichment analysis (p < 0.01) revealed that, compared to XMTE, the longissimus dorsi of YL was significantly enriched in several functional lipid subclasses (Figure 8). For example, the significant upregulation of polyunsaturated fatty acids (FA08) and specific glycerophospholipids (e.g., GP03, GP15) implies more active membrane remodeling and signal transduction in YL muscle, which are crucial for regulating adipocyte differentiation and lipid storage. Conversely, the downregulation of subclasses like cholesterol (ST04) and prenol lipids (PR01) may reflect distinct sterol metabolic pathways. The XMTE vs. AGScomparison showed a similar enrichment pattern (Figure 9), further indicating that the enrichment of specific lipid subclasses is a common feature of the marbling phenotype. In the direct comparison between YL and AGS, the differences were primarily observed in diacylglycerols (GL02), phosphatidylcholines (GP01), and sphingoid bases (SP05) (Figure 10), indicating that even among marbled breeds, YL possess a unique network for lipid metabolism fine-tuning.

3.2.3. Identification and Characterization of Major Lipid Molecules

PCA identified glycerophospholipids (GP01) and glycerolipids (GL03) as the major lipid subclasses in the longissimus dorsi muscle of all three breeds (Figure 11). Further PCA on the molecular composition within each lipid subclass allowed for the identification of characteristic major lipid molecules for each breed. To pinpoint the core differential lipids in each comparison, an integrated analysis combining the screened DLMs with the major lipid molecules was performed using Venn diagrams (Figure 12), identifying group-specific key lipid molecules.

In the XMTE vs. YL comparison, 27 major differential lipid molecules were identified, with 17 upregulated and 10 downregulated in YL (Table 1). These included significantly upregulated phosphatidylglycerols (PG), sphingomyelins (SM), and diacylglycerols (DG), as well as significantly downregulated cardiolipins (CL).

In the XMTE vs. AGS comparison, 11 major differential lipid molecules were identified, with 8 upregulated and 3 downregulated in AGS (Table 2). The upregulated molecules were primarily specific sphingolipids and bile acid-related compounds.

In the AGS vs. YL comparison, 17 major differential lipid molecules were identified, with 12 upregulated and 5 downregulated in YL (Table 3). Notably, galactosylceramide (GalCer) and certain sphingomyelins (SM) were significantly downregulated in YL compared to AGS.

3.3. Targeted Amino Acid Metabolomics Data Analysis

Based on stringent criteria (Fold Change ≥ 1, VIP ≥ 1, and p < 0.01), differential amino acid metabolites were systematically identified across the three comparison groups.

In the XMTE vs. YL comparison, Beta-Alanine (8385.83 vs. 14,609.18, p = 0.00097) and L-Aspartic acid (49,522.08 vs. 109,555.15, p = 0.0016) were identified as being significantly downregulated in YL cattle (Figure 13a, Table 4).

In the AGS vs. YL comparison, a range of amino acids, including L-Alanine, Beta-Alanine, L-Methionine, 1-Methyl-L-histidine, 3-Methyl-L-histidine, L-Arginine, and L-Carnosine, were found to be significantly downregulated in YL (p < 0.01) (Figure 13c, Table 5).

In the XMTE vs. AGS comparison, levels of L-Alanine, L-Proline, L-Lysine, L-Methionine, and L-Arginine were significantly higher in AGS than in XMTE (p < 0.01) (Figure 13b, Table 6). These differential amino acids play important roles in muscle pH regulation, energy metabolism, and antioxidant capacity, and may be key metabolic factors shaping the marbling characteristics of YL.

3.4. Integrated Analysis of Transcriptome and Lipidome

3.4.1. Identification of Key Gene Modules for Marbling-Related Lipid Formation

To identify key gene modules and networks involved in marbling formation in YL cattle, we applied WGCNA to the DEGs from all three breeds (Figure 14). Genes were clustered into distinct modules, each designated by a color. We then calculated the correlation between each module and the previously identified characteristic lipid levels to pinpoint the key modules most relevant to marbling in YL.

In the XMTE vs. YL comparison, the blue and brown modules showed broad and strong positive correlations (p < 0.01) with the majority of characteristic lipids upregulated in YL, such as sphingomyelins and diacylglycerols. This indicates that the gene sets within these two modules are key positive regulators driving marbling formation in YL (Figure 15). In the AGS vs. YL comparison, the brown and blue modules were also significantly correlated with multiple lipids that differentiate the two marbled breeds (Figure 15), suggesting that YL cattle possess unique lipid metabolism regulatory pathways even when compared to the high-quality marbled AGS breed. As a reference, the XMTE vs. AGS comparison also identified modules (e.g., green, blue) positively correlated with marbling-related lipids in AGS (Figure 15), validating the effectiveness of this analytical strategy. Modules with significant co-variation with characteristic lipids were defined as candidate key modules, providing precise targets for subsequent functional analysis.

3.4.2. GO Functional Enrichment Analysis of Key Gene Modules

To elucidate the biological functions of these key modules, we performed GO and KEGG pathway enrichment analyses. In the XMTE vs. YL comparison, the brown module was significantly enriched in the PI3K-Akt, MAPK, and Ras signaling pathways (Figure 16b,e). These classic pathways regulate cell proliferation, differentiation, survival, and metabolism, suggesting that YL cattle enhance the proliferation and differentiation of intramuscular preadipocytes and reinforce metabolic flux for lipid synthesis by augmenting responses to upstream signals like growth factors and insulin. Concurrently, genes in the blue module were concentrated in processes such as immune response, cell adhesion, and tissue repair (Figure 16a,d). These findings suggest that YL cattle actively remodel the muscle microenvironment to create a niche conducive to adipocyte infiltration, expansion, and maturation. The crosstalk between immunity and metabolism (immuno-metabolism) and the dynamic regulation of the extracellular matrix may be critical for forming a stable, well-defined marbling pattern.

In the AGS vs. YL comparison, functional annotation of the gene modules highly correlated with characteristic lipids revealed that the brown module was primarily enriched in biological processes such as regulation of axonal transport, organ growth, and mitochondrial biogenesis (Figure 17a). The blue module was significantly enriched in mRNA splicing and snRNP assembly, suggesting its potential involvement in bone tissue metabolism and immune regulation (Figure 17c).

3.5. Integrated Analysis of Transcriptome and Amino Acid Metabolome

3.5.1. Identification of Key Gene Modules for Characteristic Amino Acid Formation

By calculating the Pearson correlation coefficients between module eigengenes (MEs) and the levels of differential amino acids, we identified several key modules significantly associated with specific amino acid metabolites (Figure 18).

In the XMTE vs. YL comparison, the blue module (MEblue) showed a strong negative correlation with the two amino acids that were significantly lower in YL: Beta-Alanine (r = −0.44, p = 0.07) and L-Aspartic acid (r = −0.56, p = 0.02) (Figure 18d).

In the AGS vs. YL comparison, more complex correlation patterns emerged. The blue module (MEblue) was significantly and positively correlated with several amino acids that were lower in YL, including L-Alanine (r = 0.73, p = 6 × 10^−4^), L-Methionine (r = 0.73, p = 6 × 10^−4^), 3-Methyl-L-histidine (r = 0.88, p = 2 × 10^−6^), L-Arginine (r = 0.65, p = 0.003), and L-Carnosine (r = 0.60, p = 0.008). Additionally, the yellow module (MEyellow) showed a strong positive correlation with Beta-Alanine (r = 0.76, p = 3 × 10^−4^) and 3-Methyl-L-histidine (r = 0.61, p = 0.008) (Figure 18f).

For reference, the XMTE vs. AGScomparison identified modules regulating amino acids associated with a high marbling grade. The green module (MEgreen) was highly positively correlated with L-Alanine (r = 0.71, p = 9 × 10^−4^) and L-Methionine (r = 0.78, p = 1 × 10^−4^), while the gray module (MEgrey) was significantly positively correlated with L-Proline (r = 0.69, p = 0.002) and L-Arginine (r = 0.76, p = 3 × 10^−4^) (Figure 18e). These results clearly link specific gene expression networks to breed-specific amino acid metabolic phenotypes.

3.5.2. Functional Enrichment Analysis of Key Gene Modules

To investigate the biological functions of these key modules, we performed GO and KEGG enrichment analyses on their constituent genes.

The blue module, associated with differential amino acids (Beta-Alanine, L-Aspartic acid) in the XMTE vs. YL comparison, was significantly enriched in GO biological processes such as “positive regulation of ERK1 and ERK2 cascade,” “positive regulation of angiogenesis,” and “elastic fiber assembly” (Figure 19a). KEGG analysis indicated that these genes are primarily involved in “Focal adhesion” and the “PI3K-Akt signaling pathway” (Figure 19d).

Among the key modules related to differential amino acids in the AGS vs. YL comparison, the blue module’s functions pointed towards energy metabolism and protein synthesis. KEGG analysis revealed high enrichment in “Oxidative phosphorylation” and “Ribosome” pathways (Figure 19f), supported by GO analysis showing involvement in processes like “mitochondrial respiratory chain complex assembly” (Figure 19c). The yellow module was enriched in KEGG pathways including “Thermogenesis” and “Prion diseases.”

For context, functional analysis of key modules from the XMTE vs. AGScomparison showed that the green module was enriched in “Oxidative phosphorylation” and “Ribosome” pathways (Figure 19e), while genes in the blue module were associated with “sarcomere organization” (Figure 19b). These findings provide clear functional annotations for the molecular events underlying the differences in amino acid metabolism among the breeds.

4. Discussion

4.1. Transcriptional Regulatory Features of Marbling Formation in YL

GO enrichment analysis of differentially expressed genes (DEGs) between Simmental (XMTE) and Angus (AGS) cattle revealed that upregulated genes in AGS were primarily enriched in “cellular process,” “biological regulation,” and “metabolic process,” suggesting higher activity in these processes in AGS cattle (Figure 20b). In the Cellular Component (CC) category, DEGs were concentrated in “cellular anatomical entity” and “intracellular,” indicating their involvement in intracellular and structural functions. In the Molecular Function (MF) category, upregulated genes were significantly enriched in “binding activity” and “catalytic activity,” implying a higher metabolic capacity in AGS. KEGG enrichment analysis further elucidated the roles of these DEGs in several important signaling pathways (Figure 20e). DEGs were significantly enriched in the PI3K-Akt, MAPK, and Wnt signaling pathways, which play central roles in regulating cell proliferation, differentiation, and apoptosis [10]. This suggests that AGS cattle may promote lipid metabolism and cellular functions through these pathways. DEGs were also enriched in “oxidative phosphorylation” and lipid metabolism-related pathways, indicating unique mechanisms for energy metabolism and lipid regulation in AGS. In summary, DEGs in AGS are mainly involved in cell metabolism, signal transduction, and energy metabolism, highlighting their heightened activity in lipid metabolism and cellular functions.

In the AGS vs. YL comparison, GO enrichment analysis showed that DEGs were primarily involved in “cellular process” and “metabolic process” (Biological Process, BP), as well as “cell” and “organelle” (CC) (Figure 20c). This suggests that cellular metabolism in the adipose tissue of YL longissimus dorsi may be more active. Furthermore, DEGs related to the “immune system process” were distinct in YL, implying potentially superior disease resistance and meat quality [11]. KEGG analysis revealed that DEGs between AGS and YL were concentrated in critical signaling pathways such as PI3K-Akt and apoptosis (Figure 20f). The PI3K-Akt pathway is a key regulator of cell growth and survival [12], while apoptosis is crucial for tissue homeostasis. Therefore, differential expression of genes in these pathways may directly influence adipose development in YL, leading to its characteristic marbled fat distribution.

4.2. The Unique Lipid Molecular Basis of Marbling in YL

As shown in Figure 11, Principal Component Analysis (PCA) of lipid subclasses revealed that glycerophospholipids (GP01) and glycerolipids (GL03) are the predominant lipid subclasses in the longissimus dorsi of all three breeds. Further PCA on the lipid molecules within each subclass identified the major lipid molecules characteristic of each breed. Numerous studies have shown that neutral lipids and phospholipids, such as GPs and GLs, are essential components of cell membranes and lipid droplets, playing critical roles in regulating adipocyte differentiation [5], hypertrophy [13], and triglyceride accumulation [14].

Functional analysis of characteristic lipids in the XMTE vs. YL comparison identified 27 major differential lipid molecules, with 17 upregulated and 10 downregulated in YL. Upregulated lipids included specific phospholipids, sphingolipids, and glycerolipids such as PG(20:5(5Z,8Z,11Z,14Z,17Z)/22:0), SM(d20:0/24:1), and DG(16:0/18:1(11Z)/0:0). Phospholipids are the main components of cell membranes, and their unsaturated fatty acid content is closely related to membrane fluidity [15]. Sphingolipids are key components of lipid rafts and are involved in sensing and transducing signals from molecules like insulin [16]. This suggests that specific lipids like PG and SM may influence the differentiation, maturation, and triglyceride synthesis of YL adipocytes by altering their membrane properties and signal transduction. Among the downregulated lipids in YL compared to XMTE, CL(1′-[16:0/18:1(9Z)],3′-[18:1(9Z)/20:4(5Z,8Z,11Z,14Z)]), a cardiolipin, plays a crucial role in energy metabolism in heart and muscle tissues. Its downregulation might reduce the cell’s ability to respond to stress, thereby affecting energy production and storage in YL cattle. Meanwhile, the downregulation of 1-(O-alpha-D-glucopyranosyl)-27-keto-(1,3R,29R)-triacontanetriol, potentially involved in carbohydrate metabolism and hormone signaling, could affect glucose metabolism and endocrine regulation in YL.

In the AGS vs. YL comparison, 17 major differential lipid molecules were identified (12 upregulated, 5 downregulated in YL). While sharing some highly expressed differential lipids with the XMTE vs. YL comparison, we focused on the downregulated lipids, which included GalCer(d18:0/26:0) and SM(d18:2/18:1). Reduced levels of galactosylceramide (GalCer) and sphingomyelin (SM) could affect cell membrane integrity and signal transduction in YL. These changes might also reflect cellular adaptation to environmental stress, leading to a readjustment of metabolic pathways.

4.3. Amino Acid Metabolic Patterns in YL and Their Potential Impact on Meat Quality

Beta-Alanine, a precursor of carnosine, is critical for muscle pH buffering, antioxidant defense, and fatigue resistance [17]. Its level is closely related to the pH stability and anti-fatigue capacity of muscle tissue [18]. The significantly lower level of Beta-Alanine in YL (XMTE vs. YL: 8385.83 vs. 14,609.18, p < 0.001; AGS vs. YL: 8385.83 vs. 17,744.94, p < 0.01) may reflect a different pH regulation mechanism in their muscle. This difference could alter the rate of post-mortem pH decline and the ultimate pH, thereby affecting meat tenderness, water-holding capacity, and oxidative stability. In particular, lower Beta-Alanine levels might be associated with the unique texture and flavor of YL beef, as meat pH significantly influences protein degradation and flavor compound formation.

L-Aspartic acid, a key intermediate in the TCA cycle, participates in energy metabolism and amino acid transamination, and plays an important role in energy metabolism and nervous system function [19]. The significantly lower level of L-Aspartic acid in the longissimus dorsi of YL cattle (XMTE vs. YL: 49,522.08 vs. 109,555.15, p < 0.01) suggests that their muscle tissue may employ a different energy metabolism strategy. This metabolic divergence could affect the energy efficiency and metabolic state of muscle cells, shaping the unique muscle structure and meat quality of YL cattle.

L-Methionine, an essential amino acid and methyl donor, is involved in one-carbon metabolism, epigenetic modifications, and protein synthesis. [20] showed that L-Methionine content is positively correlated with muscle yield and protein synthesis efficiency. The lower expression of L-Methionine in YL cattle (AGS vs. YL: 15,781.16 vs. 24,967.22, p < 0.01) may reflect differences in protein synthesis and methylation regulation, affecting muscle growth patterns and protein composition, which in turn influences meat tenderness and texture. Furthermore, as a sulfur-containing amino acid, L-Methionine is a precursor to many sulfur-containing flavor compounds; its variation could directly impact the formation of the unique flavor profile of YL marbled beef [21].

Additionally, L-Arginine and L-Carnosine levels were also significantly lower in the longissimus dorsi of YL compared to AGS (L-Arginine: 66,012.22 vs. 112,079.70, p < 0.001; L-Carnosine: 8,795,215.04 vs. 13,517,703.54, p < 0.01). L-Arginine, a precursor for nitric oxide synthesis, is involved in vasodilation, immune regulation, and protein synthesis [22]. Ref. [23] found that L-Arginine significantly affects muscle growth and fat deposition by modulating the nitric oxide signaling pathway. The lower L-Arginine level in YL may be associated with its unique IMF deposition pattern, an important factor in marbling formation. The reduced level of L-Carnosine, a dipeptide with strong antioxidant activity, could affect the oxidative stability and free-radical scavenging capacity of muscle tissue [21], potentially influencing the quality stability of YL beef during storage and processing.

4.4. Integrated Analysis Reveals Mechanisms of Marbling Formation

By integrating multi-omics data, we constructed a regulatory network that not only explains the molecular basis of the superior marbling potential in YL cattle but also provides insights into the factors contributing to its inconsistent expression.

(1) Core Regulatory Role of Signal Transduction Networks. WGCNA linked key gene modules with characteristic metabolites, precisely identifying the core signaling pathways that regulate fat deposition in YL. In the comparison with Simmental, gene modules strongly and positively correlated with the accumulation of characteristic lipids (e.g., sphingomyelins, diacylglycerols) were significantly enriched in the PI3K-Akt, MAPK, and Ras signaling pathways. These pathways are central hubs for cell proliferation, differentiation, and metabolic regulation. Key genes within these modules, such as Integrin Subunit Beta 3 (ITGB3), Platelet-Derived Growth Factor Receptor Alpha (PDGFRA), and Fibronectin 1 (FN1), act as important nodes that receive and transduce signals from growth factors and the extracellular matrix. This systematically promotes the proliferation and differentiation of intramuscular preadipocytes, providing the cellular basis for lipid synthesis [10]. This indicates that YL establish an efficient pro-adipogenic regulatory network by enhancing upstream signal responses.

(2) Energy Metabolism Reprogramming Fuels Lipid Synthesis. Efficient lipid synthesis depends on an adequate supply of energy (ATP) and precursors (acetyl-CoA). This study found that several key gene modules associated with Yunling’s characteristic metabolites pointed to a systematic reprogramming of energy metabolism. In the comparison with Angus, gene modules associated with lipid profile differences were significantly enriched in pathways like oxidative phosphorylation and the TCA cycle, involving core genes such as NDUFB1 (a mitochondrial respiratory chain complex I subunit), COX7C (a complex IV subunit), and IDH3B (isocitrate dehydrogenase). Differential expression of these genes directly impacts mitochondrial function and ATP synthesis efficiency. Simultaneously, gene modules associated with the amino acid profile (particularly the lower level of L-Aspartic acid) were also enriched in the oxidative phosphorylation pathway. This reveals that Yunling muscle tissue may optimize mitochondrial function not only to supply energy and carbon sources for de novo fatty acid synthesis but also to alter the amino acid metabolic flux closely linked to the TCA cycle, collectively shaping its unique dual metabolic signature in both lipids and amino acids [24].

(3) Microenvironment Remodeling Promotes Fat Deposition. Beyond intracellular signaling and metabolism, this study also highlights the critical role of the cellular microenvironment in Yunling marbling formation. Gene modules associated with differential amino acids were significantly enriched in the “Focal adhesion” pathway, and genes like ITGB3 and FN1 were involved in both signal transduction and cell–matrix adhesion. Focal adhesions are key structures connecting the cytoskeleton to the extracellular matrix (ECM), crucial for sensing and transducing mechanical signals and regulating cell migration and tissue homeostasis. Dynamic remodeling of the ECM has been confirmed as a key factor in regulating adipose tissue development [25]. Therefore, the coordinated regulation of these genes suggests that YL may remodel the ECM structure between muscle fibers, creating a microenvironment that is physically and biochemically conducive to adipocyte infiltration, survival, and maturation, thereby promoting the formation of a diffuse marbling pattern.

In summary, in the longissimus dorsi of YL, a regulatory network centered on the PI3K-Akt and MAPK signaling pathways is activated within a remodeled, supportive cellular microenvironment. This drives a reprogramming of energy metabolism characterized by optimized oxidative phosphorylation. These synergistic events lead to the specific accumulation of characteristic lipid molecules and a systematic alteration of the amino acid profile, ultimately resulting in the superior marbling quality of Yunling beef.

4.5. Limitations and Future Perspectives

This study is based on static information from a specific developmental stage. While the correlation analysis reveals key molecular networks, it does not establish causality. The limited number of biological replicates and the sampling from a single muscle site restrict the generalizability of the findings. Therefore, future research must validate the functions of the core regulatory factors (e.g., NDUFB1, ITGB3, PDGFRA) and characteristic metabolites (e.g., SM(d20:0/24:1), Beta-Alanine) identified in this study through functional experiments, such as using bovine intramuscular preadipocyte cell lines combined with gene-editing technologies (CRISPR/Cas9). Furthermore, longitudinal studies spanning different growth stages (from calf to finishing) and including multiple muscle sites are crucial for comprehensively capturing the dynamic regulation of intramuscular fat deposition.

Future research should aim to deconstruct the cellular and molecular heterogeneity of marbling formation at higher resolutions. Integrating single-cell RNA sequencing (scRNA-seq) with spatial transcriptomics will enable the precise dissection of the unique contributions of different cell subpopulations, such as intramuscular adipocytes, muscle fibers, fibroblasts, and immune cells, and their spatial interactions during marbling formation. Additionally, incorporating proteomic and epigenetic data (e.g., DNA methylation, histone modifications) into the existing multi-omics framework will allow for the construction of a multi-layered, panoramic regulatory network, from epigenetic modifications to transcription, post-translational modifications, and metabolic products.

5. Conclusions

Through an integrated analysis of the transcriptome, lipidome, and amino acid metabolome, we provide initial insights into the mechanisms of lipid and amino acid formation in the marbled beef of YL. Key findings are as follows: (1) At the transcriptomic level, YL exhibit a unique gene expression pattern, with DEGs significantly enriched in key signaling pathways regulating adipocyte proliferation and differentiation, such as PI3K-Akt and MAPK, as well as the oxidative phosphorylation pathway, which is closely linked to energy metabolism. (2) At the lipidomic level, the longissimus dorsi of YL is rich in sphingolipids (SPs) and glycerophospholipids (GPs). Its distinct lipid profile, particularly the accumulation of characteristic molecules like PG(20:5/22:0) and SM(d20:0/24:1), forms the material basis for its superior flavor and texture. (3) At the amino acid metabolic level, YL display a metabolic pattern characterized by significantly lower levels of several functional amino acids, including Beta-Alanine and L-Aspartic acid. This may affect muscle pH, antioxidant capacity, and energy efficiency, thereby shaping their unique meat quality traits. (4) At the integrated network level, WGCNA constructed a gene-metabolite correlation network, identifying key gene modules strongly associated with characteristic lipids and amino acids. Core genes such as NDUFB1, COX7C, and ITGB3 synergistically regulate lipid and amino acid metabolism by modulating energy metabolism and cell signaling pathways, collectively driving the formation of marbling in YL.

This study not only provides a comprehensive molecular blueprint for marbling in YL cattle but also transforms our understanding from a descriptive characterization to a functional framework for improvement. The identified key genes and pathways serve as high-priority candidate targets for future molecular breeding programs and precision feeding strategies. Ultimately, this research provides a vital theoretical foundation to unlock the full genetic potential of Yunling cattle, aiming to elevate it into a consistent, high-value beef product, thereby enriching the diversity of the global high-end beef market.