Integrated multi-omics analysis of metabolome and transcriptome profiles during bovine adipocyte differentiation reveals functional divergence of FADS2 isoforms in lipid metabolism regulation

Abstract

Background

Fat metabolism plays an important role in animal health and economic benefits. However, the changes in gene expression and metabolites during fat metabolism have not been systematically studied in bovine.

Results

This study integrates transcriptomic and metabolomic strategies to delineate the metabolic and gene expression profiles during the adipogenesis of bovine preadipocytes in four different stages. Totally, we identified 328 differentially expressed metabolites (DEMs) and 5257 differentially expressed genes (DEGs) during adipogenesis. Functional enrichment of both DEMs and DEGs highlighted the important roles of fatty acid metabolic pathways. By integrating transcriptomic and metabolomic data, we identified key genes potentially regulating fatty acid metabolism, including FADS2, ACOT7 and ACOT2. We further applied comparison for the functional differences between two FADS2 isoforms (FADS2-T0 and FADS2-T2). The results proved that the lipid metabolism regulated by FADS2-2 has changed due to the loss of 46 amino acids with a transmembrane domain, which finally altering its promoting effect on bovine fat deposition.

Conclusions

In summary, our research provides important resources and key candidate genes for a systematic understanding of the changes in gene expression and lipid metabolism during the process of fat deposition.

Background

Fat metabolism is integral to a multitude of life activities, encompassing body health in animal, which has been drawn wide attention. In livestock, it plays important roles in enhancing breeding benefits, particularly with regard to meat quality, including the formation of marble patterns and the characteristics of fatty acids [1, 2]. Adipose tissue is a crucial region for fat storage and metabolism [2, 3]. It also acts as an endocrine organ, influencing organ function through metabolites. The composition and abundance of metabolites within adipose tissue, encompassing fatty acids, amino acids, aldehydes, and ketones etc. [4]. Previous studies have demonstrated a favorable link between beef's taste, nutritional profile, and economic affordability, particularly attributed to its rich unsaturated fatty acid content [5, 6].Adipocytes, the fundamental structural units of adipose tissue, undergo a complex differentiation process during lipogenesis [7, 8]. However, most research efforts have been directed towards the analysis of metabolites from various tissues of bovine, the dynamic changes during the lipogenic differentiation of bovine preadipocytes and the cellular-level mechanisms regulating metabolite homeostasis remain largely unexplored [9].

In the realm of adipose tissue metabolism, significant strides have been made in understanding the complex interplay of cellular and molecular mechanisms that regulate lipid homeostasis. By comparing the transcriptomes of bovine mammary epithelial cells from high-fat and low-fat Holstein dairy cows, it has been inferred that CD44 may serve as a candidate gene affecting milk fat synthesis [10]. The KLFs family of transcription factors, which are recognized for their inhibitory role in adipogenesis, mediate their effects through intricate interactions with C/EBPs and PPARγ [11]. Despite the advancements in the field, adipogenesis remains a highly regulated and coordinated dynamic process, with the key genes orchestrating bovine fatty acid metabolism yet to be conclusively identified. The limitations of current research underscore the necessity for a more nuanced approach, advocating for a detailed and precise multi-omics analysis at the cellular level to fully unravel the intricate regulatory mechanisms governing adipose tissue metabolism in bovine.

The fatty acid desaturase gene (FADS2) encodes the rate-limiting enzyme for endogenous production of n-3 long-chain polyunsaturated fatty acids, Delta-6 desaturase (D6D). The function of the FADS2 gene is closely related to its role in lipid and glucose metabolism. Polymorphisms in the FADS2 gene and the activity of D6D are associated with insulin resistance, type 2 diabetes, dyslipidemia, and obesity [12]. Moreover, changes in the expression of the FADS2 gene and D6D activity may alter the expression of key regulators of the fatty acid desaturase pathway, thus modifying lipid content in white adipose tissue (WAT) [13]. There is also evidence suggesting that FADS2 gene polymorphisms may affect basal metabolic rate and are related to arachidonic acid metabolism in pigs [14]. Another study expanded the known saturated fatty acid substrates of FADS2 by demonstrating its involvement in the metabolism of branched-chain fatty acids (BCFA) and non-odd-chain fatty acids (n-OCFA), thereby highlighting its pivotal role in the synthesis of monounsaturated fatty acids (MUFA) [15]. It is not difficult to infer that FADS2 plays an important role in the regulation of adipogenesis, but the regulatory effect of FADS2 during the lipogenic differentiation of bovine preadipocytes is still unclear.

