Integrated lipidomics and RNA sequencing analysis reveal novel changes during 3T3-L1 cell adipogenesis

Cell culture and 3T3-L1 cell differentiation

3T3-L1 cells were obtained from Prof. Shulin Yang (Institute of Animal Sciences, Chinese Academy of Agricultural Sciences). Dulbecco’s Modified Eagle’s Medium (DMEM) (11965; Gibco) supplemented with 10% FBS (12483, Gibco) and 1% penicillin-streptomycin (15240; Gibco) (maintaining medium) was used to culture 3T3-L1 cells. Cultures were conducted with 1 ×10⁶ cells/well in 6-well cell culture plates. These cells were grown until confluence (day 0), then differentiation was induced by adding 1 mM insulin (I5523; Sigma Aldrich), 0.25 mM dexamethasone (D175; Sigma Aldrich), and 0.5 mM 3-isobutyl-1- methylxanthine (15879, Sigma Aldrich) in the maintaining medium. After 4 days, cells were cultured in maintaining medium with 1 mM insulin. On day 10, the fully differentiated adipocytes were used for Oil Red O (O1391; Sigma Aldrich) staining. Cells at day 0 (D0) and day 10 (D10) were collected and stored at −80 °C for subsequent experiments.

Nile Red O Staining

Cells were treated with 4% PFA for 10 min. Afterwards, cells were washed with PBS and incubated with Nile Red O (O1391; Sigma Aldrich) for 30 min. After staining, images were acquired using a LEICA DMi8 microscope.

Lipid extraction and lipidomics study

Collected cells were thawed on ice, and lipids were extracted with isopropanol (IPA). A total of 120 µL precooled IPA was added, then samples were vortexed for 1 min, and incubated for 10 min at room temperature. The extraction mixture was then stored at −20 °C overnight. After centrifugation at 4,000 g for 20 min, the supernatants (one per specimen) were diluted 1:10 with IPA/acetonitrile (ACN)/H₂O (2:1:1, v:v:v) and stored at −80 °C before analysis. A total of 10 µL of each extraction mixture was removed to prepare the pooled QC samples.

All lipid samples were analyzed by liquid chromatography–mass spectrometry (LC-MS) using an ultra-high performance liquid chromatography (UPLC) system (SCIEX, UK) and a high-resolution tandem mass spectrometer TripleTOF5600plus (SCIEX, UK). The acquired LC–MS raw data were analyzed by XCMS software (SCIEX, Warrington, UK), and the retention time (RT) and M/Z data were used to identify each ion. The online Kyoto Encyclopedia of Genes and Genomes (KEGG), Human Metabolome Database (HMDB), and in-house databases were used to perform level-one and level-two identification and annotation. MS2 indicates the metabolites that not only matched with level-one fragment ions, but also matched with level-two fragment ions in in-house database. Lipid metabolites with VIP (variable importance in projection) ≥1 and fold change > 2 or < 2 were considered statistically significantly different (Li et al., 2018).

mRNA library construction and sequencing

