For this experiment, 120 healthy 2-month-old Tibetan lambs with similar body weight (16.87 ± 0.31 kg) were randomly divided into 4 treatment groups (30 lambs per group), with 5 replicates per group (6 lambs per replicate). The four treatment groups were: basal diet group (H group), basal diet + resveratrol (1.5 g/day) group (H-RES group), basal diet + β-hydroxy-β-methylbutyrate (1.25 g/day) group (H-HMB group), and basal diet + resveratrol (1.5 g/day) + β-hydroxy-β-methylbutyrate (1.25 g/day) group (H-RES-HMB group). The experimental period lasted 100 days, consisting of a 10-day adaptation period followed by a 90-day formal trial. All lambs were housed in well-ventilated, dry pens with access to an exercise yard allowing free movement. The pens were thoroughly disinfected prior to the experiment. According to the experimental design, each treatment group received feed twice daily (08:00 and 17:00) (Table 1).

Following slaughter, adipose tissue samples were collected under aseptic conditions, immediately flash-frozen in liquid nitrogen, and stored at −80 °C until further analysis. Fatty acid analysis was performed using the hydrolysis-extraction method described in the national standard “Determination of Fatty Acids in Food” (GB 5009.168-2016). Briefly, an internal standard (undecanoic triglyceride solution) was added to each adipose sample, which was then hydrolyzed with hydrochloric acid. The liberated lipids were extracted into an ether phase, followed by saponification and methylation under alkaline conditions to produce fatty acid methyl esters (FAMEs). The resulting FAMEs were quantified using an Agilent GC-6890 gas chromatograph. Medium- and long-chain fatty acid contents were calculated based on the concentrations of individual FAMEs and their corresponding conversion factors.

A 1.0 mL sample of Tibetan sheep intermuscular fat was accurately weighed and placed in a 20 mL headspace vial, then incubated at 50 °C and 500 r/min for 20 min before injection. Each sample was analyzed in triplicate, with an injection volume of 500 μL (splitless mode, syringe temperature 85 °C). Gas chromatography analysis employed a column with an initial temperature of 40 °C, using high-purity nitrogen (≥99.999%) as the carrier gas. The pressure programming was as follows: initial flow rate of 2.0 mL/min held for 2 min, increased to 10.0 mL/min over 8 min, then raised to 100.0 mL/min over 10 min and held for 10 min (total run time 20 min, injection port temperature 80 °C). Ion mobility spectrometry utilized a tritium source for ionization (drift tube length 53 mm, electric field strength 500 V/cm, temperature 45 °C), with high-purity nitrogen (≥99.999%) as the drift gas at 150 mL/min in positive ion mode. The experiment first established calibration curves for retention time and retention index using a mixed standard solution of six ketones. The retention index of target compounds was calculated based on their retention times, followed by qualitative analysis using the NIST 2020 database and IMS drift time database in VOCal software. Volatile compounds were visualized using Reporter, Gallery Plot, and Dynamic PCA plugins to generate three-dimensional spectra, two-dimensional spectra, differential spectra, fingerprint plots, and PCA plots for comparative analysis of sample differences.

Adipose tissue samples measuring approximately 3 × 3 cm were fixed in 4% paraformaldehyde for histopathological examination. Adipocyte volume was assessed using standard hematoxylin and eosin (H&E) staining. Briefly, the fixed tissue was sectioned into 5 μm-thick slices and subsequently stained with H&E. Adipose tissue morphology was observed under a light microscope at 40× magnification (Olympus DP2-BSW, Tokyo, Japan).