Altering the expression of specific genes, such as ELOVL6 and ACSL1, has been proven to significantly influence fatty acid composition in bovine adipocytes, where ELOVL6 overexpression enhances the stearic to arachidonic acid ratio and ACSL1 is crucial for unsaturated fatty acid synthesis and lipid droplet formation [6, 16]. In ruminants, adipose tissue is the primary site for fatty acid synthesis [17]. Given the established links between FADS2, D6D activity, and metabolic health, it's hypothesized that in cellular models could clarify its impact on adipocyte differentiation and lipid accumulation, particularly in livestock. As FADS2 facilitates the conversion of dietary fats into essential fatty acids, its enhanced activity might alter fatty acid profiles within adipose tissue, influencing both fat deposition and metabolic outcomes. This approach could illuminate the role of FADS2 in lipid metabolism, offering insights into improving feed efficiency and meat quality in agricultural settings.

In this study, we uncovered variations in metabolites and gene expression by integrating metabolomic and transcriptomic data during adipogenesis of bovine adipocytes on days 0, 2, 4, and 8. This study has unveiled FADS2 as a potential key player influencing the content of unsaturated fatty acids. Through the establishment of overexpression models for two distinct transcripts of FADS2 in mesenchymal stem cells, we further investigated their impact on the expression of key genes involved in adipogenesis in bovine preadipocytes and the content of fatty acids. This work sheds light on the mechanisms by which different transcripts of FADS2 regulate fatty acid metabolism and enhances our understanding of the molecular basis of fatty acid metabolism and its impact on adipose tissue development.

Results

Dataset generation for transcriptome, metabolome and lipid accumulation during bovine preadipocyte adipogenesis

In order to systematically understand the changes of gene expression, metabolites, lipid accumulation and their regulatory mechanisms in the process of fat deposition, we conducted transcriptome sequencing, evaluated the concentration of different metabolites, and recorded the changes of fat droplet size at the same time for bovine adipocytes induced on the day 0 (PA), day 2 (DA2), day 4 (DA4) and day 8 (DA8) (Fig. 1A). Generally, we saw small lipid droplets on DA2 after inducing adipocytes, and the number of small lipid droplets increased and started to get larger on DA4. On DA8, large lipid droplets were widely observed (Fig. 1A and Figure S2). We applied UPLC—MS/MS to detect 962 metabolites in 22 groups in adipocytes on PA, DA2, DA4, and DA8 (3 replicates per stage). The metabolites included acids and their metabolites, benzene and its derivatives, organic acids and their derivatives, fatty acyls (FA), nucleotides and their metabolites, glycerophosphates (GP), glycerolipids (GL), coenzymes and vitamins, and carbohydrates and their metabolites (Supplementary Fig. 2, Table S2). At the same time, the expression changes of 20,125 genes were evaluated for different adipogenesis stages using RNA-sequencing method.

Fig. 1

A Phenomenon of preadipocytes across lipid differentiation stages. Lipid droplet area accumulation across adipogenic differentiation stages. B PCA analysis of metabolomic data during lipogenic differentiation of bovine preadipocytes. C PCA analysis of the transcriptome during lipogenic differentiation of bovine preadipocytes. D Alluvial plot illustrating the distribution of DEMs across adjacent differentiation stages, with counts indicating up-/down-regulated features. E Alluvial plot illustrating the distribution of DEGs across adjacent differentiation stages, with counts indicating up-/down-regulated features. F Bubble plot of KEGG pathway enrichment analysis for DEMs and DEGs, with node sizes representing gene/metabolite counts and colors indicating significance levels

By applying PCA, cluster and correlation analysis, we evaluated the consistency among replicated samples and the quality of our datasets. Both PCA and cluster results showed that the replicated samples were clustered together (Fig. 1B, C and Supplementary Fig. 3). Intra-stage inter-replicate correlations exceeded 0.89 for both transcriptomic and metabolomic datasets at each differentiation time point (Table S3). PCA demonstrated clear temporal separation along PC1 in both omics profiles during adipocyte differentiation. Overall, we have generated effective data that can accurately exhibit changes in transcriptome and metabolome during fat deposition.

Transcriptome and metabolome variations between each stage during lipid accumulation

Both the expression levels of genes and the concentrations of metabolites showed highly diverse after the bovine preadipocyte induced to differentiation. Totally, we detected 328 DEMs (Variable Important in Projection, VIP ≥ 1) and 5,257 DEGs (False Discovery Rate, FDR < 0.05; Fold Change, FC ≥ 2 or ≤ 0.5) by comparing two different adjacent stages (Fig. 1D, E, Table S4 and Table S5). The DEMs and DEGs varied a lot among different stages. Hierarchical clustering showed stage—specific metabolite and gene expression variations (Supplementary Fig. 3). In the initial stage of adipogenesis (PA to DA2), only 51 metabolites (5.30%) enriched in FA were significantly changed (Fig. 1D and Table S4). In the middle of the adipogenesis (DA2 to DA4), largest number (259, 26.92%) of metabolites were significantly changed (Fig. 1D and Table S4). Highly active lipid synthesis was in progress within cells during this period. This involved variation of metabolites enriched in amino acid and its metabolites, nucleotide and its metabolites, organic acid and its derivatives. In the final stage of adipogenesis (DA4 to DA8), cells mainly underwent the synthesis of small lipid droplets into large lipid droplets.

