Transcriptomics and Lipid Metabolomics Analysis of Subcutaneous, Visceral, and Abdominal Adipose Tissues of Beef Cattle

Fat deposition traits are influenced by genetics and environment, which affect meat quality, growth rate, and energy metabolism of domestic animals. However, at present, the molecular mechanism of fat deposition is not entirely understood in beef cattle. Therefore, the current study conducted transcriptomics and lipid metabolomics analysis of subcutaneous, visceral, and abdominal adipose tissue (SAT, VAT, and AAT) of Huaxi cattle to investigate the differences among these adipose tissues and systematically explore how candidate genes interact with metabolites to affect fat deposition. These results demonstrated that compared with SAT, the gene expression patterns and metabolite contents of VAT and AAT were more consistent. Particularly, SCD expression, monounsaturated fatty acid (MUFA) and triglyceride (TG) content were higher in SAT, whereas PCK1 expression and the contents of saturated fatty acid (SFA), diacylglycerol (DG), and lysoglycerophosphocholine (LPC) were higher in VAT. Notably, in contrast to PCK1, 10 candidates including SCD, ELOVL6, ACACA, and FABP7 were identified to affect fat deposition through positively regulating MUFA and TG, and negatively regulating SFA, DG, and LPC. These findings uncovered novel gene resources and offered a theoretical basis for future investigation of fat deposition in beef cattle.


Introduction
Adipose tissue, known as a metabolically heterogeneous endocrine organ mostly composed of adipocytes [1], performs a crucial role in lipid metabolism, energy regulation, temperature maintenance, insulin resistance, and meat quality of farm animals [2,3]. In mammals, adipose tissue is temporal and tissue-specific [4,5]. Anatomically, adipose tissue could be classified subcutaneous adipose tissue (SAT), visceral adipose tissue (VAT), abdominal adipose tissue (AAT), intramuscular adipose tissue, intermuscular adipose tissue, heart adipose tissue, and kidney adipose tissue. These adipose tissues displayed cellular, molecular, physiological, and functional differences [6]. Notably, the larger proportion of small adipocytes in SAT resulted in its higher insulin sensitivity and a sharper absorption of circulating fatty acids (FAs) and triglycerides (TGs), which could prevent fat deposition in non-adipose tissue [7]. Studies reported higher concentrations of SAT contributed to being protective against metabolic diseases, conversely, VAT accumulation was prone to metabolic syndrome [8]. Additionally, SAT could protect the carcass from cold shortening and drip loss during cooling [9], whereas excessive abdominal fat was considered to be economic loss [10].
Currently, the RNA sequencing (RNA-seq) technique has become a powerful tool for quantitative analysis of the whole transcriptome expression profile, which is successfully applied to the complex traits of farm animals like fat deposition and fatty acid traits of beef of SAT contributed to being protective against metabolic diseases, conversely, VAT accumulation was prone to metabolic syndrome [8]. Additionally, SAT could protect the carcass from cold shortening and drip loss during cooling [9], whereas excessive abdominal fat was considered to be economic loss [10].
Currently, the RNA sequencing (RNA-seq) technique has become a powerful tool for quantitative analysis of the whole transcriptome expression profile, which is successfully applied to the complex traits of farm animals like fat deposition and fatty acid traits of beef cattle [11,12]. The differential expression level of lipogenic genes had been proven to have an effect on lipid metabolism in bovine adipocytes [13], such as SCD, FABP4, and LPL in Japanese Black and Holstein cattle [14]. As the ultimate omics level of biological systems, metabolomics has a distinct advantage in measuring the final products of genetic and environmental interactions, and then screens significant metabolic molecules as biomarkers to reflect the interrelation between the changes of biomarkers and phenotypic traits [15], which has been widely used in biomedicine [16], food nutrition [17], and crop trait selection, breeding, and evaluation [18,19]. Notably, a relevant report emphasized lipidomics was a secondary omics rather than a discipline separate from metabolomics as the units it studied were metabolites [20]. In comparison with other omics, lipidomics was still undergoing development as its relevant tools, databases, and informatics resources were not yet perfect. Previous work showed lipidomics had been applied to elucidate the lipid profile changes in brown and white adipose tissue [21,22]. Additionally, several studies reported the marker metabolites of intramuscular fat in beef cattle, such as capric acid and 3-phosphoglyceric acid in Japanese black cattle [23], oxysterol in Angus, Hereford and Wagyu  Angus [24], and glucose, choline, and acetyl groups in wagyu cross steers [25]. However, there were no reports on the metabolomics of SAT, VAT, and AAT in beef cattle.
With the development of single omics technology, multi-omics research has become a hotspot to systematically explore the genetic regulatory mechanism of the quantitative traits of livestock, which could improve the comprehensibility and reliability of the results and facilitate the progress of animal breeding. The integrated transcriptomics and metabolomics research on grass-fed and grain-fed Angus cattle reported that the meat produced by grass-fed beef was more tender and had a higher omega3/omega6 ratio than grain-fed beef [26]. Another transcriptomic and metabolomic study of adipose tissue in goats showed that L-carnitine could promote the differentiation and thermogenesis of cells in brown adipose tissue via activating the AMPK signaling pathway while reducing lipid accumulation by inducing lipolysis and thermogenesis of cells in white adipose tissue [27]. Here, this study collected three types of adipose tissues (SAT, VAT, and AAT) to identify candidate genes and key metabolites for fat deposition in beef cattle via a comprehensive analysis of transcriptomic and lipid metabolomics (Figure 1), which not only offered a data basis for further understanding the characteristics and differences of different adipose tissues of beef cattle, but provided novel insights for the gene mapping strategy in other livestock.

Ethics Statement
All the experimental procedures relating to animals were authorized by the Ethics Committee of the Institute of Animal Sciences, Chinese Academy of Agricultural Sciences (IAS, CAAS), Beijing, China (approval number: IAS2020-48).