Total RNA was extracted from liver tissue using the TRIzol reagent kit (Invitrogen, Waltham, MA, United States), and RNA quality was assessed using an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, United States) and RNase-free agarose gel electrophoresis. PolyA mRNA was enriched using Oligo(dT) magnetic beads and fragmented via ultrasonication. The fragmented mRNA was used as a template with random oligonucleotides as primers to synthesize double-stranded cDNA using M-MuLV reverse transcriptase. After end repair, adapter ligation, and size selection with AMPure XP beads, a cDNA library (~200 bp) was constructed. Clean reads were aligned to the reference sequence using HISAT 2.2.4 software (22), and gene expression levels were quantified based on the FPKM method. DEGs were identified using DESeq2 under a False Discovery Rate (FDR) of 0.05. Genes with p < 0.05 and |log2FC| ≥2 were recognized as DEGs (23). Finally, GO and KEGG enrichment analyses were performed in R based on the hypergeometric distribution (p < 0.05 considered significant). To ensure sequencing quality, stringent criteria were applied for library quality assessment: RNA integrity and DNA contamination were evaluated by agarose gel electrophoresis, RNA purity and concentration were measured via spectrophotometry (OD260/280 and OD260/230 ratios), RNA concentration was quantified using a Qubit 2.0 Fluorometer (Thermo Fisher, Waltham, MA, United States), and RNA integrity was precisely determined using the 2100 Bioanalyzer (Agilent, CA, United States).

RNA samples were reverse-transcribed into complementary DNA (cDNA) using a reverse transcription kit (Takara, Dalian, China). The standard reverse transcription reaction mixture consisted of 1 μg total RNA, 1 μL random or Oligo(dT) primers, 4 μL 5 × reverse transcription buffer, 1 μL dNTP mix, 1 μL reverse transcriptase, 0.5 μL RNase inhibitor, and RNase-free H₂O to a final volume of 20 μL. GAPDH was used as the reference gene to normalize target gene expression levels. Quantitative polymerase chain reaction (qPCR) was performed using an Applied Biosystems 7500 Fast Real-Time PCR System (Applied Biosystems, United States) with the SYBR^® Premix Ex Taq^™ Kit (TaKaRa, Dalian, China). The thermal cycling conditions were as follows: initial denaturation at 95 °C for 15 min, followed by 40 cycles of denaturation at 95 °C for 10 s and annealing at 60 °C for 30 s. The relative expression levels of target genes were calculated using the 2^−ΔΔCt method. Primer sequences are listed in Table 2.

Prior to lipid analysis, frozen adipose tissue was ground into powder using a mortar and liquid nitrogen. Lipids from subcutaneous fat were extracted via tissue homogenization using an MTBE-methanol–water three-phase solvent system. Lipid profiling was performed on a UPLC-HRMS system (Thermo Scientific, Waltham, United States) equipped with a 50 °C Accucore C30 column, using mobile phases consisting of 60% acetonitrile aqueous solution (A) containing 10 mM ammonium formate and 0.1% formic acid, and 10% acetonitrile-isopropanol solution (B). Mass spectrometry detection was conducted in both positive and negative electrospray ionization modes. Quality control (QC) samples for untargeted lipidomics were prepared by centrifuging pooled lipid extracts from equal aliquots of each sample. The analytical sequence included three experimental groups, with all extraction, identification, and quantification procedures performed by Guangzhou Gidio Biotechnology Co., Ltd.

Data processing and lipid identification were performed using LipidSearch software (Version 4.1) by matching MS2 spectra against databases of phospholipids, glycerolipids, sphingolipids, and steroids, with mass accuracies of 5 ppm for precursor ions and 5 mDa for product ions. Lipid quantification based on normalized peak area (AUC) was conducted using TraceFinder software (version 4.1). After removing duplicate lipid entries, the data were log2-transformed and subjected to orthogonal partial least squares-discriminant analysis (OPLS-DA) using SIMCA software (version 14.1; Umetrics, Umea, Sweden). Differentially accumulated lipids (DALs) were screened based on variable importance in projection (VIP) >1 and independent-sample t-test p < 0.05.

The Spearman correlation coefficient (SPCC) was calculated to analyze the association between candidate gene FPKM expression levels and metabolite contents, including key fatty acids and DALs. Gene-metabolite pairs with |SPCC| >0.5 and p < 0.05 were selected for further investigation.

The results are presented as means ± standard errors, with statistical significance defined as p < 0.05. The experimental data were statistically analyzed using SPSS 26.0 software (IBM Corp., Armonk, NY, United States) with one-way ANOVA.

Gas chromatography (GC-6890) analysis revealed significant differences (p < 0.05, Table 3) in the fatty acid composition of intramuscular fat among treatment groups (H, H-RES, H-HMB, H-RES-HMB). Among SFAs, C17:0 and C18:0 in the H and H-HMB groups were significantly higher than in the H-RES-HMB group (p < 0.05).