As for the variation in the transcriptome, there were 2,857, 2,054, 2,761 significantly DEGs between PA to DA2, DA2 to DA4, and DA4 to DA8 (Fig. 1E and Table S5). Similar with the variation of metabolite concentrations, expression patterns of genes in different processes were largely different (Supplementary Fig. 4). Only a small percentage of DEGs were shared for the three processes (Supplementary Fig. 5). Functions of DEGs were consistent with the metabolites variations (Table S6). Expression of TP53, EZH2, RB1, etc. related to negative regulation of transcription from RNA polymerase II promoter were significantly changed in the initial stage of adipogenesis (Table S6). PPAR signaling pathway and AMPK signaling pathway related genes were significantly enriched in the middle stage and the final stage of the adipogenesis, respectively (Table S6). Moreover, we employed SUPPA to analyze alternative splicing events across transcriptome data at various time points. The results indicated that the highest number of alternative splicing events occurred at DA2, followed by DA4 and PA, with the fewest events detected at DA8 (Supplementary Fig. 6). These results show significant changes in metabolite and transcriptomic profiles during adipogenic differentiation.

KEGG enrichment analysis revealed distinct temporal patterns in metabolic and transcriptional pathways during adipogenic differentiation (Fig. 1F). Metabolomic profiling showed significant shifts in fatty acid and carbohydrate metabolism, with carbohydrate pathways underrepresented at early stages (PA-DA2) and fatty acid metabolism becoming dominant during DA2-DA4, coinciding with metabolite changes. Transcriptomic data highlighted dynamic enrichment of signaling pathways (e.g., Wnt, TNF, MAPK, PI3 K-Akt), with Wnt activation occurring despite limited DEGs in DA2-DA4, and multiple pathways activated in DA2-DA8 alongside increased DEG numbers.

KEGG enrichment analysis demonstrated the distribution of differential genes and metabolites across biological pathways (Fig. 1F). Key metabolic shifts occurred in fatty acid metabolism (DA2-DA4), with transcriptomic enrichment of Wnt signaling (DA2-DA4) and MAPK/PI3 K-Akt pathways (DA4-DA8). Crucially, fatty acid metabolic changes (metabolome) temporally aligned with upstream regulatory pathway activation (transcriptome) during DA2-DA4, suggesting a synergistic role of metabolic and transcriptional regulation in adipogenic differentiation.

Integrated analysis of DEMs and DMGs

Time-course expression analysis of differential metabolites during the adipogenic differentiation process classified them into four distinct clusters (Fig. 2A). Clusters 1, 2, and 4 exhibited an overall increasing trend and were predominantly composed of amino acids and their metabolites, as well as fatty acids. In contrast, Cluster 3 showed a decreasing trend and had a relatively low composition of amino acids and their metabolites, as well as fatty acids. Instead, Cluster 3 was enriched in nucleotides and their metabolites. To further elucidate the biological significance of these metabolite clusters, we performed KEGG enrichment analysis. The results indicated that metabolites in Clusters 1, 2, and 4 were significantly enriched in pathways such as arachidonic acid metabolism, biosynthesis of unsaturated fatty acids, fatty acid biosynthesis, and insulin resistance. In contrast, metabolites in Cluster 3 were significantly enriched in pathways related to the biosynthesis of cofactors, nucleotide metabolism, and purine metabolism (Supplementary Fig. 7). These findings highlight the distinct metabolic pathways associated with different stages of adipogenic differentiation, revealing the dynamic reprogramming of metabolic processes during this developmental transition.

Fig. 2

A Time-series variation of DEMs. Pseudotime analysis divided DEMs into four clusters, showing the composition of metabolite types within each cluster. B Module-sample correlation heatmap of WGCNA. Heatmap displays correlation patterns between co-expression modules (columns) and samples (rows). C Correlation between sample and module. Heatmap illustrating pairwise correlations among different modules. D Regulatory network of partial genes of WGCNA blue module, showing interactions between hub genes in the blue module, with edge thickness representing the weight (correlation strength) between genes

