Introduction

Physical activity is one of the important ways to maintain health and can regulate the function of major organs such as heart, brain and liver1. The benefits of scientific exercise are the result of the synergy, interaction, and communication of multiple physiological systems. The positive effects of physical exercise on spatial learning and memory had been reported both in animal and human studies2. Increasing physical activity through exercise can improve fatty liver disease3. Gao et al. reported that exercise promoted lipophagy by activating AMPK/ULK1 pathway to improve nonalcoholic fatty liver disease (NAFLD)4. Hu et al. indicated that exercise alleviated the mitochondrial dysfunction to improve NAFLD by activating Sirt1-mediated Drp1 acetylation5. Regular physical exercise during pregnancy can help preventing pregnancy related disorders, such as gestational diabetes, hypertensive disorders, anxiety and prenatal depression6. Physical exercise had potential roles in prevention, treatment, and prognoses improvement of cancer patients7. The molecular mechanisms underlying exercise are still worth exploring.

Lactate, initially considered as metabolic waste products in glycolysis, not only serves as an energy source for metabolism and a major gluconeogenic precursor, but also has been shown to regulate intracellular signal transduction as a signaling molecule8. Lactate is involved in a variety of physiological and pathological processes, such as energy regulation, regulation of fatty acid metabolism, ischemic injury, inflammatory responses, and tumor growth and metastasis9. Lactylation, as a novel histone modification that can epigenetically regulate gene transcription, was formally proposed for the first time in 201910. Intracellular hypoxia induces mitochondrial protein lactate to limit oxidative phosphorylation to reduce reactive oxygen species damage11. Elevated lactylation level in colorectal cancer tissues was associated with poor prognosis12.

The integration of multi-omics data can enhance the understanding the molecular behavior underlying exercise. Multi-omics data encompasses a diverse array of biological measurements, including genomic, transcriptomics, proteomic, epigenomic, and metabolomics. For instance, transcriptomics utilizes RNA sequencing to identify differential splicing and quantification of the relative levels of gene expression13. Proteomics, on the other hand, employs mass spectrometry to estimate the abundance and alterations of proteins, while metabolomics offers a snapshot of the metabolic state of cells14. Whole genome methylation was decreased in skeletal muscle biopsies obtained from healthy sedentary individuals following one acute bout of exhaustive aerobic exercise15. Metabolomic analysis revealed increased levels of polyunsaturated free fatty acids and lysophospholipids, accompanied by reductions in ceramides, sphingolipids and phospholipids, and alterations in gut microbiome metabolites after enhanced physical activity16. An integrated analysis of metabolomic and transcriptomic data in zebrafish subjected to an exercise intervention revealed alterations in metabolic pathways, metabolites, and expression of core genes, which may contribute to exercise-induced improvements in bone mass17. To date, research has predominantly focused on one or two omics domains. Therefore, multi-omics map of the effects of exercise is needed to understand the molecular basis of exercise training-induced adaptations.

In this present study, transcriptomics, proteomic, lactylation modification and metabolomics analyses were performed on liver tissues in mice after six weeks of moderate-intensity treadmill exercise. Gene sets related to hypoxia, glycolysis, and fatty acid metabolism were used to aid in the screening of hub genes. This multi-omics integration strategy aimed to reveal the potential impact of chronic exercise on the mice liver tissue metabolism, helping to increase the understanding of the molecular mechanisms underlying the physiological changes behind exercise.

Materials and methods

Experimental mice and chronic exercise protocol

Adult C57BL/6J male mice (20–22 g, 8 weeks), purchased from Charles River Laboratory, were used in this study. After 1-week adaptation, mice were randomly divided into two groups: the chronic exercise (Experimental group) and the sedentary control (Control group). All animal procedures were conducted during rodents’ nocturnal phase (6:00 pm to 12:00 pm). The animal experiments were complied with the institutional ethical guidelines, ARRIVE guidelines, and approved by the Ethics Committee on Animal Care and Use of the Capital Institute of Pediatrics (DWLL2021015). The moderate-intensity treadmill exercise was performed as described previously18,19. Mice from each group were anesthetized with isoflurane and sacrificed by cervical dislocation within 24 h after the last exercise. The liver tissues were collected following the same procedures as in the previous study18. Principal component analysis (PCA) was performed with “princomp” function in the R package “stats” (version 4.3.1). Transcriptomics, proteomics, and lactylation modification analysis was performed on liver tissue samples from 6 mice (n = 3, both groups). For metabolomics analysis, liver tissue samples from 12 mice (n = 6, both groups) were used.