For MUFAs, C15:1n5 in the H-RES-HMB group was significantly higher than in the H, H-RES, and H-HMB groups (p < 0.05), while C18:1n9 in the H-RES, H-HMB, and H-RES-HMB groups was significantly higher than in the H group (p < 0.05). Among PUFAs, C18:2n6 in the H-RES and H-RES-HMB groups was significantly higher than in the H and H-HMB groups (p < 0.05), with the H-HMB group also exhibiting higher C18:2n6 than the H group (p < 0.05). C18:3n6 in the H-RES and H-RES-HMB groups was significantly higher than in the H and H-HMB groups (p < 0.05). C18:3n3 in the H-RES and H-RES-HMB groups was significantly higher than in the H group (p < 0.05). C20:3n6 and C20:3n3 in the H-RES-HMB group were significantly higher than in the H and H-HMB groups (p < 0.05), while C20:3n3 in the H-RES group was also significantly higher than in the H group (p < 0.05).

In this study, gas chromatography-ion mobility spectrometry (GC-IMS) was successfully employed for the accurate identification of volatile flavor compounds in Tibetan sheep meat. As shown in Table 4, the contents of 2-hexanone, 3-hexanone, and methyl acetate in all experimental groups (H-RES and H-HMB groups) were significantly higher than those in the control group (H group) (p < 0.05), with the H-RES-HMB group exhibiting significantly higher levels than the other two groups (p < 0.05). For 3-pentanone, the H-HMB and H-RES-HMB groups showed significantly higher contents than the H and H-RES groups (p < 0.05). The control group (H group) had significantly higher levels of 3-methyl butanal, heptanal, and n-octanal compared to the experimental groups (H-RES, H-HMB, and H-RES-HMB groups) (p < 0.05).

Histological examination of HE-stained sections under 40× light microscopy revealed that the intramuscular adipocyte areas in the H, H-RES, and H-HMB groups were significantly larger than those in the H-RES-HMB group (p < 0.01) (Figure 1A). The adipocyte diameter in the H group was markedly greater than that in both the H-RES and H-RES-HMB groups (p < 0.01), while being significantly larger than the H-HMB group (p < 0.05) (Figure 1C). The adipocyte densities in the H-RES, H-HMB, and H-RES-HMB groups were all significantly higher than that in the H group (p < 0.01) (Figure 1D). Microscopic examination revealed that all four experimental groups (Figures 1A–D) maintained normal adipose tissue morphology. Remarkably, the H-RES-HMB group demonstrated the smallest adipocyte area and diameter, along with the highest cell density values, with all these differences showing statistical significance.

To investigate the regulatory mechanisms of RES and HMB on intermuscular fat deposition in Tibetan sheep, this study performed transcriptome analysis of intermuscular adipose tissue from experimental and control groups using RNA-seq technology. After quality control filtering, the H, H-RES, H-HMB, and H-RES-HMB groups obtained 17.12, 15.79, 15.19, and 15.22 million high-quality clean reads, respectively, with mapping rates >91% for all samples (reference genome: sheep Oar_v4.03, Supplementary Table 1). Differential gene analysis results in Figure 2A showed that the H, H-RES, H-HMB, and H-RES-HMB groups contained 102, 92, 146, and 265 specifically expressed genes, respectively, while 12,562 genes were co-expressed across all four groups. The UpSet plot (Figure 2B) indicated that the H vs. H-RES-HMB group had the highest number of differentially expressed genes (2,643), with 525 shared differentially expressed genes among the three comparison groups. PCA analysis revealed that PC1 and PC2 explained 82.1 and 10.9% of the variance, respectively, with the two principal components covering 93% of the variance. The scatter points of the H-RES-HMB group were relatively more concentrated within the same quadrant (Figure 2C). Volcano plots further demonstrated that the H vs. H-RES, H vs. H-HMB, and H vs. H-RES-HMB comparisons detected 1,411 (828 upregulated, 583 downregulated), 2,560 (1,589 upregulated, 971 downregulated), and 2,643 (1,844 upregulated, 799 downregulated) differentially expressed genes, respectively (Figure 2D).