To further explore the coordinated gene expression patterns during adipogenic differentiation, we conducted Weighted Gene Co-expression Network Analysis (WGCNA) on the transcriptomic data. This analysis identified a total of 13 distinct gene modules, each representing a cluster of highly correlated genes (Fig. 2B and C). Notably, the blue module exhibited a high correlation with the DA8 group, suggesting that it may contain genes that play important roles in lipid accumulation. We constructed a regulatory network of HubGenes and other lipid—accumulation related genes within the blue module, with PC as the hub gene. PC, which participates in metabolic pathways related to glucose homeostasis, has strong co—expression relationships with genes like FADS2, SCD, ACSL4, CPT1 C, and PLIN1 (Fig. 2D). FADS2 is crucial for long-chain polyunsaturated fatty acid synthesis, and SCD encodes an enzyme for fatty acid biosynthesis. These genes are closely associated with fat—related pathways such as fatty acid production, suggesting their potential roles in lipid accumulation during the late stages of adipogenic differentiation. The constructed gene interaction network reveals potential lipid—accumulation related genes, providing a basis for further exploration of the regulatory mechanisms of lipid accumulation.

Integration analysis of DEGs and DEMs during lipid accumulation

We applied a comprehensive analysis to explore the relationship between gene expression and metabolite concentration during bovine preadipocyte lipogenic differentiation. We totally detected 5257 genes and 328 metabolites that significantly changed in at least one comparisons between two different adjacent stages. There were 9,309, 25,891 and 28,620 DEG-DEM pairs that showed a consistent changing trend, and 8,806, 19,138 and 22,233 gene- metabolite pairs that showed an opposite changing trend for the comparisons according to the nine-quadrant plot analysis (Supplementary Fig. 8, Table S7). We next performed correlation analysis for each DEG-DEM pair. The result showed that 93,543 pairs were highly correlated (|r|> 0.8 and p < 0.05) with each other for the three comparisons of different adpogeneisis stages (Fig. 3A). Among them, there were 11,691 correlation events between DEGs and DEMs that were shared in at least two comparisons.

Fig. 3

A Correlation clustering heatmap of DEGs and DEMs. Hierarchical clustering analysis showing co-variation patterns between DEGs and DEMs. Color intensity represents Pearson correlation coefficients, with red indicating positive correlations and green negative correlations. B Expression patterns of unsaturated fatty acid biosynthesis components. C CCA of differential genes and differential metabolites of unsaturated fatty acid biosynthetic pathways. When genes and metabolites are located in the same region, those farther from the origin and closer to each other indicate stronger associations

It’s of noted that most of correlated gene-metabolite pairs were saw in the final stage of adipogenesis. They were mostly enriched in the unsaturated fatty acid biosynthesis pathway including accumulation of arachidonic acid (AA), eicosapentaenoic acid (EPA), linoleic acid (LA), oleic acid, 11,14-eicosadienoic acid, and docosahexaenoic acid (DHA). This was accompanied by the temporal upregulation of FADS2, SCD, ACOT2, and ACOT7 on DA4 and DA8, contrasting with the downregulation of ACOX1, ELOVL6, HACD2, HACD4, and SCP2 (Fig. 3B). To further explore the key genes regulating unsaturated fatty acid content in the unsaturated fatty acid biosynthetic pathway, we utilized Canonical Correlation Analysis (CCA) to study the biosynthetic pathway of unsaturated fatty acids. The results indicated that FADS2, ACOT7, and ACOT2 exhibited strong associations with polyunsaturated fatty acids such as AA, EPA, LA, DHA, and FFA (20:2), as well as monounsaturated fatty acids like FFA(18:1) (Fig. 3C).

Structural and functional analysis of FADS2 alternative splicing isoforms

FADS2 displayed the highest correlation with AA, EPA, LA, and FFA (18:1) in the above results. We next focused on FADS2 to apply further research. In our previous study, two different isoforms produced by alternatively spliced events of FADS2 were identified [18]. The two FADS2 isoforms were named as FADS2-T0 and FADS2-T2, separately. The 6 th and 7 th exon of FADS2-2 were alternatively skipped, which generated an isoform lost 46 amino acids. This event was consistently observed in species from the orders Pholidota, Cetartiodactyla, and Carnivora, but was absent in some Artiodactyla and Primates (Supplementary Fig. 9). We conducted a structural and functional analysis to further explore potential different functions between the five FADS2 isoforms. By annotating the functional domain of FADS2-T0 using UniProt database, we discovered that the lost amino acid sequences generated by exon skipping events in FADS2-T2 were just located in the fatty acid desaturase domain (Fig. 4A). The three-dimensional conformation of the two FADS2 isoforms were predicted using AlphaFold3. The results showed that the amino acid sequence variation between FADS2-T0 and FADS2-T2 caused the two alpha-helices and one random coil structure differences (Fig. 4B). Further function prediction using DeepTMHMM and Deeploc2 proved that this alternative protein structure was just changed the transmembrane domain in the endoplasmic reticulum (Fig. 4C, D). These results implied that the FADS2-T2 isoform might exhibits different functions during adipogenesis comparing to the FADS2-T0 isoform.

Fig. 4