Transcriptomics analysis

Total RNA extraction and RNA sequencing procedures were same as in the previous study18. RNA-seq count data were converted to transcripts per million (TPM), and TPM data were log2 (data + 1) processed to normalize the data. Differentially expressed genes (DEGs) were obtained with p < 0.05 & (fold change (FC) > 1.2 or FC < 0.83) using the limma (v 3.46.0), visualized by volcano maps and bar chart with ggplot2 package (version 3.5.1). In total, 200 hypoxia-related genes (HRGs) in HALLMARK_HYPOXIA, 200 glycolysis-related genes (GRGs) in HALLMARK_GLYCOLYSIS, and 158 fatty acid metabolism-related genes (FRGs) in HALLMARK_FATTY_ACID_METABOLISM were obtained from MSigDB (https://www.gsea-msigdb.org/gsea/msigdb/). We identified the intersection of DEGs with 200 HRGs, 200 GRGs, and 158 FRGs, respectively. Subsequently, we took the union of these intersected genes to obtain the intersection genes for further downstream analysis. David (version 2024q1, https://david.ncifcrf.gov/)20 was used to perform Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis for the intersection genes with p < 0.05. The protein-protein interaction (PPI) network of intersection genes was constructed with STRING (version 12.0, https://cn.string-db.org/)21.

Proteomics and lactylation modification analyses

Liquid chromatography-tandem mass spectrometry was performed at PTM Bio lab (Hangzhou, China) as previously described18. The abundance of lactylated protein have been normalized by the protein levels. Differentially expressed proteins (DEPs), differentially lactylation modification proteins (DLPs), and differentially lactylated sites (DLSs) were determined with p < 0.05 & (FC > 1.2 or FC < 0.83) based on a Student’s t-test. The results were visualized by volcano plots (ggplot2, version 3.5.1) and heat maps (pheatmap, version 1.0.12).

Identification of hub genes and transcription factors (TFs) regulatory network

The hub genes were obtained by overlapping the intersection genes with DEPs and DLPs. The gene expression, protein, and lactylation modification levels of hub genes were displayed in box plots with ggpubr package (version 0.6.0). The TFs regulating the hub genes were predicted using the Cistrome DB (http://genemania.org) with RP score > 0.7. Subsequently, the TFs regulatory network was visualized using Cytoscape (version 3.8.2)22.

Immune microenvironment infiltration analysis

The murine Microenvironment Cell Populations counter (mMCP-counter) was employed to evaluate infiltration of immune cells in murine liver tissues. It is a transcriptome-based tool for accurately quantifying tissue-infiltrated immune and stromal cell populations in mouse samples23. mMCP-counter can be accessed as an R package (https://github.com/cit-bioinfo/mMCP-counter). It takes a gene expression profiles matrix as input and returns the abundance of RNA originating from 16 defined cell populations present in the heterogeneous sample.

Correlation analysis

In addition, the correlation analysis between gene expression, protein, and lactylation modification levels of hub genes and immune cells was performed. Pearson correlation analysis was performed using psych package (version 2.4.6) in R to calculate the correlation coefficient and its corresponding p-value. The visual presentation of their correlations was shown by a heat map using the ggplot2 package (version 3.5.1).

Metabolomics analysis

Similarly, the metabolomics analysis was performed as previously described18. The normalization of metabolites was performed using MetaboAnalyst 5.0 (https://www.metaboanalyst.ca/MetaboAnalyst/docs/RTutorial.xhtml) to make metabolites concentration distribution is close to normal distribution. Differential metabolites (DEMs) were determined by variable importance value (VIP) > 1 with p < 0.05 based on a Student’s t-test. KEGG enrichment analysis was performed using Metabolite Biological Role (MBROLE) 2.0 (http://csbg.cnb.csic.es/mbrole2/index.php) with p < 0.05. In addition, the correlation analysis between hub genes and DEMs was performed. The psych package (version 2.4.6) in R was used to calculate the Pearson correlation coefficient and its corresponding p-value.

Statistical analysis

The PCA using the “princomp” function in the R package “stats” (version 4.3.1) was performed to evaluate the accuracy of grouping. The “limma” (3.46.0) package was used to obtain the DEGs. David (https://david.ncifcrf.gov/) was used to perform functional enrichment analysis. DEPs, DLPs, DLSs, and DEMs were determined with p < 0.05 & (FC > 1.2 or FC < 0.83) based on a Student’s t-test. mMCP-counter was used to perform immune infiltration analysis. Pearson correlation analysis was performed using psych package (version 2.4.6) in R to calculate the correlation coefficient and its corresponding p-value. Unless otherwise stated, a p-value < 0.05 was deemed statistically significant. The software packages used for each omics analysis with details were summarized in Table 1.

Table 1 The software packages used for each omics analysis in this study.
Full size table

Results

Identification of intersection genes

The PCA plot showed that there was a significant separation between the experimental and control groups (Figure S1A). Compared to the control group, 2034 DEGs, including 1456 up-regulated mRNAs and 578 down-regulated mRNAs, were detected in the experimental group (Fig. 1A, Table S1). Among them, 40 genes overlapped with HRGs, 23 genes overlapped with GRGs, and 36 genes overlapped with FRGs, resulting in a total of 82 intersection genes (Fig. 1B). GO analysis indicated that these intersection genes were enriched in lipid metabolic process, fatty acid metabolic process, and negative regulation of apoptotic process (Fig. 2A, Table S2). KEGG enrichment analysis revealed several significantly enriched pathways, such as metabolic pathways, fatty acid metabolism, PPAR signaling pathway, and HIF-1 signaling pathway (Fig. 2B, Table S2). PPI network was displayed in Fig. 3, including 72 nodes and 274 edges. The number of proteins interacting with Fasn was the highest (up to 21 proteins).

Fig. 1
figure 1

Identification of intersection genes. (A) Volcano and Venn plots of the DEGs. (B) Venn diagram displaying the DEGs and HRGs/GRGs/FRGs. DEG, differentially expressed genes; HRG, hypoxia-related gene; GRG, glycolysis-related gene; FRG, fatty acid metabolism-related gene.

Full size image
Fig. 2
figure 2

Functional enrichment of intersection genes. (A) GO analysis for intersection genes. (B) KEGG enrichment analysis for intersection genes. GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes.

Full size image
Fig. 3
figure 3

PPI network of intersection genes. Node size and color depth are proportional to degree.

Full size image

Proteomics and lactylation modification analyses

The PCA plot showed that there was a significant separation between the experimental and control groups (Figure S1B). Compared to the control group, 577 DEPs, including 253 up-regulated proteins and 324 down-regulated proteins, were detected in the experimental group (Fig. 4A-B, Table S3). Lactylation modification analysis identified 264 DLSs of 141 DLPs, including 7 up-regulated lactylated sites and 257 down-regulated lactylated sites (Fig. 4C-D, Table S4).

Fig. 4
figure 4

Proteomics and lactylation modification analyses. (A) Volcano plot of the DEPs. (B) Heatmap of the DEPs. (C) Volcano plot of the DLSs. (D) Heatmap of the DLSs. DEP, differentially expressed protein; DLS, differentially lactylated site.

Full size image

Identification of hub genes

Venn diagram shows the overlap among 577 DEPs, 141 DLPs, and 82 intersection genes (Figure S2A). There were 24 shared genes between DEPs and DLPs, 19 shared genes between DEPs and intersection genes, 10 shared genes between DLPs and intersection genes. Finally, three hub genes (Aldoa, Acsl1, and Hadhb) were obtained by overlapping 82 intersection genes with 577 DEPs and 141 DLPs. Among them, aldolase A, fructose-bisphosphate (Aldoa) is a gene related to glycolysis and hypoxia, and also a gene related to fatty acid metabolism. The gene expression levels of three hub genes were up-regulated in the experimental group (Fig. 5A). The protein expression levels of acyl-CoA synthetase long-chain family member 1 (Acsl1) and hydroxyacyl-CoA dehydrogenase trifunctional multienzyme complex subunit beta (Hadhb) were up-regulated in the experimental group, while Aldoa was down-regulated (Fig. 5B). Lactylation modification analysis revealed significantly decreased lactylation levels of Aldoa_K42 and Hadhb_K273 in the experimental group, as well as decreased lactylation level of Acsl1_K676 without significance (Fig. 5C). The TF regulatory network, including 3 hub genes and 55 TFs, was constructed (Figure S2B). Aldoa was regulated by 45 TFs, Acsl1 was regulated by 14 TFs, and Hadhb was regulated by 6 TFs. In particular, peroxisome proliferator activated receptor gamma (Pparg) regulated all three hub genes. Among the 55 TFs, the expression levels of 12 TFs (Brd2, Btaf1, Chd4, Creb1, Esr1, Ets1, Fos, Gata4, Max, Myc, Rxrg, and Stat1) were significantly different between the experimental and control groups (Figure S2C).

Fig. 5
figure 5

Characterization of hub genes. (A) Gene expression levels of hub genes. (B) Protein levels of hub genes. (C) Lactylation modification levels of hub genes. ns indicated not significant, *indicated p < 0.05, **indicated p < 0.01, ***indicated p < 0.001, ****indicated p < 0.0001.

Full size image

Immune microenvironment infiltration analysis

The immune infiltration analysis revealed a decreased score for monocytes/macrophages and an increased score for endothelial cells in the experimental group (Fig. 6A). The correlation analysis indicated that the gene expressions of three hub genes were negatively correlated with monocytes/macrophages and positively correlated with endothelial cells (Fig. 6B). The protein expressions of Acsl1 and Hadhb were negatively correlated with monocytes/macrophages and positively correlated with endothelial cells, while Aldoa was positively correlated with monocytes/macrophages and negatively correlated with endothelial cells (Fig. 6C). The lactylation levels of Aldoa_K42 and Hadhb_K273 were positively correlated with monocytes/macrophages and negatively correlated with endothelial cells (Fig. 6D).

Fig. 6
figure 6

Analysis of immune microenvironment. (A) The immune infiltration levels of two immune cells. (B) The correlation analysis between gene expression levels of hub genes and two immune cells. (C) The correlation analysis between protein expression levels of hub genes and two immune cells. (D) The correlation analysis between lactylation modification levels of hub genes and two immune cells. * indicated p < 0.05, ** indicated p < 0.01.

Full size image

Metabolomics analysis

The PCA score plot clearly showed a complete separation between the experimental and control groups in positive (Figure S1C) and negative (Figure S1D) mode. In positive mode, 459 DEMs were identified, among which 255 metabolites were significantly increased and 204 metabolites were significantly decreased, in the experimental group (Fig. 7A-B, Table S5). Of these 459 DEMs, KEGG enrichment analysis was performed for 81 metabolites with KEGG_ID using MBROLE 2.0. These 81 DEMs were significantly enriched in metabolic pathways, purine metabolism, alpha-linolenic acid metabolism, caffeine metabolism, and sphingolipid metabolism (Fig. 7C, Table S6). The “metabolic pathways” was a common pathway between intersection genes-enriched pathways and DEMs-enriched pathways in positive mode. The correlation analysis was performed between 26 metabolites in “metabolic pathways” and the gene and protein expression of hub genes (Fig. 7D-E). The gene expressions of Acsl1 was the most negatively correlated with PTM_3508 [(-)-Jasmonic acid] (r = -0.994) and the gene expressions of Aldoa was the most positively correlated with PTM_43 (Uridine 5’-diphosphate) (r = 0.973) (Fig. 7D). The protein expression of Acsl1 was the most negatively correlated with PTM_3508 [(-)-Jasmonic acid] (r = -0.996) and the protein expressions of Aldoa was the most positively correlated with PTM_3508 [(-)-Jasmonic acid] (r = 0.998) (Fig. 7E).

Fig. 7
figure 7

Metabolomics analysis in positive mode. (A) Volcano plot of the DEMs in positive mode (B) Heatmap of the DEMs in positive mode. (C) KEGG pathway annotation for DEMs in positive mode. (D) The correlation analysis between gene expression levels of hub genes and DEMs in Metabolic pathways. (E) The correlation analysis between protein expression levels of hub genes and DEMs in Metabolic pathways. DEM, Differential metabolite.

Full size image

In negative mode, 181 DEMs were identified, among which 84 metabolites were significantly increased and 97 metabolites were significantly decreased, in the experimental group (Fig. 8A-B, Table S7). Of these 181 DEMs, KEGG enrichment analysis was performed for 65 metabolites with KEGG_ID using MBROLE 2.0. These 65 DEMs were significantly enriched in metabolic pathways, purine metabolism, arachidonic acid metabolism, primary bile acid biosynthesis and pantothenate and CoA biosynthesis (Fig. 8C, Table S8). The “metabolic pathways” was a common pathway between intersection genes-enriched pathways and DEMs-enriched pathways in negative mode. The correlation analysis was performed between 33 metabolites in “metabolic pathways” and the gene and protein expression of hub genes (Fig. 8D-E). The gene expressions of Acsl1 was the most negatively correlated with PTM_104 (s-adenosylhomocysteine (r = -0.999) and the gene expressions of Aldoa was the most positively correlated with PTM_631 (4-pyridoxic acid) (r = 0.984) (Fig. 8D). The protein expression of Acsl1 was the most negatively correlated with PTM_104 (s-adenosylhomocysteine (r = -0.994) and the protein expressions of Aldoa was the most positively correlated with PTM_104 (s-adenosylhomocysteine (r = -0.995) (Fig. 8E).

Fig. 8
figure 8

Metabolomics analysis in negative mode. (A) Volcano plot of the DEMs in negative mode. (B) Heatmap of the DEMs in negative mode. (C) KEGG pathway annotation for DEMs in negative mode. (D) The correlation analysis between gene expression levels of hub genes and DEMs in Metabolic pathways. (E) The correlation analysis between protein expression levels of hub genes and DEMs in Metabolic pathways. DEM, Differential metabolite.

Full size image

Discussion

The importance of physical exercise in the treatment of various diseases and improvement of organ functions has been extensively reported24,25,26. The bioinformatics analysis of omics data helps to deepen the understanding of the molecular mechanism underlying exercise. A proteomics study indicated that exercise training could alter the hepatic protein expression profiles of aged rat, thereby regulating the key liver metabolic pathways27. Analysis of transcriptomic and metabolomics data revealed that early active exercise had significant effects on bone development and metabolism in mice28. In the present study, we identified three hub genes (Hadhb, Aldoa, and Acsl1) whose mRNA and protein expression levels and lactylation levels were significantly altered during chronic exercise based on mutil-omics data (transcriptome, proteome and lactatomics data).

Mitochondrial trifunctional protein (MTP), consisting of two α subunits (HADHA) and two β subunits (HADHB), locates on the inner mitochondrial membrane and plays important roles in fatty acid β-oxidation to participate in and regulate fatty acid metabolism29. ACSL1 (acyl-CoA synthetase long chain family member 1) contributes the majority of total acyl-CoA synthetase (ACS) activity in human skeletal muscle and also plays important roles in fatty acid metabolism30. Acsl1 deficiency in skeletal muscle of mice impairs fatty acid oxidation31,32. Elevated Acsl1 induced by exercise training has potential to increase the utilization of fatty acids through β-oxidation30. It has been reported that one genome-wide significant polymorphism (rs116143768) within the ACSL1 gene was associated with fat loss efficiency in response to aerobic exercise33. Stierwalt et al. suggested that the expression of Acsl1 may undergo regulation by exercise training34. In this study, both the gene and protein expression level of Acsl1 and Hadhb were increased in the liver tissues of experimental mice, suggesting that chronic exercise may elevate the expressions of enzymes involved in fatty acid β-oxidation, thereby enhancing the utilization of fatty acids as an energy source. This underscores the importance of exercise in promoting healthy fatty acid metabolism and potentially enhancing fat loss efficiency.

In addition, up to 21 proteins, including Acsl1 and Hadhb, interacted with Fasn in the PPI network. In this study, both the gene and protein expression levels of Fasn were increased in the liver tissues of experimental mice. Fatty acid synthase (FASN) is a key lipogenic enzyme catalysing the terminal steps in the de novo biogenesis of fatty acids35. In the short-term high-fat diet experiment of female C57BL/6 mice, the expression of Fasn was increased in the setting of exercise36. Hippocampal-specific up-regulation of FASN was induced by voluntary running in 20 week-old male mice, which plays an important role in exercise-mediated cognitive enhancement37. Taken together, these findings may provide new perspectives for understanding the health benefits of exercise.

ALDOA (aldolase A, fructose-bisphosphate) is a highly conserved enzyme that plays a role in the fourth step of glycolysis and participates in a variety of cellular functions and biological processes38,39. The genome resequencing on buffaloes revealed that ALDOA was involved in pathways associated with the response to exercise39. Compared to sedentary men, ALDOA was expressed at lower levels in the left vastus lateralis muscles of sprinters40. In this study, the gene expression level of Aldoa were up-regulated and protein expression level of Aldoa were down-regulated in the liver tissues of experimental mice. Taken together, these findings suggested the multifaceted role of Aldoa in cellular metabolism, particularly its involvement in glycolysis and the response to exercise. Further research into the specific mechanisms underlying the regulation and function of Aldoa may yield significant insights into human health and disease.

In the TF regulatory network, Pparg regulated all three hub genes (Hadhb, Acsl1, and Aldoa). Peroxisome proliferator activated receptor gamma (PPARG/PPARγ) is a member of the nuclear hormone receptor superfamily. Tong et al. reported that PPARγ could be the direct target of IRF6 in high-fat diet (HFD)-induced fatty liver41. Swimming exercise partly suppresses lipid accumulation in HFD-induced non-alcoholic fatty liver disease in mice by down-regulating PPARγ-mediated pathways42. Based on function enrichment analysis for the intersection genes, three hub genes were enriched in Metabolic pathways (mmu01100), which was a common pathway between intersection genes-enriched pathways and DEMs-enriched pathways. The genes, Hadhb and Acsl1, were enriched in Fatty acid degradation (mmu00071) and Fatty acid metabolism (mmu01212). These findings further illustrated that chronic exercise has a certain effect on mice liver tissue metabolism. In addition, the lactylation levels of hub genes significantly different between the control and experimental groups, which need to be further explored.

The immune cell infiltration analysis revealed a decreased infiltration of monocytes/macrophages and an increased infiltration of endothelial cells in the experimental group. The high level of CD68+ monocytes/macrophages was detected in muscles of muscular dystrophy, X-linked mouse model (mdx) mice, and low-intensity training of mdx mice further increased the CD68 level43. Intense aerobic exercise training increased glycolysis and reduced the formation of reactive oxygen species to alter the muscle microvascular endothelial cells properties44. The specific role of immune cells in the response to exercise needs to be further explored to add new insights into the effects of exercise on the immune microenvironment.

Some limitations exist in this study. First, the sample size used for sequencing of each omics was small. Second, there was a lack of multi-omics level validation and in-depth functional exploration of the three hub genes. Third, the association study among different omics was relatively superficial and not deep enough.

In summary, the findings of this study highlight the pivotal roles of three hub genes in the glycolysis and fatty acid metabolism, particularly in the context of chronic exercise. The moderate-intensity exercise may modulate the mice liver tissue metabolism by regulating the expression of Hadhb, Aldoa, and Acsl1. Future research exploring the molecular mechanisms underlying these findings may provide further insights into the role of exercise in promoting metabolic health. In the future, we will consider more in-depth studies on the correlation among different omics.