Animals and Sampling
The 73 experimental individuals slaughtered in this study were all from Huaxi cattle resource population established in Inner Mongolia Aokesi Livestock Breeding Co., Ltd. (Inner Mongolia, China) After weaning, the cattle were raised under the same fattening strategies (Feed composition: corn silage, corn straw, wine tank, soybean meal, corn tablet, concentrated feed, and concentrate supplement) and unified standardized management according to GB/T 27642-2011, and then slaughtered in Zhongao Food Co., Ltd. (Aohan Banner, Chifeng, China). The SAT, VAT, and AAT of six healthy bulls with an average age of 26 months and an average weight of 700 kg were immediately collected and frozen in liquid nitrogen for the extraction of total RNA, fatty acid, and lipid, respectively.

Functional Enrichment Analysis
ClusterProfiler (v4.0.0) package was used to implement Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses of DEGs and significant modules obtained by WGCNA, respectively [33]. GO terms and pathways with p < 0.05 and q < 0.05 were considered to be significantly enriched. Particularly, the PPAR signaling pathway and biological processes of lipid metabolism were worthy of great concern [34]. The protein-protein interaction (PPI) network of candidate genes was constructed via the STRING database (https://string-db.org/ (accessed on 4 February 2022)) to explore the interaction relationship between DEGs.

Fatty Acid Species and Content Determination
The fatty acid in each sample was extracted in the Central Laboratory Platform of the IAS of CAAS following their standard procedures. The fatty acid content (X) was expressed as grams per 100 g (g/100 g) and its calculation formula was as follows: X i : the content of each fatty acid in each sample; A i : the peak area of each fatty acid methyl ester in each sample; C si : the concentration of fatty acid methyl ester standard in the standard determination fluid in milligrams per kilogram (mg/kg); C i : the concentration of each fatty acid methyl ester in the sample liquid calculated from the standard curve, in mg/kg; A si -the peak area of each fatty acid methyl ester standard contained in the standard determination fluid; f : the proportion of each fatty acid in the total standard fluid; F FAME−FA : the conversion coefficient of each fatty acid methyl ester to fatty acid; d: the internal standard factor; 0.0001: the coefficient of fatty acid content per gram converted to fatty acid content per 100 g of sample.

Lipid Species and Content Determination
The lipid components in SAT, VAT, and AAT were extracted by a three-phase (MTBEmethanol-water) solvent system combined with tissue homogenization technology. Lipids were detected using ultra-high performance liquid chromatography-high resolution mass spectrometry (UPLC-HRMS) via ultimate TM 3000 ultra-performance liquid chromatographytandem quadrupole-electrostatic field orbital trap high-resolution mass spectrometry (Thermo Scientific, Waltham, MA, USA). An Accucore C30 core-shell column was utilized for lipid molecules' chromatographic separation under 50 • C which were eluted with 60% acetonitrile in water (A) and 10% acetonitrile in isopropanol (B) both containing 10 mM ammonium formate and 0.1% formate. Positive and negative ion electrospray ionization was used to detect mass spectrometry (MS). The same amount of lipid extract from each sample was taken to be mixed and centrifuged to prepare the pooled quality control (QC) sample for data quality control in the untargeted lipidomics testing. In this study, four QC samples were evenly inserted into the analysis sequence. The sample extract, lipid identification, and quantification analyses were conducted at Beijing BioMiao Biotechnology Co., Ltd. (Beijing, China).

Lipidomics Raw Data Preprocessing
LipidSearch software (v4.1) was conducted to process untargeted lipidomics data including peak picking and lipid identification. The detected MS2 spectrum was searched against in silico predicted spectra of a diverse phospholipid, neutral glycerolipid, sphingolipid, glycosphingolipids, steroids, etc. The mass accuracy for precursor and MS/MS product ions searching were set at 5 ppm and 5 mDa, respectively. TraceFinder software (v4.1; Thermo Scientific, Waltham, MA, USA) was used to extract the area under curve (AUC) of all lipid molecules for lipid quantitative, and then the AUC data were transferred to Excel software (Microsoft, Redmond, WA, USA) for data normalization using a linear regression algorithm.

Lipid Metabolomic Profiling Analysis
The lipidomic data deriving from different measurements were merged and those detected with multiple methods were excluded to ensure the uniqueness of lipid, and then Log2 was transformed for final statistical analysis. SIMCA (v14.1; Umetrics, Umea, Sweden) and OmicShare tools (https://www.omicshare.com/tools/ (accessed on 24 February 2022)) were employed to process orthogonal partial least squares discriminant analysis (OPLS-DA) on three adipose tissues. The differently accumulated lipids (DALs) were selected with variable importance in projection (VIP) > 1 and the independent sample t-test p-value < 0.05 [35]. Lastly, using MetaboAnalyst (https://www.metaboanalyst.ca (accessed on 26 February 2022)) for lipidomic pathway analysis of DALs [36].

Transcriptomic and Lipid Metabolomic Integration Analysis
The Spearman correlation coefficients (SPCC) were calculated for the integrative analysis between the FPKM expression of candidate genes and the content of metabolites including important fatty acids and DALs. Gene-metabolite correlation with the criteria |SPCC| > 0.5 and p < 0.05 will be selected for further investigation.

The Evaluation of Sequencing Data
In the present study, 14 samples (five SAT, four VAT, and five AAT) were finally selected for transcriptome analysis according to the sequencing quality assessment of adipose tissue. After quality control (Table S1), 93.91G clean sequence was obtained. The average values of Q20 and Q30 were 97.81% and 94.00%, respectively, and the average GC content was 52.19%. The alignment of clean data to the BTA reference genome (ARS-UCD1.2) showed that the total mapping rate of each sample was over 93%, further indicating the quality of clean data was high and could be used for transcriptome analysis.

Principal Component Analysis among Three Groups
Principal component analysis was conducted on the gene expression level of samples in SAT-VAT ( Figure S1a), SAT-AAT ( Figure S1b), and VAT-AAT ( Figure S1c) groups. In the three groups, the variation of gene expression explained by PC1 and PC2 was more than 50% and 15%, respectively, illustrating the samples were grouped reasonably.

Identification of Differentially Expressed Genes among Three Groups
With the classical criteria p-adj < 0.05 and |log 2 FoldChange| ≥ 1, a total of 3642, 1926, and 1061 DEGs were screened, respectively, in SAT-VAT ( Figure 2a and Table S2), SAT-AAT ( Figure 2b and Table S3), and VAT-AAT ( Figure 2c and Table S4) groups. Hierarchical clustering analysis was performed on 4370 DEGs identified from the three groups ( Figure 3a). The result showed the expression patterns of genes were similar within the three groups. Notably, compared with SAT, the gene expression patterns of VAT and AAT were more consistent. Additionally, 124 overlapped DEGs were identified in the three groups, including 119 coding genes, three lncRNA, and two pseudogenes, among which 110 annotated coding DEGs located in bovine autosomes were mainly focused on ( Figure 3b and Table S5).

Transcriptomic and Lipid Metabolomic Integration Analysis
The Spearman correlation coefficients (SPCC) were calculated for the integrative analysis between the FPKM expression of candidate genes and the content of metabolites including important fatty acids and DALs. Gene-metabolite correlation with the criteria |SPCC| > 0.5 and p < 0.05 will be selected for further investigation.

The Evaluation of Sequencing Data
In the present study, 14 samples (five SAT, four VAT, and five AAT) were finally selected for transcriptome analysis according to the sequencing quality assessment of adipose tissue. After quality control (Table S1), 93.91G clean sequence was obtained. The average values of Q20 and Q30 were 97.81% and 94.00%, respectively, and the average GC content was 52.19%. The alignment of clean data to the BTA reference genome (ARS-UCD1.2) showed that the total mapping rate of each sample was over 93%, further indicating the quality of clean data was high and could be used for transcriptome analysis.

Principal Component Analysis among Three Groups
Principal component analysis was conducted on the gene expression level of samples in SAT-VAT ( Figure S1a), SAT-AAT ( Figure S1b), and VAT-AAT ( Figure S1c) groups. In the three groups, the variation of gene expression explained by PC1 and PC2 was more than 50% and 15%, respectively, illustrating the samples were grouped reasonably.

Identification of Differentially Expressed Genes among Three Groups
With the classical criteria p-adj < 0.05 and |log2FoldChange| ≥ 1, a total of 3642, 1926, and 1061 DEGs were screened, respectively, in SAT-VAT ( Figure 2a and Table S2), SAT-AAT ( Figure 2b and Table S3), and VAT-AAT ( Figure 2c and Table S4) groups. Hierarchical clustering analysis was performed on 4370 DEGs identified from the three groups ( Figure 3a). The result showed the expression patterns of genes were similar within the three groups. Notably, compared with SAT, the gene expression patterns of VAT and AAT were more consistent. Additionally, 124 overlapped DEGs were identified in the three groups, including 119 coding genes, three lncRNA, and two pseudogenes, among which 110 annotated coding DEGs located in bovine autosomes were mainly focused on ( Figure 3b and Table S5).

Functional Enrichment Analysis of Differentially Expressed Genes
GO annotation performed on the DEGs identified in the SAT-VAT, SAT-AAT, and VAT-AAT groups demonstrated these DEGs were significantly enriched in 194, 133, and 20 GO terms, respectively (Table S6). The DEGs screened in the SAT-VAT group were associated with the organic acid metabolic process (GO:0006082), oxoacid metabolic process (GO:0043436), and carboxylic acid metabolic process (GO:0019752). For the SAT-AAT group, the DEGs were mainly involved in vasculature and blood vessel development (GO:0001944, GO:0001568), wound healing (GO:0042060), and cardiovascular system development (GO:0072358). Cell components and biological processes such as cell adhesion, movement, and migration were primarily enriched by the DEGs identified in the VAT-AAT group. Additionally, Figure 4a these DEGs screened in the three groups. The details of pathway enrichment were listed in Table S7. A total of 14 overlapped pathways were enriched in the three groups ( Figure 4d and Table 1). Notably, the PPAR signaling pathway played an important role in fat deposition. Therefore, in the three groups, three overlapped DEGs (PCK1, PLIN4, and FABP7) enriched in this pathway should be focused on.

Functional Enrichment Analysis of Differentially Expressed Genes
GO annotation performed on the DEGs identified in the SAT-VAT, SAT-AAT, and VAT-AAT groups demonstrated these DEGs were significantly enriched in 194, 133, and 20 GO terms, respectively (Table S6). The DEGs screened in the SAT-VAT group were associated with the organic acid metabolic process (GO:0006082), oxoacid metabolic process (GO:0043436), and carboxylic acid metabolic process (GO:0019752). For the SAT-AAT group, the DEGs were mainly involved in vasculature and blood vessel development (GO:0001944, GO:0001568), wound healing (GO:0042060), and cardiovascular system development (GO:0072358). Cell components and biological processes such as cell adhesion, movement, and migration were primarily enriched by the DEGs identified in the VAT-AAT group. Additionally, Figure Table S7. A total of 14 overlapped pathways were enriched in the three groups ( Figure 4d and Table 1). Notably, the PPAR signaling pathway played an important role in fat deposition. Therefore, in the three groups, three overlapped DEGs (PCK1, PLIN4, and FABP7) enriched in this pathway should be focused on.    The hierarchical clustering on the FPKM expression levels of 4370 DEGs identified in three groups showed SAT clustered separately, whereas VAT and AAT clustered more mixed, further implying that VAT and AAT had more similar structures or functions ( Figure 5a). Therefore, VAT and AAT were merged into one group and renamed "VAAT" in this study. The "pickSoftThreshold" function was then used to calculate the β value. Finally, an appropriate value, β = 14, was determined for the scale-free network construction (Figure 5b). The dynamic-cutting method was used for functional module division, and a total of 12 functional modules were obtained (Figure 5c), six of which were significant modules, including darkolivegreen (|r| = 0.81, p = 4 × 10 −4 ), mediumpurple3 (|r| = 0.67, p = 0.009), purple (|r| = 0.82, p = 3 × 10 -4 ), darkmagenta (|r| = 0.65, p = 0.01), floralwhite (|r| = 0.79, p = 9 × 10 -4 ), and lightcyan (|r| = 0.63, p = 0.02) (Figure 5d and Table S8). Except for the darkolivegreen module that had a strong negative correlation with SAT, the other five modules were significantly positively correlated with SAT and negatively correlated with VAAT. Of note, the genes in the grey module did not belong to any module. mixed, further implying that VAT and AAT had more similar structures or functions (Figure 5a). Therefore, VAT and AAT were merged into one group and renamed "VAAT" in this study. The "pickSoftThreshold" function was then used to calculate the β value. Finally, an appropriate value, β = 14, was determined for the scale-free network construction ( Figure 5b). The dynamic-cutting method was used for functional module division, and a total of 12 functional modules were obtained (  Table S8). Except for the darkolivegreen module that had a strong negative correlation with SAT, the other five modules were significantly positively correlated with SAT and negatively correlated with VAAT. Of note, the genes in the grey module did not belong to any module.

Significant Module Genes Analysis
Functional enrichment analysis for the six significant module genes was firstly conducted with the "clusterProfile" package (Table S9). As illustrated in Figure 6a,b, genes in the lightcyan module were mainly involved in biological processes related to fat metabolism, including the fatty acid metabolic process (GO:0006631), organic acid metabolic process (GO:0006082), PPAR signaling pathway (bta03320), fatty acid metabolism (bta01212), etc. There were no significant GO terms and pathways enriched by the darkmagenta and medumpurple3 module genes, respectively. Genes in the darkolivegreen, purpl, and floralwhite modules were mainly enriched to pathways associated with signal transduction

Significant Module Genes Analysis
Functional enrichment analysis for the six significant module genes was firstly conducted with the "clusterProfile" package (Table S9). As illustrated in Figure 6a,b, genes in the lightcyan module were mainly involved in biological processes related to fat metabolism, including the fatty acid metabolic process (GO:0006631), organic acid metabolic process (GO:0006082), PPAR signaling pathway (bta03320), fatty acid metabolism (bta01212), etc. There were no significant GO terms and pathways enriched by the darkmagenta and medumpurple3 module genes, respectively. Genes in the darkolivegreen, purpl, and floralwhite modules were mainly enriched to pathways associated with signal transduction and disease occurrence. Taken together, the module functional enrichment analysis showed that the genes in the lightcyan module were mainly related to fat deposition; thus, this module was selected for the subsequent analysis. Table S10 listed the significantly enriched GO terms and pathways regarding fat deposition and their involved genes in the lightcyan module. Deleting duplicated genes, a total of 33 DEGs were determined in this module. Furthermore, according to the filter criteria with |gene significance (GS)| > 0.4 and |module membership (MM)| > 0.8, 90 hub genes were screened in this module. Combined with the above two conditions, a total of 17 genes such as SCD, ELOVL6, ACACA, SLC27A6, ADIPOQ, OXCT1, PCK2, ACADSB, and EHHADH were selected in the lightcyan module ( Table 2). Among them, RGN and PDHA1 would not be considered as key candidate genes affecting fat deposition as they were located on the sex chromosome. Therefore, 15 genes were initially screened as important candidates for further analysis. Combined with the above two conditions, a total of 17 genes such as SCD, ELOVL6, ACACA, SLC27A6, ADIPOQ, OXCT1, PCK2, ACADSB, and EHHADH were selected in the lightcyan module (Table 2). Among them, RGN and PDHA1 would not be considered as key candidate genes affecting fat deposition as they were located on the sex chromosome. Therefore, 15 genes were initially screened as important candidates for further analysis.   Figure 6. GO annotation (a) and KEGG pathway (b) of genes enriched in significant modules.

Identification of Candidate Genes Affecting Fat Deposition
Combined with DEGs and WGCNA analysis strategies among three adipose tissues, 18 potential candidate genes were screened to affect fat deposition, and their average expression levels in three adipose tissues were analyzed (Table S11). The results showed that SCD was abundantly expressed in SAT, whereas PCK1 expression level was the highest in AAT. TDH and SHC4 could be approximated as non-expressed genes as their expression levels in the three tissues were less than or close to 1. Therefore, through gene expression, literature investigation, and gene function, 15 genes including SCD, ELOVL6, ACACA, OXCT1, ACADSB, EHHADH, FABP7, PLIN4, PCK1, PCK2, MBOAT7, PM20D1, CPT2, CHPT1, and PRKAR2A were selected to construct PPI network for further determining the pivotal candidate genes ( Figure S2). The thicker the interaction line, the higher the interaction score between the two genes. Those with confidence greater than 0.4 were regarded as moderately important interaction networks, and those with confidence greater than 0.9 were defined as important interaction networks [37]. There was a stronger interaction between ACACA and ELOVL6 as the interaction confidence was as high as 0.838 (Table 3). PLIN4, MBOAT7, PM20D1, and PRKAR2A had no interaction with other genes. Therefore, 11 genes including SCD, ELOVL6, ACACA, OXCT1, ACADSB, EHHADH, FABP7, PCK1, PCK2, CPT2, and CHPT1 were finally filtrated as candidates affecting fat deposition (Figure 7).  In this study, twenty-four fatty acid species were detected in all three adipose tissues (Table S12), including 11, 5, and 8 kinds of saturated fatty acids (SFA), monounsaturated fatty acids (MUFA), and polyunsaturated fatty acids (PUFA), respectively, which correspondingly accounted for 59.49-76.91%, 20.34-37.94%, and 2.52-2.79% of the total fatty acid content. The fatty acid content analysis represented that there was no significant difference in SFA content between VAT and AAT; however, both of them were significantly higher than those in SAT (Figure 8a-c). Among them, palmitic acid (C16:0) and stearic acid (C18:0) were the two most important SFA constitutions. Conversely, the MUFA content in SAT was significantly higher than that in VAT, but there was no significant difference with that in AAT. Oleic acid (C18:1n-9) made up the highest proportion of MUFA, followed by palmitoleic acid (C16:1). The PUFA contents in the three adipose tissues were stable, all within 3% of the total content. The desaturation index (DI) of SAT was significantly higher than that of VAT and AAT, while the elongation index (EI) was opposite, which might be related to the lipogenic gene activity of adipose tissue.

Fatty Acid Composition and Content
In this study, twenty-four fatty acid species were detected in all three adipose tissues (Table S12), including 11, 5, and 8 kinds of saturated fatty acids (SFA), monounsaturated fatty acids (MUFA), and polyunsaturated fatty acids (PUFA), respectively, which correspondingly accounted for 59.49-76.91%, 20.34-37.94%, and 2.52-2.79% of the total fatty acid content. The fatty acid content analysis represented that there was no significant difference in SFA content between VAT and AAT; however, both of them were significantly higher than those in SAT (Figure 8a-c). Among them, palmitic acid (C16:0) and stearic acid (C18:0) were the two most important SFA constitutions. Conversely, the MUFA content in SAT was significantly higher than that in VAT, but there was no significant difference with that in AAT. Oleic acid (C18:1n-9) made up the highest proportion of MUFA, followed by palmitoleic acid (C16:1). The PUFA contents in the three adipose tissues were stable, all within 3% of the total content. The desaturation index (DI) of SAT was significantly higher than that of VAT and AAT, while the elongation index (EI) was opposite, which might be related to the lipogenic gene activity of adipose tissue.

Lipid Composition and Content
Untargeted lipidomics technology was used to detect lipids from SAT, VAT, and AAT in the positive and negative ion modes. A total of 418 lipid species classified into 12 lipid classes were detected, including triglyceride (TG), glycerophosphocholine (PC), glycerophosphoethanolamine (PE), diacylglycerol (DG), ceramide (Cer), sphingomyelin (SM), glycerophosphoinositol (PI), lysoglycerophosphocholine (LPC), neutral glycosphingolipid, acyl carnitine (AcCA), sphingoid base, and sitosterol ester (SiE), among which TG was the main component of lipids making up 50.48% ( Figure S3). QC was used to standardize the lipid data. As shown in Figure 8d and Table S13, the TG content in SAT was significantly higher than those of VAT and AAT, while the relative contents of DG and LPC were significantly lower than those of VAT and AAT. There was no significant difference in the contents of 12 lipids classes between VAT and AAT, implying VAT and AAT were more similar.

Assessment of Lipid Data Quality
Plotting the time series diagram of the first principal component score of the QC sample. The small variation of PC1 of samples from QC1 to QC4 meant that the analysis quality of the platform was stable ( Figure S4a). The relative standard deviation (RSD) of more

Lipid Composition and Content
Untargeted lipidomics technology was used to detect lipids from SAT, VAT, and AAT in the positive and negative ion modes. A total of 418 lipid species classified into 12 lipid classes were detected, including triglyceride (TG), glycerophosphocholine (PC), glycerophosphoethanolamine (PE), diacylglycerol (DG), ceramide (Cer), sphingomyelin (SM), glycerophosphoinositol (PI), lysoglycerophosphocholine (LPC), neutral glycosphingolipid, acyl carnitine (AcCA), sphingoid base, and sitosterol ester (SiE), among which TG was the main component of lipids making up 50.48% ( Figure S3). QC was used to standardize the lipid data. As shown in Figure 8d and Table S13, the TG content in SAT was significantly higher than those of VAT and AAT, while the relative contents of DG and LPC were significantly lower than those of VAT and AAT. There was no significant difference in the contents of 12 lipids classes between VAT and AAT, implying VAT and AAT were more similar.

Assessment of Lipid Data Quality
Plotting the time series diagram of the first principal component score of the QC sample. The small variation of PC1 of samples from QC1 to QC4 meant that the analysis quality of the platform was stable ( Figure S4a). The relative standard deviation (RSD) of more than 94.0% of lipids was less than 30% ( Figure S4b), further proving the quality of lipid data was excellent. Generally, R2 and Q2 should be higher than 0.5, and the difference between them should be less than 0.3. The R2 and Q2 in SAT-VAT and SAT-AAT groups demonstrated that the model had good interpretation and prediction ability ( Figure S5a,b). However, for the VAT-AAT group (Figure S5c), the R2 and Q2 of the model represented that the difference between the VAT and AAT was smaller than that of the other two groups.

Screening of Differentially Expressed Lipids among Three Groups
Using the criteria with VIP > 1 and p < 0.05 (t-test) to screen DALs (Table S14). In SAT-VAT (Figure 9a), SAT-AAT (Figure 9b), and VAT-AAT (Figure 9c) groups, 97, 116, and 29 DALs, respectively, were identified. TG accounted for the highest proportion of the DALs. Additionally, two common DALs were identified among the three groups, which were TG (49:3) (monolinolenic triglyceride) and TG (54:5) (monostearic triglyceride). Therefore, TG played a decisive role in the development of fat deposition.

Pathway Enrichment of Differentially Expressed Lipids
MetaboAnalyst was used to perform pathway enrichment on the DALs identified in the three groups. A total of seven metabolic pathways were enriched (Figure 10), one of which was Glycosylphatidylinostiol (GPI)-anchor biosynthesis; three of which were enriched in pathways associated with lipid metabolism, including glycerophospholipid metabolism, glycerolipid metabolism, and sphingolipid metabolism; and three of which were enriched in pathways involved in unsaturated fatty acid metabolism, such as linoleic acid metabolism, α-Linoleic acid metabolism, and arachidonic acid metabolism. This study focused on pathways related to lipid metabolism and unsaturated fatty acid metabolism.

Pathway Enrichment of Differentially Expressed Lipids
MetaboAnalyst was used to perform pathway enrichment on the DALs identified in the three groups. A total of seven metabolic pathways were enriched (Figure 10), one of which was Glycosylphatidylinostiol (GPI)-anchor biosynthesis; three of which were enriched in pathways associated with lipid metabolism, including glycerophospholipid metabolism, glycerolipid metabolism, and sphingolipid metabolism; and three of which were enriched in pathways involved in unsaturated fatty acid metabolism, such as linoleic acid metabolism, α-Linoleic acid metabolism, and arachidonic acid metabolism. This study focused on pathways related to lipid metabolism and unsaturated fatty acid metabolism.

Pathway Enrichment of Differentially Expressed Lipids
MetaboAnalyst was used to perform pathway enrichment on the DALs identified in the three groups. A total of seven metabolic pathways were enriched (Figure 10), one of which was Glycosylphatidylinostiol (GPI)-anchor biosynthesis; three of which were enriched in pathways associated with lipid metabolism, including glycerophospholipid metabolism, glycerolipid metabolism, and sphingolipid metabolism; and three of which were enriched in pathways involved in unsaturated fatty acid metabolism, such as linoleic acid metabolism, α-Linoleic acid metabolism, and arachidonic acid metabolism. This study focused on pathways related to lipid metabolism and unsaturated fatty acid metabolism.

Transcriptomic and Lipid Metabolomic Integration Analysis
The Spearman correlation coefficients between 11 candidate genes and important metabolites including FAs and DALs were calculated for transcriptomic and lipid metabolomic integration analysis (Table S15). In contrast to PCK1, the other candidates had the same regulatory mode on FAs and DALs (Figure 11), namely, these genes had positive regulations with MUFA and TG, whereas negatively regulated SFA, PUFA, DG, and LPC.

Transcriptomic and Lipid Metabolomic Integration Analysis
The Spearman correlation coefficients between 11 candidate genes and important metabolites including FAs and DALs were calculated for transcriptomic and lipid metabolomic integration analysis (Table S15). In contrast to PCK1, the other candidates had the same regulatory mode on FAs and DALs (Figure 11), namely, these genes had positive regulations with MUFA and TG, whereas negatively regulated SFA, PUFA, DG, and LPC.
Genes 2023, 14, x FOR PEER REVIEW 14 of 21 Figure 11. Correlation analysis between candidate genes related to fat deposition and differential fatty acids and lipids. The legend on the right shows the correlation coefficient. Red and blue represent the positive and negative correlation, respectively. The asterisk indicates significance: ** p < 0.01, * p < 0.05.

Discussion
Fat deposition occurs with the increase and differentiation of adipocytes, which was mediated by multiple adipogenic genes and catalytic enzymes [38,39]. Prior studies showed that the functional differences in different adipose tissues were caused by intrinsic gene expression [40]. The study of gene functions related to bovine fat deposition has become a prevalent research hotspot. In this study, through the transcriptomics and metabolomics analysis of SAT, VAT, and AAT, 11 candidate genes were identified to affect fat deposition via their involvement in fatty acid synthesis, metabolism and oxidation, pyruvate metabolism, and carbohydrate synthesis and digestion (Figure 12), including fatty acid synthesis genes (ACACA, SCD, and ELOVL6), fatty acid oxidation genes (EHHADH and CPT2), ketoacid metabolism genes (ACADSB and OCXT1), glycolytic genes (PCK1 and PCK2), fatty acid transporter genes FABP7, and CHPT1. Additionally, it was further confirmed that PPAR signaling pathway was a pivotal pathway associated with fat deposition [41,42]. Figure 11. Correlation analysis between candidate genes related to fat deposition and differential fatty acids and lipids. The legend on the right shows the correlation coefficient. Red and blue represent the positive and negative correlation, respectively. The asterisk indicates significance: ** p < 0.01, * p < 0.05.

Discussion
Fat deposition occurs with the increase and differentiation of adipocytes, which was mediated by multiple adipogenic genes and catalytic enzymes [38,39]. Prior studies showed that the functional differences in different adipose tissues were caused by intrinsic gene expression [40]. The study of gene functions related to bovine fat deposition has become a prevalent research hotspot. In this study, through the transcriptomics and metabolomics analysis of SAT, VAT, and AAT, 11 candidate genes were identified to affect fat deposition via their involvement in fatty acid synthesis, metabolism and oxidation, pyruvate metabolism, and carbohydrate synthesis and digestion (Figure 12), including fatty acid synthesis genes (ACACA, SCD, and ELOVL6), fatty acid oxidation genes (EHHADH and CPT2), ketoacid metabolism genes (ACADSB and OCXT1), glycolytic genes (PCK1 and PCK2), fatty acid transporter genes FABP7, and CHPT1. Additionally, it was further confirmed that PPAR signaling pathway was a pivotal pathway associated with fat deposition [41,42]. sic gene expression [40]. The study of gene functions related to bovine fat deposition has become a prevalent research hotspot. In this study, through the transcriptomics and metabolomics analysis of SAT, VAT, and AAT, 11 candidate genes were identified to affect fat deposition via their involvement in fatty acid synthesis, metabolism and oxidation, pyruvate metabolism, and carbohydrate synthesis and digestion (Figure 12), including fatty acid synthesis genes (ACACA, SCD, and ELOVL6), fatty acid oxidation genes (EHHADH and CPT2), ketoacid metabolism genes (ACADSB and OCXT1), glycolytic genes (PCK1 and PCK2), fatty acid transporter genes FABP7, and CHPT1. Additionally, it was further confirmed that PPAR signaling pathway was a pivotal pathway associated with fat deposition [41,42]. Fatty acid composition is a complex quantitative trait that determined the hardness and oiliness of adipose tissue and then affected the sensory quality and nutritional value of meat [43]. In mammalian cells, FAs undergo elongation, desaturation, β-oxidation, and peroxidation to synthesize complex lipids such as TG, phospholipids, and sphingolipids Fatty acid composition is a complex quantitative trait that determined the hardness and oiliness of adipose tissue and then affected the sensory quality and nutritional value of meat [43]. In mammalian cells, FAs undergo elongation, desaturation, β-oxidation, and peroxidation to synthesize complex lipids such as TG, phospholipids, and sphingolipids for mediating the stabilization of intracellular lipids [44]. Acetyl-CoA carboxylase (ACC) is a rate-limiting enzyme essential for de novo fatty acid synthesis. Its cytoplasmic isoform ACCα is encoded by acetyl-CoA carboxylase α (ACACA) gene [45]. The C-terminal carboxytransferase domain of ACACA catalyzed the carboxylation of acetyl-CoA to malonyl-CoA and, therefore, was involved in the first step of fatty acid synthesis [46]. This study presented that ACACA was positively correlated with total MUFA and TG, which roughly coincided with the study carried out by Zhang et al. (2010), who found that SNP3 genotype of ACACA was significantly correlated with TG in hybrid beef cattle [47]. Moreover, ACACA was significantly related to fat-related traits in beef cattle, such as backfat thickness [47] and intramuscular fat deposition [48]. Interestingly, in our previous studies, ACACA had been viewed as the most promising candidate affecting subcutaneous fat deposition in Chinese Simmental beef cattle [32]. In this study, ACACA was significantly enriched in multiple GO terms and pathways related to fat deposition and, thus, it was forecasted to be a molecular marker of fat deposition in beef cattle.
Two isoforms of stearoyl-CoA desaturase (SCD), namely SCD1 and SCD5, have been identified in cattle, of which SCD1 mRNA was highly expressed in adipose tissue [49]. Lehnert et al. (2007) reported that increased SCD activity in SAT would lead to high MUFA content [50]. Correspondingly, the current study showed that SCD was significantly positively related to MUFA content. Additionally, the C18:0 content in SAT was significantly lower than that of VAT and AAT, but the expression level of SCD was higher, which confirmed the increase of C18:0 inhibited SCD expression [51]. Therefore, the high expression of SCD in SAT was mainly caused by the lower SFA and higher MUFA contents in SAT. Previous studies had shown that SCD could regulate bovine FA content and TG synthesis [52,53]. SCD activity in bovine adipose tissue was positively correlated with the fatty acid desaturation index [54], which emerged as a regulator of the transformation from SFA to MUFA. MUFA is the major substrate for TG synthesis, hence it could be speculated that SCD exerted effects on fat deposition via regulating MUFA content.
Elongation of long-chain fatty acids family member 6 (ELOVL6) was another fatty acid synthesis gene identified in this study, which encoded an enzyme that was a key regulator of the first reaction in the elongation cycle of C12-16 long-chain fatty acids (LCFAs) in adipose tissue, such as catalyzing the extension of C16:0 to C18:0 and palmitoleic acid (C16:1n-7) to vaccenic acid (C18:1n-7) [55]. Interfering with ELOVL6 increased the content of C16:0 but decreased the content of C18:1n9 [56]. Increased PUFA content significantly inhibited ELOVL6 expression [57]. Consistent with these findings, this study demonstrated that ELOVL6 had a positive correlation with the contents of a series of MUFA, such as C14:1, C16:1, C17:1, and C18:1n9, but negatively correlated with SFA and PUFA. Junjvlieke et al. (2020) reported that ELOVL6 could regulate adipogenesis and lipolysis of bovine adipocytes through FoxO, cAMP, and Wnt signaling pathways [58]. It was found that the bovine ELOVL6 gene activity was related to the expression of peroxisome proliferator activated receptor γ (PPARγ) in intramuscular fat [59] and subcutaneous fat [60]. In the study of Qinchuan cattle, ELOVL6 activity was positively correlated with PPARγ expression, which promoted adipocyte proliferation in SAT and regulated lipid metabolism via regulating the expression of cell cycle genes [61]. Additionally, studies had also found that ELOVL6 was not only an upstream inhibitor of lipid synthesis in perinatal dairy cows [62], but also a regulator of intramuscular fat deposition in yaks [63]. Consequently, ELOVL6 could be recognized as a potentially important effector of fat deposition.
Fatty acid oxidation has also been recognized as the major biological process to regulate fat deposition. Intracellular fatty acid β oxidation begins with the regulation of the translocation of activated fatty acid acyl-CoA from the cytoplasm to the mitochondrial matrix, and this process was mediated by carnitine palmitoyltransferases 1 and 2 (CPT1 and CPT2) [64]. The expression of CPT2 could be reduced when fatty acid oxidation stimulator PPARγ was silenced [65], indicating CPT2 could participate in fatty acid oxidation by interacting with PPARγ. Enoyl-CoA hydratase and 3-hydroxyacyl CoA dehydrogenase (EHHADH) is essential for the production of medium-chain dicarboxylic acids. As a regulator, it was involved in the second and third steps of fatty acid β-oxidation [66]. Overexpression of EHHADH compensated for the functional defect of HSD17B4; therefore, it could replace HSD17B4 to regulate the oxidation process of lauric acid (C12:0) [67]. EHHADH binding to PPARα provided a positive feedback loop to regulate the expression level of PPARα according to metabolic demands and then maintained lipid homeostasis [68]. 3-Oxoacid CoA-transferase 1 (OXCT1) is another rate-limiting enzyme indirectly involved in fatty acid β-oxidation through ketone body metabolism. In extrahepatic tissues, ketone bodies were decomposed by OXCT1 and eventually participated in the TCA cycle in the form of acetyl-CoA to produce ATP [69]. This process provided substrates for the biosynthesis of lipids and sterols in many tissues. Presently, research on OXCT1 mainly focused on human cancers and ketoacidosis [69]; however, its functional study in adipose tissue of domestic animals was limited. Differential expression of OXCT1 gene was identified in back adipose tissue of the F2 generation of Korean soil pigs and Yorkshire pigs [70]. Downregulated OXCT1 in ovine adipocytes promoted lipid accumulation, further emphasizing the importance of OXCT1 in adipogenesis [71]. More interestingly, OXCT1 expression was down-regulated in adipose tissue of obese subjects compared with healthy lean subjects [72]. Therefore, OXCT1 might be involved in the TCA cycle through ketone body metabolism to affect fat deposition.
Acyl-coenzyme A dehydrogenase, short/branched-chain (ACADSB), found in mitochondria, regulates the catabolism of branched-chain amino acids and FAs [73]. ACADSB was used as a target for the diagnosis and treatment of malignant diseases. Down-regulation of ACADSB inhibited fatty acid catabolism and then led to lipid accumulation [74]. In bovine mammary epithelial cells (bMECs), ACADSB participated in the metabolism of TG and FA [75]. The contents of TG and FA in bMECs were reduced after the knockdown of this gene, demonstrating that ACADSB was a key factor affecting milk fat synthesis and lipid metabolism [76].
In vertebrates, two phosphoenolpyruvate carboxykinase (PEPCK) isoenzymes are detected, in which the cytosolic isoenzyme PEPCK-C and mitochondrial isoenzyme PEPCK-M are encoded, respectively, by phosphoenolpyruvate carboxykinase 1 (PCK1) and phosphoenolpyruvate carboxykinase 2 (PCK2) [77]. PCK2 was recognized as a potential regulator of fat deposition in broilers [78]. PCK1 catalyzed oxaloacetate to phosphoenolpyruvate via GTP and then released GDP for energy metabolism [79]. This biological process could produce glycerol to regulate adipogenesis [80]. Additionally, PCK1 was a target gene of PPARγ and regulated phosphoenolpyruvate through the PPARγ signaling pathway to affect the production of FAs [81], indicating PCK1 could interact with PPARγ to participate in the processes concerning fat deposition [82]. Moreover, this gene had been screened to be associated with fat deposits in porcine [83] and chicken [37]. However, up to now, the regulatory mechanism of PCK1 gene in bovine fat deposition was still unclear, and its function should be further explored.
Fatty acid binding proteins (FABPs) affect the uptake, transport, and metabolism of fatty acids. FABP7 played an important role in the uptake of LCFAs in cell membranes [84]. Its high expression promoted cell proliferation by activating PPARγ signaling pathway [85]. FABP3 and FABP4, members of the same family as FABP7, had been reported to be significantly associated with marbling, subcutaneous fat thickness, and fatty acid composition in cattle [86][87][88]. However, there were no studies reported that FABP7 was related to fat deposition. In this study, FABP7 was an overlapped DEG identified by transcriptome differential analysis between three groups, and it was significantly enriched in PPAR signaling pathway, which highlighted the importance of FABP7 in fat deposition. Cholinephosphotransferase 1 (CHPT1) regulated choline metabolism by catalyzing phosphatidyl choline synthesis [89]. Inactivation of CHPT1 could restrain phospholipid synthesis [90]; thus, this gene played a central role in bubble membrane formation, while its role in fat deposition was inconclusive.

Conclusions
Collectively, the integrative analysis of transcriptomics and lipid metabolomics in SAT, VAT, and AAT revealed 11 candidates including fatty acid synthesis genes (ACACA, SCD, and ELOVL6), fatty acid oxidation genes (EHHADH and CPT2), ketoacid metabolism genes (ACADSB and OCXT1), glycolytic genes (PCK1 and PCK2), fatty acid transporter genes FABP7, and CHPT1 were all identified to be involved in fatty acid synthesis, metabolism and oxidation, pyruvate metabolism, and carbohydrate synthesis and digestion, which directly or indirectly participated in TG synthesis and thus affected fat deposition. Among them, FABP7 and PCK1 were differentially expressed in three adipose tissues. These findings could provide deeper insights into the regulation mechanisms of fat deposition in beef cattle.  Figure S5 OPLS-DA score plots in SAT-VAT (a), SAT-AAT (b), and VAT-AAT (c) comparative groups. R2X, R2Y, and Q2Y were the three parameters describing the quality of the OPLS-DA model, where the first two parameters represented respectively the explanatory ability of X variations and Y variations explained by all extracted components, and the third parameter denoted the predictive ability of the model. Table S1 Statistical analysis of transcriptome sequencing of different adipose tissues; Table S2 The differentially expressed genes identified in SAT-VAT group; Table S3 The differentially expressed genes identified in SAT-AAT group; Table S4 The differentially expressed genes identified in VAT-AAT group; Table S5 The overlaps of differentially expressed genes among SAT-VAT, SAT-AAT, and VAT-AAT groups; Table S6 GO annotation of differentially expressed genes identified in SAT-VAT, SAT-AAT, and VAT-AAT groups; Table S7 Pathway enrichment of differentially expressed genes identified in SAT-VAT, SAT-AAT, and VAT-AAT groups; Table S8 Weighted gene co-expression network analysis; Table S9 The enriched GO terms and KEGG pathways for the genes in significant modules; Table S10 The significant GO terms and KEGG pathways related to fat deposition and their involved genes in the lightcyan module; Table S11 The average FPKM expression levels of candidate genes in different adipose tissues; Table S12 The compositions and contents of fatty acid in different adipose tissues; Table S13 Comparison of lipid contents in different adipose tissues; Table S14 The differentially expressed lipis screened in SAT-VAT, SAT-AAT, and VAT-AAT groups; Table S15 The Spearman correlation coefficients between the FPKM expression of 11 candidate genes and the content of important metabolites.  Informed Consent Statement: Not applicable.

Data Availability Statement:
The raw sequencing reads of subcutaneous, visceral, and abdominal adipose tissues have been submitted to Sequence Read Archive (SRA) with accession number PR-JNA824614. Given that the paper has not yet been published, the data will be accessible on 1 May 2023, but can be obtained in advance from the corresponding author upon reasonable request.