A The transcript structure and domains of FADS2 gene. B The tertiary structure diagrams of FADS2-T0 and FADS2-T2. proteins and the structure of the exon deleted in FADS2-T2 is represented in magenta. C Transmembrane region prediction results. The second transmembrane region in FADS2-T2 exhibits a shortened length compared to FADS2-T0. D Distinct features of FADS2-T0 compared to FADS2-T2 in endoplasmic reticulum. Exon deletion in FADS2-T2 induces amino acid changes within and adjacent to transmembrane domains. E Lipogenic differentiation of C3H10 T1/2 cells and oil red O staining. F qPCR for adipogenesis-related genes

We overexpressed FADS2-T0 and FADS2-T2 isoforms in C3H10 T1/2 cells separately to examine their different roles during adipogenesis. Infection efficiency assessed by Flow cytometry assessed and overexpression levels assessed RT-qPCR for FADS2-T0 and FADS2-T2 were both supported that overexpression of the two kinds of isoforms were comparable (Supplementary Fig. 10, Supplementary Fig. 11). Apoptosis assays revealed that overexpression of either FADS2 T0 or T2 primarily would lead to minimal apoptosis (Supplementary Fig. 12). Additionally, cell proliferation assays showed that both FADS2 T0 and T2 would promote cell proliferation (Supplementary Fig. 13). To evaluate the impact of the two FADS2 isoforms on adipogenesis, we induced differentiation in these cells and performed Oil Red O staining at 8 th day. The results indicated that overexpression of FADS2-T0 significantly promote the lipid formation than FADS2-T2 (Fig. 4E). We compared the expression differences of nine important adipogenesis related genes (FASN, PCNA, GPDH, MAT2 A, SCD1, PPARG, ACC1, LPL and FABP4) between overexpressed-FADS2-T0 and overexpressed-FADS2-T2 cells. The expression patterns of most genes were similar in overexpressed-FADS2-T0 and overexpressed-FADS2-T2 cells, while the expression of SCD1 and PPARG exhibited significantly different variation trends compared to the controls (Fig. 4F). The elevated SCD1 expression in FADS2-T2 cells may reflect compensatory upregulation of unsaturated fatty acid desaturases due to impaired polyunsaturated fatty acid biosynthesis capacity caused by the FADS2-T2 mutation. Concurrently, the heightened PPARG expression in FADS2-T2 cells likely corresponds to active lipid droplet formation during the DA8 metabolic phase, where enhanced lipogenic activity drives transcriptional activation of this master adipogenic regulator.

Impact of FADS2-T0 and FADS2-T2 isoform on lipid metabolism

To explore the impact of different FADS2 isoforms on lipid metabolism, we performed a comprehensive lipidomics analysis on cells overexpressing FADS2-T0 and FADS2-T2. A total of 604 lipid metabolites were identified (Fig. 5A, Table S8). PCA revealed distinct separation between the FADS2-T0 and FADS2-T2 groups, indicating potential differences in their lipid profiles (Fig. 5B). We totally identified 112 differentially produced lipid metabolites, with 109 showing higher levels in the FADS2-T2 group (Fig. 5C, D, Supplementary Fig. 14). It’s of noted that 91 (83.48%) the 109 increased lipid metabolites in FADS2-T2 belong to glycerophospholipids with different carbon chain lengths, and others were sphingolipids, and glycerolipids. Comparing to the FADS2-T0 cells, the FADS2-T2 cells downregulated Carnitine C3:1-2OH, Carnitine C4-OH, and PE(P-18:1_20:3) (Fig. 5D). KEGG enrichment analysis of the differential lipid metabolites highlighted pathways such as general metabolic pathways, glycerophospholipid metabolism, and ether lipid metabolism (Fig. 5E). We checked the FFA (free fatty acids) metabolites that had strong associations with FADS2. Notably, FFA(20:4), which is closely related to FADS2 function, was significantly higher in the FADS2 T2 group (Supplementary Fig. 15). Overall, our result proved that overexpression of FADS2-T0 and FADS2-T2 would cause different lipid metabolism processes during adipogenesis.

Fig. 5

A The composition of lipid metabolites. B PCA analysis of lipid metabolites between FADS2-T0 and FADS2-T2. C Volcano plot of differential lipid metabolites between FADS2-T0 and FADS2-T2. D Network diagram of differential lipid correlations. E KEGG enrichment analysis of metabolites between the FADS2-T0 and FADS2-T2

Discussion

In this study, we conducted a comprehensive metabolomic analysis across different stages of adipogenesis, and identified 328 differential metabolites by comparing two different adjacent stages. Temporal analysis showed increases in amino acids, fatty acids, glycerophospholipids, and glycerolipids, mirroring lipid droplet formation as observed previously [19]. KEGG enrichment analysis of these metabolites highlighted pathways like amino acid biosynthesis, general metabolism, and unsaturated fatty acid synthesis, consistent with earlier studies on adipocyte metabolic reprogramming [20]. The link between the biosynthesis of unsaturated fatty acids and adipogenesis has been well-established in previous research [21, 22]. Our findings underscore the intricate metabolic reconfigurations that underpin adipocyte differentiation and highlight the importance of lipid metabolism in this process. However, metabolomics alone cannot provide a comprehensive analysis of the regulatory mechanisms underlying adipose deposition in bovine preadipocytes, additional methodologies are required to further elucidate this process.