To further elucidate the functional characteristics of differentially expressed genes (DEGs), we performed systematic functional enrichment analysis. KEGG pathway analysis (Supplementary Table 2) revealed that the H vs. H-RES comparison was significantly enriched in 320 pathways, H vs. H-HMB in 335 pathways, and H vs. H-RES-HMB in 326 pathways. In the differential KEGG enrichment analysis (Figure 3), four significant pathways (p < 0.05) were commonly identified in both the H vs. H-HMB and H vs. H-RES-HMB comparisons, including the Calcium signaling pathway (ERBB4, P2RX7, ERBB3, P2RX3, and SLC8A1), Hippo signaling pathway (WNT9A, WNT10B, WNT6, and WNT2B), Estrogen signaling pathway (HSP90AA1, TGFA, and RARA), and Arachidonic acid metabolism (PLA2G4A, ALOX12, and PTGDS). Notably, these pathways were also significantly enriched in the H vs. H-RES comparison, and all four are closely associated with the regulation of lipid metabolism.

Gene set enrichment analysis (GSEA) results demonstrated that the Calcium signaling pathway (Figure 4A) and Arachidonic acid metabolism (Figure 4D) (with negative ES values) exhibited higher activity in the control group (H), whereas the Hippo signaling pathway (Figure 4B) and Estrogen signaling pathway (Figure 4C) (with positive ES values) were significantly activated in the treatment groups (H-RES, H-HMB, and H-RES-HMB). Differential expression analysis further revealed that 15 lipid metabolism-related genes, including ERBB4 and P2RX7, were significantly upregulated (p < 0.05) in the treatment groups (H-RES, H-HMB, and H-RES-HMB) compared to the control group (H).

To validate the reliability of the RNA-seq results, we selected nine genes for expression level verification using quantitative reverse transcription polymerase chain reaction (qRT-PCR). The mRNA expression levels of all genes were consistent with the RNA-seq data, confirming the reliability of the RNA-seq results (Figure 5).

Based on the previous findings that intermuscular fat deposition significantly affects the meat quality of Tibetan sheep, this study systematically identified lipid species and content through lipidomics: first, lipid molecules were characterized using UPLC-HRMS, followed by multivariate statistical analysis to reveal differences in lipid molecules among different treatment groups. This study examined the lipid metabolic characteristics of intermuscular adipocytes and explored the metabolic-related networks of the H vs. H-RES, H vs. H-HMB, and H vs. H-RES-HMB groups. Based on high-throughput LC-MS analysis with an adjusted p-value threshold of less than 0.05, a total of 499 POS-positive ion lipid metabolites and 690 NEG-positive ion lipid metabolites were detected (Supplementary Table 3). As shown in Figure 6A, in the POS ion lipid composition statistics, except for the “Others” category (10.77%), TG accounted for the largest proportion (37.47%), while LPC accounted for the smallest proportion (1.10%). Similarly, Figure 6B shows that in the NEG ion lipid composition statistics, except for the “Others” category (21.80%), PC accounted for the largest proportion (15.36%), while PG accounted for the smallest proportion (2.82%).