We used TRIzol reagent (Invitrogen, Waltham, MA) to exact total RNA from cells in accordance with the manufacturer’s procedure. Then we analyzed the RNA quantity and purity. Poly-T oligo-attached magnetic beads were used to purify Poly(A) RNA from total RNA with two rounds of purification. After purification, the mRNA was divided into small pieces at a high temperature using divalent cations. Then, according to the procedures of the mRNA Seq sample preparation kit, the cut RNA fragments were reverse transcribed to generate the final cDNA library (Illumina, San Diego, CA, USA). FastQC (version 0.11.2) was used to evaluate the quality of sequenced data. HISAT2 (version 2.0.4) was used to get clean data compared to the Mus musculus genome (Ensembl v101) (Kim, Langmead & Salzberg, 2015). StringTie (version 1.3.4) and Gffcompare (version 0.9.8) were used to assemble and quantify the transcripts, respectively (Pertea et al., 2015; Pertea et al., 2016). Differential gene expression between two groups (three biological replicates per condition) was performed using the DESeq2 R package available from Bioconductor (Love, Huber & Anders, 2014). Genes with an adjusted p-value < 0.05 and fold change > 2 or < −2 found by DESeq were assigned as differentially expressed. The statistical power, calculated using RNASeqPower (https://doi.org/doi:10.18129/B9.bioc.RNASeqPower) was 0.9829191. RNA sequencing data can be accessed on the SRA database, using accession number PRJNA795061.

Quantitative real time PCR (qRT-PCR)

qRT-qPCR was performed using TaKaRa SYBR Premix EX Taq (TaKaRa RR420A, JAPAN) on a QuantStudio5 Real-Time PCR System (Applied Biosystems, Waltham, MA, USA). All experiments contained three biological replicates, and each sample was quantified in duplicate. SPSS software and Excel were used to analyze the data. Relative gene expression was analyzed as 2^−ΔΔCt, and students t -test was used to calculate statistical difference. Data was shown as mean ± standard error (SE). Primer sequences are shown in Table S1.

Correlation analysis

Three sets of transcriptome data and three sets of lipidomics data were used for correlation analysis. Correlation analysis was performed using the OmicStudio tools at https://www.omicstudio.cn/tool.

3T3-L1 adipogenesis

Oil Red O staining showed that accumulation of lipid droplets within the cells increased during induction of differentiation (Figs. 1A–1F). In order to identify the changes between undifferentiated and differentiated 3T3-L1 cells, we collected samples at D0 and D10 for subsequent testing.

LC-MS-based untargeted lipidomics

We analyzed a total of 12 3T3-L1 preadipocyte and mature adipocyte samples (n = 6 per group) by LC-MS-based untargeted lipidomics to reveal differences in lipid composition. A total of 5,335 features were detected in positive ion mode, of which 3,187 could be annotated, with 850 matches to our in-house database. A total of 1,554 features were detected in negative anion mode, 873 of which could be annotated, with 378 matches to the in-house database. A total of 1,228 MS2 metabolites were identified from positive and negative ion modes (Table S2). Among all of the detected features, there were more up-regulated metabolites than down-regulated metabolites in mature fat cells (Fig. 2A).

Lipid profiles were compared using principal component analysis (PCA) (Fig. 2B). In unsupervised mode, the samples clustered together by cell type. Partial least squares discriminant analysis (PLS-DA) was used to identify the altered metabolites, and we found significant differences between the D0 and D10 groups (Figs. 2C, 2D). These results demonstrated significant lipid metabolite changes during differentiation of 3T3-L1 cells.

The q-values obtained by Benjamini–Hochberg (BH) correction using univariate fold-change analysis, the t-test, and the VIP value obtained by PLS-DA were analyzed with multivariate statistics and used to screen for differentially expressed lipid metabolites. Based on VIP values and relative abundance, we identified 454 differentially expressed MS2 metabolites between D10 and D0, among which 301 were up-regulated and 153 were down-regulated (Fig. 3A). These differentially expressed metabolites included 214 glycerophospholipids (GPs), 152 glycerolipids (GLs), 73 sphingolipids (SPs), ten fatty acyls (FAs), three sterol lipids (ST), and two prenol lipids (PLs) (Fig. 3B, Table S3). Among the 214 differentially expressed GPs, there were 50 PCs (26 decreased, 24 increased), 36 LPCs (16 decreased, 20 increased), 29 EtherPEs (17 decreased, 12 increased), 23 PEs (4 decreased, 19 increased), 17 PIs (one decreased, 16 increased), 11 EtherPCs (eight decreased, three increased), 11 LNAPEs (one decreased, 10 increased) and nine LPEs (two decreased, seven increased). Among the 152 differentially abundant GLs, there were 66 TGs (two decreased, 64 increased), 27 EtherTGs (four decreased, 23 increased), 25 DGs (all increased), and 23 EtherMGDGs (18 decreased, five increased). The 73 differentially expressed SPs included 28 SMs (16 decreased, 12 increased), 12 Cer_NS (three decreased, nine increased), nine HexCer_NS (4 decreased, five increased) and eight SHexCer (5 decreased, three increased).

The top 10 down-regulated lipids were PS 18:0_22:6, LPC 30:0-SN1, CAR 18:1, CAR 16:0, PE O-18:2_22:5, LPC 36:2-SN1, PC O-18:1_18:1, LDGCC 40:6, PE 38:4, and SMGDG O-11:0_28:6. The top 10 up-regulated lipids were all TGs: TG 16:1_16:1_17:1, TG 15:0_16:0_16:1, TG 14:1_16:1_17:1, TG 10:0_16:0_16:1, TG 16:0_16:1_16:2, TG 14:1_16:1_16:1, TG 14:0_16:1_16:1, TG 14:0_15:0_16:1, TG 8:0_16:0_16:1, and TG 14:0_14:1_16:1 (Table 1 and Table S3).

The differentially expressed lipids were analyzed for biochemical pathway enrichment using the KEGG database (Fig. 4). We found that differentially regulated lipid molecules in the mature adipocytes were broadly related to metabolism (e.g., ether lipid metabolism and SP metabolism) and organismal systems (e.g., adipocytokine signaling pathway and regulation of lipolysis in adipocytes).

Differential gene expression between undifferentiated and differentiated 3T3-L1 cells

Cell samples at D0 and D10 were analyzed via RNA-sequencing (RNA-seq). The main characteristics of the libraries are shown in Table S4. The libraries contained 49,519,629 raw reads on average. After removing adaptors and low-quality/ambiguous sequences, an average of 47,209,124 valid clean reads remained. Among all samples, 96.45% of the valid reads mapped to the mouse genome database, including 72.76% unique mapped reads and 23.69% multi-mapped reads (Table S5). The gene expression profiles of D0 and D10 samples were analyzed with PCA, which revealed significant differences in gene expression patterns between D10 and D0 (Fig. 5A). Volcano plots were used to visualize the distribution of differentially expressed genes (DEGs) between D10 and D0 cells. There were 6,193 DEGs between D10 and D0, including 1,878 up-regulated and 4,315 down-regulated genes (Fig. 5B). The heatmap comparison represents some gene expression changes (Fig. 5C).

We next validated the RNA-seq data with qRT-PCR for randomly selected genes. The relative expression levels of acyl-coenzyme A dehydrogenase medium chain (Acadm), acyl-CoA synthetase long-chain family member 1 (Acsl1), angiopoietin-like 4 (Angpt14), fatty acid binding protein 5 (Fabp5), lipase hormone sensitive (Lipe), patatin-like phospholipase domain containing 2 (Pnpla2), stearoyl-Coenzyme A desaturase (Scd1), and lipoprotein lipase (Lpl) were significantly increased in differentiated 3T3-L1 cells (D10, p < 0.05, Fig. 6). The expression trends for all genes validated by qRT-PCR were consistent with the results from RNA-seq analysis, demonstrating the high quality of the sequencing data.

The 6193 DEGs identified with RNA-seq were analyzed for functional enrichment using the Gene Ontology (GO) database. There were significantly enriched biological processes (e.g., lipid metabolic process, fatty acid biosynthetic and lipid transport), cellular components, and molecular functions (Fig. S1). Enrichment analysis was also conducted using the KEGG database (Fig. 7). The integrated DEGs of undifferentiated and differentiated 3T3-L1 cells were mainly enriched in thermogenesis, the PPAR signaling pathway, and regulation of lipolysis in adipocytes.

Joint analysis of RNA-seq and lipid metabolome results

We compared the enriched KEGG pathways between DEGs and differentially expressed lipid metabolites. There were seven pathways enriched in both datasets (p < 0.05), namely adipocytokine signaling pathway, AGE-RAGE signaling pathway in diabetic complications, glycerolipid metabolism, insulin resistance, oxidative phosphorylation, regulation of lipolysis in adipocytes, and retrograde endocannabinoid signaling (Fig. 8 and Table S6). Genes with a fold change >5 and FPKM >1 (Table 2) were strongly correlated with metabolites in these pathways (Fig. S2 and Table S7). All of the DEGs were screened and a protein interaction network was constructed using String. The key hub genes were identified with Cytohubba in Cytoscape. The top ten hub genes were Pnpla2, diacylglycerol O-acyltransferase 1 (Dgat1), diacylglycerol O-acyltransferase 2 (Dgat2), fatty acid binding protein 4 (Fabp4), adiponectin C1Q and collagen domain containing (Adipoq), Lipe, Acsl1, solute carrier family 2, member 4 (Slc2a4), Lpin1, and Lpl (Fig. 9).

Expression patterns of hub genes and key genes that regulate lipid metabolism during adipogenesis

Expression patterns of the top ten hub genes were analyzed in cell samples collected at D0, D2, D4, D6, D8, and D10 during adipogenesis. In addition, expression of the key adipogenesis genes CCAAT enhancer binding protein alpha (Cebp/ α) and peroxisome proliferator activated receptor gamma (Ppar γ) were also quantified; both were significantly activated during differentiation (Fig. S3). Pnpla2 (Fig. 10A), Dgat1 (Fig. 10B), Fabp4 (Fig. 10D), and Lipe (Fig. 10F) were significantly upregulated during differentiation, with the highest expression at D4. Dgat2 was highly expressed at D2 and D4 during differentiation, but the expression level was lower than in undifferentiated cells at D10 (Fig. 10C). Adipoq (Fig. 10E), Acsl1 (Fig. 10G), Lpin1 (Fig. 10H), Lpl (Fig. 10I), and Plin1 (Fig. 10J) were significantly upregulated during differentiation, with the highest expression at D8.