Integrating transcriptomic and metabolomic data, as previous research suggests, is crucial for understanding adipocyte differentiation [23]. In our study, transcriptome sequencing of bovine preadipocytes identified 20,125 genes. 5,257 DEGs were enriched in the unsaturated fatty acid biosynthesis pathway, consistent with metabolome results and previous studies [24, 25]. Canonical Correlation Analysis of this pathway showed strong associations between FADS2, ACOT7, ACOT2 and polyunsaturated fatty acids. The identification of critical genes and their correlation with metabolites offers a foundation for future research aimed at enhancing adipose tissue development and lipid metabolism in livestock.

FADS2 is recognized as the rate-limiting factor in mediating unsaturated fatty acid biosynthesis [26,27,28]. Research has shown that FADS2 gene polymorphisms influence the composition of fatty acids [29, 30]. To explore the critical function of the FADS2 gene in the intricate metabolic fluxes associated with adipogenic differentiation in bovine preadipocytes, we have established a C3H10 T1/2 cell line with stable FADS2 overexpression, thereby creating a robust in vitro model for this purpose. Our research confirms that FADS2 overexpression markedly upregulates key lipogenic genes essential to lipid metabolism, such as ACC1, FABP4, FASN, GPDH, LPL, PPARG, SCD1, and HACD2, findings consistent with existing literature [24, 25, 31].

Our study further investigated the impact of different FADS2 isoforms on lipid metabolism. The FADS2-T0 exhibited stronger lipid droplet accumulation in the later stages of adipogenesis, while the FADS2-T2 showed higher levels of metabolites in glycerophospholipids, sphingolipids, and glycerolipids. These observations were linked to the interactions of FADS2 with related proteins. The protein–protein interaction (PPI) network analysis via the STRING database (confidence score ≥ 0.7) identified FADS2 interactions with elongases (ELOVL5, ELOVL6), phospholipases (PLA2G3, PLA2G16), and stearoyl-CoA desaturase (SCD), suggesting a multifaceted regulatory axis in bovine lipid metabolism (Supplementary Fig. 16).The ELOVL family enzymes are well-documented to catalyze fatty acid elongation, particularly for polyunsaturated fatty acids (PUFAs) [32]. PLA2G family members, such as PLA2G3 and PLA2G16, are involved in phospholipid hydrolysis and membrane remodeling [33]. Fatty acid β-oxidation is mediated by carnitine, with acyl-CoAs serving as substrates for this process [34]. The expression of SCD1 were detected differentially expressed when overexpressed FADS2-T0 and FADS2-T2 in our study. The elevated SCD1 expression in FADS2-T2 cells may reflect compensatory upregulation of unsaturated fatty acid desaturases due to impaired polyunsaturated fatty acid biosynthesis capacity caused by the FADS2-T2 mutation. We propose that the observed changes in lipid metabolites could be related to the presence of different FADS2 splice variants.

Conclusion

This study firstly constructed a network for dynamic integration between gene expression and metabolite contents during bovine adipogenesis. Multi-omics identified FADS2 as a key regulator of UFA biosynthesis, with overexpression experiments confirming its role in activating lipogenic genes. Study further uncovered that distinct FADS2 isoforms would drive divergent lipid metabolic outcomes, which might be responsible for different adipogenesis regulation through the variation of a 46-amino acid sequence. These findings advance mechanistic insights into lipid metabolism regulation, offering targets for optimizing livestock adipose development.

Methods

Cell culture and processing

The bovine preadipocytes from subcutaneous adipose tissue were prepared as following. Briefly, all excised adipose tissues were cut into about 1 mm3 small sections, which were then digested in DMEM containing 1 mg/ml Type I collagenase for 1 h at 37 °C with shaking at 15- cycles/min. The cell suspension was sequentially filtered through a 250 μm sterilized nylon mesh. The filtrate was centrifuged at 1000X g for 5 min. The pelleted cells with fibroblast morphology were resuspended in growth medium. Bovine preadipocytes and C3H10 T1/2 cells were used to perform different experiments in our study. Both of them were cultured in a high glucose DMEM medium (4.5 g/mL Glucose, 4.0 mM L-glutamine, Cytiva) containing 10% FBS supplemented with 1% penicillin–streptomycin at 37 °C and 5% CO2. The culture medium was changed every 48 h. Cells were passaged at 80% density in medium dishes and passed to the fourth generation for induced differentiation. After cells reached confluence in 6-well plates, a medium containing 1.0 μmol/L dexamethasone (DEX), 0.5 mmol/L 3-isobutyl-1-methylxanthine (IBMX), 1.0 μmol/L rosiglitazone, and 10 mg/L insulin were used to induce cells for 48 h. Then, fresh medium containing 10 mg/L insulin and 1.0 μmol/L rosiglitazone was used to maintain the differentiation [35]. Cells were washed twice using phosphate buffer saline (PBS) before collection. Cells were collected using a cell scraper (LABSELECT) at the 0 th, 2nd, 4 th, and 8 th days of induced differentiation, snap frozen in liquid nitrogen for 20 min, and stored in a −80°C freezer, with three replicates for each period.