PCA results revealed clear separation among the three groups, including quality control samples, in the principal component (PC1 × PC2) score plots (Figures 7A,B). A total of 277 differential metabolites were identified (260 in positive mode and 17 in negative mode, Figure 7C). Multivariate analysis of subcutaneous fat metabolites showed that in the H-RES group, 19 differential metabolites exhibited significant changes (12 upregulated and 7 downregulated). In the H-HMB group, 69 differential metabolites were significantly altered (62 upregulated and 7 downregulated), while in the H-RES-HMB group, 189 differential metabolites showed significant changes (186 upregulated and 3 downregulated). OPLS-DA plots and permutation test plots confirmed substantial differences in metabolites among the three groups of Tibetan sheep intermuscular fat, demonstrating the accuracy and reproducibility of the LC-MS/MS method (Figures 8A,B). Generally, R² and Q² should exceed 0.5, with a difference of less than 0.3 between them. Under type/POS and type/NEG ion modes, t1 accounted for the largest proportions (22 and 49%, respectively). In the type/POS mode, the highest number of significantly altered metabolites in the H vs. H-RES group were TGs (45 species, 42.0%), while in the H vs. H-HMB group, TGs were also predominant (79 species, 38.5%). Similarly, in the H vs. H-RES-HMB group, TGs showed the most significant changes (135 species, 34.4%). Under the Type/NEG mode, CL and PC were the most altered metabolites in the H vs. H-RES group (11 species each, 20.0%), whereas PC and PS were predominant in the H vs. H-HMB group (24 species each, 14.1%). In contrast, Hex1Cer exhibited the highest number of changes in the H vs. H-RES-HMB group (102 species, 21.9%). Thus, it was concluded that TGs were the most abundant in the Type/POS mode, whereas Hex1Cer dominated in the type/NEG mode (lipid KEGG statistics, Supplementary Table 4). In Figure 9 and Supplementary Table 5, lipid metabolism was enriched in pathways such as glycerophospholipid metabolism, Arachidonic acid metabolism, linoleic acid metabolism, alpha-linolenic acid metabolism, and ether lipid metabolism. Global and overview maps were enriched in metabolic pathways and biosynthesis of secondary metabolites, while the nervous system was associated with retrograde endocannabinoid signaling. Cancer: overview was linked to choline metabolism in cancer. Metabolites in these pathways were predominantly LPCs and PCs, with five key metabolites identified in lipid metabolism: LPC (20:0), LPC (20:1), PC (21:2), PC (37:0), and PC (38:5).

Spearman correlation analysis was performed to calculate correlation coefficients between the expression levels of 15 candidate genes and medium- and long-chain fatty acid content as well as key lipid metabolites, enabling integrated analysis of transcriptomic and lipid metabolomic data. In Figure 10A, multiple genes exhibited significant negative correlations with saturated fatty acids (SFAs). ERBB4, RARA, PLA2G4A, SLC8A1, HSP90AA1, ERBB3, TGFA, P2RX7, WNT10B, ALOX12, WNT9A, P2RX3, WNT2B, WNT6, and PTGDS showed inverse associations with specific SFA chain lengths (C12:0–C24:0). Notably, HSP90AA1 correlated with the broadest range (C13:0–C24:0), while WNT9A and P2RX3 were linked to fewer SFAs (e.g., C13:0–C15:0 and C15:0/C20:0/C24:0, respectively). In Figure 10B, Several genes displayed significant positive correlations with polyunsaturated fatty acids (PUFAs). ERBB4 was linked to C18:2N6, while ALOX12 correlated with C18:2N6 and C18:3N3. P2RX7, ERBB3, SLC8A1, and PLA2G4A showed broad associations, spanning C18:2N6 to C22:5N6. WNT10B, WNT6, HSP90AA1, RARA, and PTGDS shared correlations with C18:2N6–C20:3N3, whereas P2RX3 and WNT9A were linked to longer-chain PUFAs (e.g., C20:4N6, C22:5N6). TGFA and WNT2B also exhibited strong PUFA associations, particularly with C18:2N6–C22:2N6.

In Figure 10C, genes (P2RX7, PTGDS, WNT2B, ERBB3, TGFA, SLC8A1, HSP90AA1, WNT6, WNT10B, PLA2G4A, RARA, ERBB4, and ALOX12) showed significant positive correlations with MUFAs (C15:1N5), while genes (P2RX7, PTGDS, WNT2B, SLC8A1, HSP90AA1, WNT6, and WNT10B) were positively correlated with MUFAs (C18:1N9), and genes (P2RX7, PTGDS, and WNT2B) exhibited positive correlations with MUFAs (C20:1N9). In Figure 10D, genes (WNT9A and WNT2B) demonstrated significant negative correlations with LPC (20:1), whereas genes (ERBB3 and ALOX12) showed positive correlations with the lipid metabolite PC (21:2). Furthermore, genes (ERBB4, PLA2G4A, HSP90AA1, RARA, ERBB3, SLC8A1, TGFA, WNT10B, P2RX7, WNT6, and ALOX12) were positively correlated with PC (37:0), and genes (ERBB4, PLA2G4A, HSP90AA1, ERBB3, SLC8A1, TGFA, WNT10B, P2RX7, WNT6, ALOX12, WNT9A, PTGDS, and WNT2B) showed positive associations with PC (38:5).