Sample preparation and LC–MS analysis

The sample stored at liquid nitrogen was thawed on ice. A 500 μL solution (Methanol:Water = 4:1, V/V) containing internal standard (L-2-chlorophenylalanine, [2H5]-Hippuric Acid and [2H5]-Phenoxy acetic Acid) was added into the cell sample and vortexed for 3 min [36]. The sample was placed in liquid nitrogen for 5 min and on dry ice for 5 min, and then thawed on ice and vortexed for 2 min. This freeze–thaw circle was repeated three times in total. The sample was centrifuged at 12,000 rpm for 10 min (4 °C). A 300 μL of supernatant was collected and placed at −20°C for 30 min. The sample was then centrifuged at 12,000 rpm for 3 min (4 °C). A 200 μL aliquots of supernatant were transferred for the UPLC-MS/MS analysis.

T3 UPLC Conditions: The sample extracts were analyzed using an LC–ESI–MS/MS system (UPLC, ExionLC AD, https://www.sciex.com.cn/; MS, QTRAP® System, https://www.sciex.com/). The analytical conditions were as follows, UPLC: column, Waters ACQUITY UPLC HSS T3 C18 (1.8μm, 2.1 mm*100 mm); column temperature, 40 °C; flow rate, 0.4 mL/min; injection volume, 2μL or 5μL; solvent system, water (0.1% formic acid): acetonitrile (0.1% formic acid); gradient program, 95:5 V/V at 0 min, 10:90 V/V at 10.0 min, 10:90 V/V at 11.0 min, 95:5 V/V at 11.1 min, 95:5 V/V at 14.0 min.

The Triple TOF mass spectrometer was used for its ability to acquire MS/MS spectra on an information-dependent basis (IDA) during an LC/MS experiment. In this mode, the acquisition software (TripleTOF 6600, AB SCIEX) continuously evaluates the full scan survey MS data as it collects and triggers the acquisition of MS/MS spectra depending on preselected criteria. In each cycle, 12 precursor ions whose intensity is greater than 100 were chosen for fragmentation at collision energy (CE) of 30 V (12 MS/MS events with product ion accumulation time of 50 ms each). ESI source conditions were set as follows: Ion source gas 1 as 50 Psi, Ion source gas 2 as 50 Psi, Curtain gas as 25 Psi, source temperature 500 °C, Ion Spray Voltage Floating (ISVF) 5500 V or −4500 V in positive or negative modes, respectively.

LIT and triple quadrupole (QQQ) scans were acquired on a triple quadrupole-linear ion trap mass spectrometer (QTRAP), QTRAP® LC–MS/MS System, equipped with an ESI Turbo Ion-Spray interface, operating in positive and negative ion modes and controlled by Analyst 1.6.3 software (Sciex). The ESI source operation parameters were as follows: source temperature 500 °C; ion spray voltage (IS) 5500 V (positive), −4500 V (negative); ion source gas I (GSI), gas II (GSII), curtain gas (CUR) were set at 50, 50, and 25.0 psi, respectively; the collision gas (CAD) was high. Instrument tuning and mass calibration were performed with 10 and 100 μmol/L polypropylene glycol solutions in QQQ and LIT modes, respectively. A specific set of MRM transitions was monitored for each period according to the metabolites eluted within this period.

Quantification of lipid droplets

To investigate the cellular changes at different time points, nuclei were stained with DAPI (Solarbio, Beijing, China), while lipid droplets were labeled with BODIPY 493/503 (MedChemExpress, New Jersey, USA). Subsequently, the stained cells were imaged using the Opera Phenix high-content confocal imaging system (PerkinElmer, Massachusetts, USA). The acquired images were processed with the Harmony software, and phenotypic data related to lipid droplets were analyzed using ImageJ.

Oil Red O staining was performed using the Oil Red O staining kit (Solarbio, Beijing, China). Images of Oil Red O-stained adipocytes were acquired using a Nikon imaging system (Nikon, TE2000-U, Tokyo, Japan) at 10 × and 20 × magnification for cells induced.

RNA extraction and qPCR assay

Total RNA was extracted from cultured the fourth generation bovine preadipocytes using the FastPure Cell/Tissue Total RNA Isolation Kit V2 (Vazyme, Nanjing, China). Complementary DNA (cDNA) was reverse-transcribed from total RNA (1 µg) (Vazyme, China). Quantitative analysis of gene expression was performed by qPCR according to the SYBR Green I chimeric fluorescence kit (Vazyme, China) on the CFX Connect Real-Time PCR platform (BIO-RAD). Primers used in this study are shown in Table S1. The qPCR conditions were as follows: 5 min at 95 °C, 10 min at 95 °C, 40 cycles, 10 s at 95 °C, 10 s at 60 °C, and 15 s at 72 °C. The mRNA expression was normalized by comparison with the bovine cytoskeleton actin gene ACTB (NM_173979.3). The 2^−ΔΔCT method was employed to analyze the qPCR result.

RNA sequencing

Total RNA from the induced differentiated cells in the 0 th, 2nd, 4 th, and 8 th days were extracted using TRIzol (Invitrogen, Carlsbad, CA, USA) according to the Manufacturer’s instructions. We measured the quantity and purity of RNA using a NanoDrop 8000 Spectrophotometer (NanoDrop Technologies, Wilmington, DE, USA) and Agilent 2100 Bioanalyzer System (Agilent Technologies, Santa Clara, CA, USA). Libraries were constructed and sequenced using the DNBSeq-T7 platform with paired-end (150 bp) reads.

NGSQCToolkit v2.3.3 was used to trim adapter sequences, and low-quality reads [37]. The clean reads were aligned on the reference genome (ARS-UCD1.2) along with gene annotation from the Ensembl website (https://www.ensembl.org/index.html, June 25, 2021) using the HISAT2 v2.1.0 with the default parameters [38]. The spliced reads were assembled into transcripts using the StringTie v1.3.3 software. DEG were measured using DESeq2 v3.15 [39]. Gene functional annotation analyses were applied using the online DAVID software (https://david.ncifcrf.gov/, July 30, 2021) [40]. Gene enrichment in annotation terms was assessed using Fisher’s exact test. The p-values were subjected to FDR correction for multiple comparisons. We analyzed alternative splicing events including alternative 3’ splice site (A3), alternative 5’ splice site (A5), alternative first exon (AF), alternative last exon (AL), mutually exclusive exons (MX), retained intron (RI), and skipping exon (SE) using SUPPA [41, 42]. Local splicing events were generated from Ensembl annotations (ARS-UCD1.2), and transcript abundance estimates (TPM) from StringTie were processed to compute PSI values. Differential splicing analysis across differentiation time points was performed, using default parameters.

DEMs and DEGs

VIP values were extracted from the OPLS-DA results, which also contained score plots and permutation plots generated using the R package MetaboAnalystR. The data were log-transformed (log2) and mean-centered before OPLS-DA. To avoid overfitting, a permutation test (200 permutations) was performed. DEMs between groups were filtered by VIP ≥ 1 and absolute Log2 FC ≥ 1. Significant DEGs between groups were filtered by absolute Log2 FC ≥ 1 and FDR < 0.05.

Lentiviral packaging and infection

The pCDH-CMV-MCS-EF1-copGFP-T2 A-Puro was selected as the lentiviral core plasmid, and pMD2.G and psPAX2 were selected as the packaging helper plasmid for the three-plasmid expression system for lentiviral packaging. pCDH-CMV-MCS-EF1-copGFP-T2 A-Puro has all accessory genes of the packaging plasmid except for vif, vpu, vpr, and nef, which increased the safety of the vector without affecting the titer and transfection ability of the virus. Among them, pCDH-CMV-MCS-EF1-copGFP-T2 A-Puro contains a separate EF1 start sequence to activate the expression of green fluorescent protein (COPGFP), and it carries puromycin (PURO) resistance gene and ampicillin (Ampicillin) resistance gene, offering flexibility and options during the screening process. The 293 T cell was inoculated in a six-well cell culture plate. Three different plasmids (pMD2.G, psPAX2, core plasmid) were mixed according to the ratio of 2:3:4. The mixed plasmids were then introduced into the 293 T cells to package lentivirus using the jetPRIME®. The lentivirus was collected at 48 h and used to infect C3H10 T1/2 cells.

Protein functional prediction

To predict the functions of proteins translated from transcripts, we utilized AlphaFold for three-dimensional structure prediction [43]. We employed the UniProt database to analyze the structural domains of proteins encoded by FADS2 transcripts [44]. Additionally, we utilized DeepTMHMM for predicting protein transmembrane regions and Deeploc2 for determining the subcellular localization of these proteins [45, 46]. We collected translated protein sequences from various species in the NCBI database, performed multiple sequence alignments using the muscle software, and visualized the results with ENDscript2 [47, 48]. To explore potential protein–protein interactions, we employed the STRING database [49].