Exercise-responsive skeletal muscle genes mechanistically linked to metabolic dysfunction-associated steatotic liver disease

Abstract

BACKGROUND

Metabolic dysfunction-associated steatotic liver disease (MASLD) is a common chronic liver disease that progresses from simple steatosis to inflammation, fibrosis, and cirrhosis. Currently, no effective targeted therapy is available. Exercise is a well-recognized non-pharmacological intervention with clear benefits. However, the biological mechanisms by which skeletal muscle responds to regular exercise and contributes to MASLD improvement remain poorly understood.

AIM

To identify exercise-responsive biomarkers in skeletal muscle associated with MASLD and explore their diagnostic and therapeutic potential.

METHODS

We analyzed skeletal muscle transcriptomic datasets from the gene expression omnibus. Differentially expressed genes (DEGs) were detected and then analyzed using Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment methods. To identify key genes, we employed weighted gene co-expression network analysis (WGCNA) and least absolute shrinkage and selection operator (LASSO) regression. Correlation with diagnostic efficacy was performed utilizing a validation group and receiver operating characteristic (ROC) analysis. Finally, an obese mouse model was established and subjected to endurance aerobic training. Gastrocnemius muscle tissue was validated at the messenger RNA, protein, and secretion levels to confirm the identified biomarkers.

RESULTS

Transcriptomic analysis identified 61 DEGs between pre-exercise and post-exercise samples, with 40 upregulated and 21 downregulated genes. GO enrichment analysis showed that extracellular matrix (ECM) organization and collagen fibril formation were significantly enriched. KEGG pathway analysis further highlighted cytoskeleton dynamics in muscle cells and ECM-receptor interactions. Integrated DEGs, WGCNA and LASSO analysis identified 12 hub genes. Validation cohort and ROC analysis demonstrated strong diagnostic performance for nine hub genes (COL3A1, COL1A2, BGN, LAMB1, PECAM1, LAMA4, THBS4, PXDN and THY1). In the mouse model, three hub genes (Lama4, Pecam1 and Pxdn) were significantly upregulated, while Thbs4 was downregulated after exercise in skeletal muscle tissue.

CONCLUSION

This study identified four exercise-responsive skeletal muscle-expressed genes (LAMA4, PECAM1, PXDN and THBS4). These genes are mechanistically associated with MASLD and may serve as myokine-like candidates. Our research offers new perspectives on the pathophysiology of MASLD and suggest possible strategies for precision diagnosis and therapy.

Key Words: Metabolic dysfunction-associated steatotic liver disease; Myokines; Weighted gene co-expression network analysis; LAMA4; PECAM1; PXDN; THBS4

Core Tip: This study identified four exercise-responsive skeletal muscle-expressed genes (LAMA4, PECAM1, PXDN and THBS4), which were mechanistically linked to metabolic dysfunction-associated steatotic liver disease (MASLD). These genes may serve as myokine-like candidates, shedding light on the molecular pathways by which exercise mediates metabolic adaptation in MASLD. The findings provided novel insights into the pathophysiology of MASLD and suggested these biomarkers as promising candidates for precision diagnostics and therapeutic interventions.

INTRODUCTION

Metabolic dysfunction-associated steatotic liver disease (MASLD), formerly known as non-alcoholic fatty liver disease, is characterized by hepatic steatosis in the absence of hepatocellular injury[1]. MASLD has emerged as the leading chronic liver condition globally, impacting > 38% of adults and 7%-14% of children, with projections indicating a 55% increase by 2040[2]. MASLD is closely linked to obesity, insulin resistance, and type 2 diabetes mellitus (T2DM), with up to 70% of T2DM patients developing MASLD[3,4]. Alarmingly, the rising incidence among adolescents not only leads to earlier disease onset but also prolongs the duration of exposure, thereby amplifying the risk of progressive liver disease and extrahepatic complications[5-7]. MASLD can progress from simple steatosis to inflammation, fibrosis, and cirrhosis[8], which underscores the urgency of effective interventions.

A growing body of research highlights chronic hepatic inflammation as a central factor in MASLD progression, primarily driven by Kupffer cell activation and elevated levels of pro-inflammatory cytokines such as interleukin (IL)-6 and tumor necrosis factor (TNF)-α[9,10]. Non-invasive markers, including uric acid to high-density lipoprotein cholesterol ratio, C-reactive protein, and systemic immune inflammation index have been proposed as diagnostic tools for MASLD[8,11]. Although several pharmacological treatments have been proposed, many clinical trials have not yielded positive results, and others are still ongoing. Consequently, lifestyle modifications remain the cornerstone of MASLD management[12]. Notably, exercise has been recognized as an effective modulator of both systemic inflammation and liver health. Physical activity has been demonstrated to lower pro-inflammatory cytokines such as TNF-α and IL-6, while restoring immune homeostasis, offering potential therapeutic benefits for MASLD[13,14]. Therefore, it is essential to further explore the mechanisms by which exercise improves MASLD.

Skeletal muscle, comprising about 40% of body weight in non-obese individuals, is crucial for maintaining systemic energy homeostasis and inflammation[15]. Myokines, such as IL-6[16], irisin[17], growth differentiation factor 11[18], and fibroblast growth factor 21 (FGF21)[19], secreted by skeletal muscle, have been shown to communicates with the other organs, such as liver[20] and adipose tissue[21], to regulate lipid and glucose metabolism. For example, irisin levels in circulation are lower in T2DM patients[22,23], while exercise-induced increase in fibronectin type III domain-containing protein 5 promotes irisin release, improving glucose tolerance and reducing insulin resistance[24,25]. These results emphasize the essential contribution of skeletal muscle to inflammation and metabolic regulation.

However, existing studies predominantly focus on individual myokines such as IL-6, irisin, and FGF21, leaving gaps in our understanding of the broader landscape of exercise-responsive muscle-derived factors involved in metabolic health. Notably, the molecular mechanisms by which skeletal muscle-derived factors mediate the protective effects of exercise in MASLD remain underexplored. In particularly, comprehensive transcriptional profiling of skeletal muscle response to exercise in MASLD remain unexamined.

Thus, the current study systematically analyzed skeletal muscle transcriptomic data from five gene expression omnibus (GEO) datasets, utilizing weighted gene co-expression network analysis (WGCNA), differentially expressed genes (DEGs) analysis, and least absolute shrinkage and selection operator (LASSO) methods. These analyses, conducted on pre-exercise (Con) and post-exercise (Train) muscle biopsies, identified key hub genes that were further validated in animal models to assess their relevance in MASLD.

By integrating multi-omics bioinformatics with experimental validation, this study aimed to identify exercise-responsive biomarkers in skeletal muscle that are mechanistically linked to MASLD. The novelty of our work lies in identifying candidate molecular targets involved in muscle-liver crosstalk, with focus on the exercise-induced adaptations of skeletal muscle and its mechanistic links to liver health. This study enhances our understanding of how muscle-derived factors orchestrate the benefits of exercise in MASLD and highlights the potential for precision diagnostic and therapeutic strategies targeting myokine-like candidates.

MATERIALS AND METHODS

Data collection

Sequencing data were obtained from the GEO database, including muscle biopsies from individuals with obesity and T2DM, both Con and Train. To ensure uniformity, we mapped the probes to their respective gene symbols using the platform annotations. The discovery groups were derived from the following datasets: GSE161749, GSE48278, GSE156247, and GSE53598. The discovery group was the GSE58249 dataset. Supplementary Table 1 provides comprehensive details of all datasets. These datasets served to construct co-expression networks and pinpoint key genes. All data procurement and usage adhered to the GEO database guidelines. Figure 1 illustrates the study design and corresponding analytical workflow.

Figure 1 Flowchart for identifying biomarkers related to exercise in the population with metabolic dysfunction-associated steatotic liver disease, including data extraction, procession and analysis. PCA: Principal component analysis; DEGs: Differentially expressed genes; WGCNA: Weighted gene co-expression network analysis; KEGG: Kyoto Encyclopedia of Genes and Genomes; GO: Gene Ontology; PPI: Protein-protein interaction; LASSO: Least absolute shrinkage and selection operator; ROC: Receiver operating characteristic; ELISA: Enzyme-linked immunosorbent assay; qPCR: Quantitative polymerase chain reaction.

Data preprocessing

The muscle biopsies datasets were merged into a single expression matrix, followed by adjustment for batch effects utilizing the ComBat function in the R sva package[26,27]. To assess the effectiveness of batch effect correction, we conducted principal component analysis using the prcomp function in R, followed by data visualization with the ggplot2 and RColorBrewer packages.

Differential expression analysis was performed utilizing limma and edgeR packages in R[28,29]. Due to the temporal and spatial variability of the samples in discovery groups, we set the threshold for differential gene expression to [|log₂ fold change (FC)| ≥ 0.26, adjusted P < 0.05] to maximize the identification of DEGs across the datasets. This threshold was selected based on standard practices for publicly available gene expression data, and aims to capture both biologically relevant changes and moderate expression shifts. Given the complexity and inherent variability in large-scale transcriptomic datasets, we felt that this threshold would be more appropriate for ensuring reproducibility while still capturing important exercise-responsive gene changes. The R package pheatmap was used to generate heatmaps for visual representation of the DEGs[30].

WGCNA

WGCNA was performed using the corresponding R package on the discovery groups to define key gene modules[31]. Initially, sample clustering was used to detect and exclude outlier samples, ensuring data quality for subsequent analysis. Trait data distinguishing Con and Train samples were processed to examine module-trait relationships. The co-expression network was established with an optimal soft-threshold power to achieve scale-free topology. Gene correlations were assessed through adjacency and topological overlap matrices. We then identified gene modules through hierarchical clustering and dynamic tree cutting, merging those with eigengene correlations exceeding 0.25. Associations between module eigengenes and clinical traits were evaluated, and significant modules were analyzed through calculations of gene significance (GS) and module membership (MM). In the analysis of module-trait correlations, genes exhibiting high hub modularity (GS > 0.20 and MM > 0.50) were identified as hub genes. Results were visualized using dendrograms and heatmaps. To ensure the robustness of these correlations, we performed subsampling analyses and examined their stability across datasets. The observed correlation coefficients (CC) remained relatively stable, supporting the consistency of module-trait relationships. To enable further analysis, all relevant data were compiled, including the refined gene matrix, module eigengene-trait relationships, and the calculated GS and MM values.

Identification of common hub genes and prognostic modeling

Hub genes shared between the DEG set and the WGCNA hub gene set were identified. A Venn diagram was drawn using the Venn package in R. Prognostic models for exercise-related hub genes were developed using LASSO regression analysis via the glmnet package[32]. K-fold cross-validation (k = 10) was employed to tune model parameters, and the optimal model was chosen based on the minimum partial likelihood deviance.

Functional enrichment analysis

Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis[33] and Gene Ontology (GO) analysis[34] were conducted to explore the biological significance of the identified genes and modules by identifying metabolic pathways. Additionally, the metascape platform[35] was employed to further characterize the distinctive attributes of the hub genes.

Construction and analysis of the protein interaction network

Using the online STRING platform, protein-protein interaction (PPI) networks were generated[36]. A cutoff value of interaction score > 0.70 was set in the STRING database. The gene with the highest score within subnetworks was chosen as a hub gene.

Validation of exercise biomarkers

Gene expression box plots for key biomarkers were created using the ggpubr package in R. Using the pROC package in R[37], we plotted receiver operating characteristic (ROC) curves and quantified the area under the curve (AUC) values to assess the diagnostic performance of the model in terms of its sensitivity and specificity[38].

Exercise and in vivo studies in mice

Male C57BL/6 mice (4 weeks old, body weight: 18 ± 2 g) were obtained from Changzhou Cavens Model Animal Co. Ltd. All animals were housed under specific pathogen-free (SPF) conditions at the Institute of Translational Medicine, Nanjing Medical University, within the SPF-grade transgenic animal center. The mice were kept in controlled conditions (22 ± 2 °C, 50%-60% humidity), with a 12 hours/12 hours light-dark cycle, and ad libitum access to food and water.

Following a 7-day acclimatization period, the mice received a high-fat diet (HFD) (60% kcal fat, D12492, research diets) for 14-weeks period to induce metabolic dysfunction. The minimum animal number was calculated using power analysis and in accordance with the 3Rs principles. Mice were randomly divided into two experimental groups (n = 6 per group): HFD with no exercise (HFD-Sed) and HFD with exercise (HFD-Train). Before the formal training regimen began, the mice were acclimated to the treadmill exercise over a 3-day period (60 minute/day at 8 m/minute). Following the acclimatization period, the mice underwent 10 weeks of running training on a treadmill at 12 m/minute, 1 hour/day, 5 days/week. Mice in the HFD-Sed group remained sedentary in their cages. During exercise training, mice were gently encouraged to continue running by tapping their tails and backs. All animal experiments were approved by the Ethics Committee for Animal Experiments at Affiliated Children’s Hospital of Jiangnan University (Ethical Review Consent No. WXCH2024-05-101). Mice were grouped using a random number table, and sample collection and evaluation were performed using a double-blind method. The supervisor was the only one aware of the group assignments at various stages of the experiment.

Metabolism measurement

Body weight was monitored regularly throughout the study. All mice were euthanized at 30 weeks of age, after 10 weeks of the experimental period. Serum, liver and gastrocnemius muscle tissues were harvested for subsequent analysis.

After collection, mouse blood was kept at 4 °C for 2 hour, then centrifuged at 3000 × g for 5 minutes (4 °C) to separate serum. The serum was used for both total cholesterol (TC) (E1005; Applygen, Beijing, China) and triglyceride (TG) (E1003; Applygen, Beijing, China) measurements. For TC measurement, the serum was diluted fourfold with anhydrous ethanol. Cholesterol standard solutions with concentrations of 2500 μmol/L, 625 μmol/L, 156 μmol/L, 39 μmol/L, and 0 μmol/L were prepared. Ten microliters of samples were added to the reaction tubes, followed by 190 μL of the working solution. After incubating at 37 °C for 20 minutes, the absorbance of the reaction mixtures was measured at 550 nm, and the TC concentration was determined using the standard curve. For TG measurement, the same procedure was followed, except that the serum was diluted with physiological saline.

For liver tissue, weigh 50 mg of tissue and homogenize it with lysis buffer. The mixture was initially incubated at room temperature (RT) for 10 minutes, and transfer the supernatant. Next, the supernatant was heated at 70 °C for 10 minutes, followed by centrifugation at 2000 rpm for 5 minutes. The final supernatant was collected for subsequent analysis. For TG measurement (E1013; Applygen), dilute the glycerol standard to final concentrations of 1000 μmol/L, 250 μmol/L, 62.5 μmol/L, 15.625 μmol/L, and 0 μmol/L using physiological saline. For TC measurement (E1015; Applygen), dilute the cholesterol standard to final concentrations of 2500 μmol/L, 625 μmol/L, 156 μmol/L, 39 μmol/L, and 0 μmol/L using anhydrous ethanol. Mix 10 μL of sample with 190 μL of working solution, incubate at 37 °C for 15-20 minutes, optical density at 550 nm was measured. We determined the supernatant protein concentration with a bicinchoninic acid (BCA) Protein Quantification Kit (E112-02; Vazyme, Nanjing, Jiangsu Province, China) and used to normalize hepatic TC and TG levels.

Real-time quantitative polymerase chain reaction

Total RNA was isolated from the gastrocnemius muscles of mice using the Universal RNA Purification Kit (EZB-RN4; EZBioscience, Roseville, MN, United States). One microgram of total RNA was reverse-transcribed into complementary DNA with the HiScript III RT SuperMix (R323-01; Vazyme) and then subjected to quantitative polymerase chain reaction (qPCR) on the SLAN-96S PCR system (HongShi Medical Tech, Shanghai, China) using ChamQ SYBR qPCR Master Mix (Q321-02; Vazyme). The 2^-ΔΔCt method was employed to analyze gene expression, using Actb for normalization. Primer sequences are listed in Supplementary Table 2.

Histological staining

Following fixation in 4% neutral paraformaldehyde for 24 hours, liver tissues were embedded in paraffin, sectioned to 4 μm, and subjected to hematoxylin and eosin (HE) staining for histological analysis. For routine HE staining, deparaffinized sections (using eco-friendly dewaxing solution I and II for 20 minutes each) were dehydrated through an ethanol gradient. After rinsing, the nuclei were stained with hematoxylin for 3-5 minutes, followed by cytoplasmic counterstaining with eosin for 15 seconds. The sections underwent differentiation, bluing, and dehydration before being mounted.

For liver lipid detection, frozen liver samples were sectioned into 5 μm slices, equilibrated at RT for 5 minutes, and then fixed for 15 minutes. Sections were incubated with freshly prepared oil red O solution for 15 minutes, and rinsed with 60% isopropanol. Hematoxylin counterstaining was performed for 1 minute, followed by alcohol-acetic acid differentiation for 1-5 seconds. Sections were rinsed with running water, then distilled water, and sealed with glycerol and gelatin. All slides were digitized using a pathology slide scanner (KFBIO KF-PRO-120; Ningbo, Zhejiang Province, China) and visualized under 20 × magnification using K-Viewer software.

Western blotting

Gastrocnemius muscle tissues were lysed in radio immunoprecipitation assay buffer (P0013B; Beyotime) with protease and phosphatase inhibitors (P1049; Beyotime). Then tissues were homogenized using a pre-cooled grinding machine (MM400; RETSCH, Germany), followed by incubation (20 minutes, 4 °C) and centrifugation (13000 × g, 30 minutes, 4 °C). After BCA quantification (E112-02; Vazyme), twenty milligrams of protein per sample was separated by 4%-20% sodium dodecyl sulfate-polyacrylamide gel electrophoresis gel (Tanon, Shanghai, China) for electrophoresis, and transferred to polyvinylidene difluoride membranes (ISEQ00010; Merck Millipore, Burlington, Germany). Membranes were blocked with 5% non-fat milk in tris buffered saline tween (TBST) (1 hour, RT) and incubated overnight (4 °C) with primary antibodies specific to LAMA4 (10465-1-AP; 1:1000, Proteintech, IL, United States), PECAM1 (85898-4-RR; 1:1000, Proteintech), THBS4 (ab263898; 1:1000, abcam, MA, United States), PXDN (abs153157; 1:1000, Absin, Shanghai, China) and β-actin (81115-1-RR; 1:30000, Proteintech). Following three TBST washes, membranes were incubated with horseradish peroxidase-conjugated secondary antibodies (1:10000 in TBST) for 1 hour at RT. Chemiluminescent signals were detected using Tanon 5200 multi-imaging system (Tanon), with β-actin as the internal loading control.

Enzyme-linked immunosorbent assay

Mouse PECAM1 enzyme-linked immunosorbent assay (ELISA) kit (EM0155; Fine Test, Wuhan, Hubei Province, China) and THBS4 ELISA kit (MM-46905M1; Meimian; Yancheng, Jiangsu Province, China) were used to measure PECAM1 and THBS4 levels in mouse serum. All analyses were conducted following the manufacturers’ instructions.

Statistical analysis

Statistical analysis and data visualization were performed using GraphPad Prism 10.1.2. Results are presented as the mean ± SEM. Before applying statistical tests, the normality of the data was assessed using the Shapiro-Wilk test. Statistical significance was evaluated using a two-tailed independent Student’s t-test. For experiments involving multiple comparisons, the Benjamini-Hochberg procedure was applied to correct for false discovery rates. Significance levels are indicated as follows: ^aP < 0.05, ^bP < 0.01, ^cP < 0.001, and not significant.

RESULTS

Bioinformatics workflow

As outlined in Figure 1, we gathered and normalized four GEO datasets to create a combined dataset. We conducted differential expression analysis to identify DEGs, whose biological functions and pathway associations were subsequently interrogated through GO and KEGG enrichment analyses. Additionally, a PPI network was constructed based on the identified DEGs. Subsequently, we applied WGCNA to build a co-expression network and pinpoint potential hub genes, followed by an integrated analysis combining both DEGs and hub genes to select candidate genes. To validate the findings, expression validation was carried out using a validation cohort to assess the diagnostic efficacy of the hub genes. Finally, an obese mouse model was established and subjected to endurance aerobic training, with gastrocnemius muscle tissue selected for gene and protein expression verification.

Identification of DEGs

To profile exercise-induced transcriptional changes in muscle tissue of individuals with metabolic disorders before and after exercise, we integrated four independent datasets (GSE161749, GSE48278, GSE156247 and GSE53598) from the GEO database, encompassing exercise intervention studies across diverse metabolic disorder populations. Comprehensive details for all datasets are available in Supplementary Table 1.

To mitigate batch effects across platforms, we applied the ComBat algorithm for data normalization. Principal component analysis confirmed the normalization of GEO samples across the four datasets (Figure 2A and B). We identified DEGs via the limma package in R. A threshold of |log2FC| ≥ 0.26 and P < 0.05 was applied, resulting in the identification of 61 DEGs, including 40 upregulated genes and 21 downregulated genes, whose distribution is visualized in volcano plots and heatmaps (Figure 2C and D). The log₂FC, average expression across all samples, and adjusted P values after multiple testing correction are detailed in Supplementary Table 3. These DEGs effectively distinguished Train from Con muscle samples, as illustrated in the figures.

Figure 2 Identification of differentially expressed genes in muscle tissue samples. A and B: Principal component analysis of different gene expression omnibus datasets. The gene expression profiles were visualized using scatter plots based on the first two principal components with (B) and without (A) batch effect correction; C: Volcano plots of the differentially expressed genes (DEGs); D: Heatmap showing DEGs between pre-exercise and post-exercise samples. PC1: Principal components 1; PC2: Principal components 2; FC: Fold change; Sig: Significance; Down: Down-regulation; Not: Not significant; Up: Up-regulation; Con: Pre-exercise; Train: Post-exercise.

Functional enrichment analyses of DEGs

To identify the biological roles of the 61 DEGs, we carried out both GO functional and KEGG pathway analyses. GO analysis revealed significant enrichment in biological processes, such as “extracellular structure organization”; cellular components, including “collagen-containing ECM”; and molecular functions such as “ECM structural constituent” (Figure 3A). KEGG analysis highlighted enrichments in “cytoskeleton in muscle cells” and “focal adhesion” (Figure 3B). See Supplementary Tables 4 and 5 for the complete GO terms and KEGG pathways. These results suggest the DEGs are involved in ECM organization and collagen remodeling. Additionally, a PPI network analysis identified 17 DEGs with significant interactions (Figure 3C), including COL1A1, COL1A2, COL3A1, COL4A1 and COL4A2, and others, underscoring their relevance in MASLD pathophysiology.

Figure 3 Functional enrichment analysis of differentially expressed genes. A: Bubble charts showing the enrichment of Gene Ontology terms among differentially expressed genes (DEGs); B: Bubble charts depicting the enrichment of Kyoto Encyclopedia of Genes and Genomes pathways among DEGs; C: Protein-protein interaction network diagram for 17 DEGs. BP: Biological processes; CC: Cellular components; MF: Molecular functions; ECM: Extracellular matrix; PI3K: Phosphatidylinositol 3-kinase; Akt: Protein kinase B; AGE-RAGE: Advanced glycation end products-receptor for advanced glycation end products.

WGCNA

After confirming the absence of missing values in the dataset, we conducted WGCNA to construct co-expression modules, with a soft threshold (β) to 7, which ensured a scale-free network (fit index > 9) (Figure 4A). Dynamic tree cut method was applied to amalgamate modules with similar patterns (Figure 4B). The heatmap of module-trait correlations (Con and Train) is shown in Figure 4C. The MEmediumpurple3 module showed the highest positive correlation with Train muscle samples (CC = 0.15, P = 0.02), while the MEpowderblue module exhibited a strong negative correlation (CC = 0.2, P = 0.003) with Train muscle samples (Figure 4D and E).

Figure 4 Weighted gene co-expression network analysis of differentially expressed genes. A: Scale-free topology fitting index (left) and mean connectivity (right) across soft-thresholding powers (β), with the red line indicating the threshold correlation coefficient of 0.9; B: Gene clustering dendrogram constructed using the topological overlap matrix, with color bars representing an co-expression modules; C: Heatmap illustrating the correlation between module eigengenes and clinical traits in metabolic dysfunction-associated steatotic liver disease subjects before and after exercise, with positive correlations in red and negative in blue; D: Scatter plot of gene significance vs module membership in the mediumpurple3 module; E: Scatter plot of gene significance vs module membership in the powder blue module. Con: Pre-exercise; Train: Post-exercise.

Identification and validation of common hub genes

We intersected the DEGs with WGCNA results, identifying 14 key genes (Figure 5A). LASSO regression screening yielded potential 12 hub genes: COL3A1, COL1A2, COL14A1, BGN, LAMA4, THY1, THBS4, PAMR1, PXDN, LAMB1 and ASPN (Figure 5B and C). GO analysis performed using the metascape platform, revealed enrichment in pathways such as the “Naba core matrisome”, and “ECM proteoglycans” (Figure 5D). DisNET analysis further highlighted strong associations with “familial thoracic aortic aneurysm and aortic dissection” and “aortic aneurysm abdominal” (Figure 5E). Additionally, molecular complex detection analysis enriched for collagen chain trimerization (Figure 5F). Finally, the PPI network revealed significant interactions among seven genes (BGN, COL1A2, COL3A1, ASPN, COL14A1, LAMA4 and LAMB1) (Figure 5G).

Figure 5 Identification and functional enrichment analysis of key hub genes. A: Venn diagram showing the overlap of hub genes from differentially expressed genes and weighted gene co-expression network analysis; B and C: Identification of twelve genes through least absolute shrinkage and selection operator regression; curves represent individual genes, with the dashed line marking the optimal lambda; D: Gene Ontology (GO) enrichment analysis of hub genes, with -log₁₀ P values on the X-axis; E: DisGeNET enrichment of hub genes, X-axis showing -log₁₀ P values; F: Top molecular complex detection clusters of hub genes from the protein-protein interaction (PPI) network, with GO analysis highlighting key biological functions; G: PPI network of the 7 key hub genes. WGCNA: Weighted gene co-expression network analysis; GO: Gene ontology; MCODE: Molecular complex detection.

To assess the relationship between 12 hub genes with exercise, we validated these genes using the GSE58249 cohort (GSE58249). The expression patterns of COL3A1, COL1A2, BGN, LAMB1, PECAM1, LAMA4, THBS4, PXDN, and THY1 were consistent with the discovery cohort (Figure 6 and Supplementary Figure 1). ROC curve analysis was used to assess the diagnostic potential of the 12 hub genes, the area under the AUC for these nine genes (THY1, BGN, COL1A2, COL3A1, LAMA4, LAMB1, PECAM1, PXDN, and THBS4) demonstrated strong diagnostic performance (Figure 7A-I and Supplementary Figure 2A-I). In contrast, COL14A1, ASPN and PAMR1 showed no significant diagnostic ability (Figure 7J-L and Supplementary Figure 2J-L). In conclusion, these findings suggest that 9 hub genes demonstrate strong predictive accuracy and may serve as valuable biomarkers.

Figure 6 Expression of key hub genes and validation in validation group. A-L: Box plots illustrating the expression levels of genetic biomarkers in the validation cohort (GSE58249): THY1 (A); BGN (B); COL1A2 (C); COL3A1 (D); LAMA4 (E); LAMB1 (F); PECAM1 (G); PXDN (H); THBS4 (I); ASPN (J); COL14A1 (K); PAMR1 (L). ^bP < 0.01. ^cP < 0.001. NS: Not significant; Con: Pre-exercise; Train: Post-exercise.

Figure 7 Diagnostic prediction efficacy analysis of key hub genes and validation in validation group. A-L: Receiver operating characteristic curve analysis for the validation cohort (GSE58249): THY1 (A); BGN (B); COL1A2 (C); COL3A1 (D); LAMA4 (E); LAMB1 (F); PECAM1 (G); PXDN (H); THBS4 (I); ASPN (J); COL14A1 (K); PAMR1 (L). AUC: Area under the curve; CI: Confidence interval.

Establishment of a metabolic dysfunction model combined with endurance exercise and validation of genetic biomarker expression

To validate the bioinformatics analysis results, we performed an in vivo animal experiment. Obesity is a major driver of MASLD; therefore, mice were administered an HFD for 14 weeks to induce MASLD. After 10 weeks of treadmill training, mice in the HFD-Train group exhibited a significant reduction in body weight compared to the HFD-Sed group (Figure 8A). Exercise markedly decreased serum TC and TG levels, as well as hepatic TC and TG levels (Figure 8B-E). Histological analysis revealed markedly reduced liver inflammation and lipid deposition following exercise (Figure 8F).

Figure 8 Expression of genetic biomarkers in high-fat diet-induced metabolic dysfunction-associated steatotic liver disease combined with endurance exercise. A: Body weight changes in mice subjected to 14 weeks of high-fat diet (HFD) followed by 12 weeks of endurance exercise. Measurements were taken before and after the exercise intervention, n ≥ 5; B and C: Serum total cholesterol (TC) (B) and triglyceride (TG) (C) levels in HFD mice with no exercise (HFD-Sed) and HFD with exercise (HFD-Train), n ≥ 5; D and E: Liver TC (D) and TC (E) levels in mice with HFD-Sed and HFD-Train group, n ≥ 5; F: Representative images of hematoxylin and eosin staining and oil red O staining of liver sections. Scale bars, 100 μm; G: The messenger RNA levels of Lama4, Pecam1, Pxdn and Thbs4 in gastrocnemius muscle of HFD-Sed and HFD-Train mice were measured by quantitative polymerase chain reaction, Actb was used as an internal control, n ≥ 5; H-K: Western blot analysis to assess the protein expression levels of LAMA4 (H), PECAM1 (I), PXDN (J), THBS4 (K), n ≥ 4; L and M: The serum level of PECAM1 (L) and THBS4 (M) were measured using an enzyme-linked immunosorbent assay, n ≥ 5. ^aP < 0.05. ^bP < 0.01. ^cP < 0.001. HFD: High-fat diet; HFD-Sed: High-fat diet with no exercise; HFD-Train: High-fat diet with exercise; TC: Total cholesterol; TG: Triglyceride; HE: Hematoxylin and eosin; mRNA: Messenger RNA.

Total RNA was extracted from gastrocnemius muscle and subjected to qPCR analysis. Among the candidate genes identified, exercise intervention markedly increased the expression of Lama4, Pecam1 and Pxdn, however, Thbs4 messenger RNA expression was significantly reduced after exercise (Figure 8G). However, the messenger RNA expression of Bgn, Lamb1, Thy1, Col3a1, and Col1a2 showed no significant differences (Supplementary Figure 3). The HFD-Train group exhibited significantly higher protein expression levels of LAMA4, PECAM1 and PXDN compared to HFD-Sed group, whereas THBS4 expression was lower (Figure 8H-K). ELISA results further supported these findings, showing elevated serum PECAM1 levels and reduced THBS4 levels in the HFD-Train group (Figure 8 L and M). In summary, the expression validation confirmed that LAMA4, PECAM1 and PXDN were upregulated after exercise, consistent with their predicted roles in exercise-mediated metabolic regulation, while THBS4 was downregulated.

This differential expression profile indicates that while most candidate genes may be involved in the benefits of exercise for MASLD through upregulation, THBS4 may downregulated, implying a distinct role in the transcriptional response to exercise adaptation.

In conclusion, these findings demonstrate that endurance exercise promotes the upregulation of LAMA4, PECAM1, and PXDN, reinforcing their relevance as biomarkers in exercise-induced metabolic adaptation and MASLD improvement. By contrast, the suppression of THBS4 highlights an alternative transcriptional program, suggesting its potential involvement in a separate regulatory axis of exercise-mediated protection against MASLD.

DISCUSSION

Exercise is a well-recognized non-pharmacological intervention with significant benefits for improving MASLD. However, the biological mechanisms through which skeletal muscle responds to regular exercise and contributes to MASLD improvement are still not well understood. By integrating multi-omics bioinformatics with experimental validation, this study identified four exercise-responsive muscle-expressed genes (LAMA4, PECAM1, PXDN and THBS4). These genes are mechanistically associated with MASLD and our detection of secreted PECAM1 and THBS4 in circulation suggests that they function as myokine-like candidates. The results elucidate the molecular adaptations in skeletal muscle induced by exercise and highlight their links to liver health.

We first conducted differential expression analysis to identify DEGs and performed functional enrichment using GO and KEGG pathway analyses to uncover their associated biological functions. Next, a co-expression network was constructed to identify hub genes, which were further refined through integration with DEGs. The diagnostic potential of these hub genes was validated in an external cohort. In addition, we established an obese mouse model subjected to endurance aerobic training, collecting gastrocnemius muscle tissue to validate gene expression.

Our transcriptomic analysis identified 61 DEGs between muscle samples taken before and after exercise, including 40 upregulated and 21 downregulated genes. GO analysis revealed significant enrichment in pathways related to collagen fibril organization and extracellular matrix (ECM), emphasizing the role of ECM in muscle development[39], growth, repair[40] and the transmission of contractile forces[41].

KEGG pathway analysis further highlighted the enrichment of pathways associated with cytoskeleton dynamics in muscle cells and ECM-receptor interactions. Combining WGCNA with LASSO, we identified 12 hub genes. Subsequent ROC analysis revealed that nine of these hub genes, including COL3A1, COL1A2, BGN, LAMB1, PECAM1, LAMA4, THBS4, PXDN and THY1, demonstrated strong diagnostic performance. In our mouse model, we observed upregulation of Lama4, Pecam1 and Pxdn, along with downregulation of Thbs4 following exercise.

Laminin isoforms, a key component of the basement membrane, are heterotrimeric structures consisting of α, β, and γ subunits[42]. LAMA4, belonging to the laminin family, is predominantly localized in the basement membrane[43]. These isoforms crucial roles in forming membrane structures and signaling pathways, which vary depending on tissue type and disease state[42,44]. LAMA4 has been shown upregulated during adipogenesis and regulates the basement membrane in adipose tissue[45]. Knockout studies in mice have revealed improved phenotypes and metabolism in the absence of LAMA4[44,46,47]. However, no association has been found between low LAMA4 expression and metabolically healthy obesity[48]. Previous studies have investigated the role of LAMA4 in collagen composition and ECM remodeling in adipose tissue. There is limited information on whether LAMA4 response to exercise in muscle tissue under MASLD conditions. Our study contributes to this gap by demonstrating that LAMA4 expression in skeletal muscle responds to endurance exercise under MASLD conditions. These findings suggest a potential muscle-liver crosstalk mechanism mediated by ECM remodeling, providing a possible target for therapeutic intervention in MASLD.

PECAM1 is a transmembrane protein in the immunoglobulin superfamily. It acts as both an endothelial cell adhesion molecule and mechanoreceptor, facilitating interactions between neighboring endothelial cells. PECAM1 is involved in regulating inflammation, the migration of leukocytes, and vascular responses, especially during sepsis[49]. Additionally, PECAM1 can be cleaved to release a soluble form in response to shear stress. This soluble form is secreted into the circulation and contributes to amplifying inflammatory responses[49-51]. Transcriptomic analysis has shown reduced PECAM1 expression in the adipose tissue of obese individuals[52]. Exercise has been found to enhance PECAM1 expression in circulating angiogenic cells[53]. In our MASLD model with exercise intervention, we observed significant upregulation in both the expression and secretion PECAM1. These findings suggest that PECAM1 functions as a myokine-like candidate, contributing to the vascular and metabolic benefits of exercise in MASLD. Further research is required to elucidate its role and clinical utility.

PXDN, a secreted heme peroxidase, facilitates the oxidative cross-linking of collagen IV in the ECM, generating hypobromous acid through the reaction of hydrogen peroxide and bromide[54]. PXDN is mainly expressed in cardiovascular tissues and is released into the bloodstream, where it contributes to the progression of cardiovascular diseases[55]. In endothelial cells, PXDN has been shown to promote cell death by inducing apoptosis[56] and programmed necrosis[57]. Elevated PXDN expression levels have been noted in the aorta of T2DM rats[58], and its deletion has been linked to early onset of obesity[59,60]. Our bioinformatics analysis and experimental validation further reveal that PXDN is significantly overexpressed in the HFD-Train group. This upregulation after exercise suggests a potential role for PXDN in muscle-derived regulation of systemic metabolism in MASLD. However, the precise role of PXDN in lipid and glucose metabolism remains unclear and warrants further investigation to fully understand its impact.

In contrast to the upregulation of LAMA4, PECAM1 and PXDN, our analysis revealed a significant downregulation of Thbs4 expression after exercise intervention. THBS4 encodes thrombospondin 4, a matricellular glycoprotein involved in ECM organization, cell–matrix interactions, and tissue remodeling[61,62]. THBS4 has been implicated in wound healing and tissue repair by modulating matrix remodeling[63]. It has also been linked to the pathogenesis of cardiovascular disorders, tumor progression, and neurodegenerative diseases, highlighting its broad relevance in disease biology[64,65]. In liver cancer, THBS4 promotes tumor growth, metastasis and epithelial mesenchymal transition through interactions with integrin β1, activating the focal adhesion kinase/phosphatidylinositol 3-kinase/protein kinase B signaling pathway[66]. Elevated THBS4 levels in hepatocellular carcinoma (HCC) correlate with poor prognosis[67]. Additionally, THBS4 has been suggested as an early-stage biomarker for HCC[68]. As a secreted factor, THBS4 plays a role in inter-organ communication and may influence systemic metabolic balance. Its downregulation after exercise could be beneficial, highlighting THBS4 as a potential therapeutic target for improving MASLD.

Our findings highlight the integration of bioinformatics and experimental validation in identifying candidate biomarkers and therapeutic targets for MASLD. We acknowledge that the relatively permissive |log₂FC| ≥ 0.26 threshold, chosen to maintain sensitivity in detecting relevant gene expression changes, may influence the number of DEGs identified. While applying stricter thresholds (e.g., |log₂FC| ≥ 0.5) would reduce the number of DEGs, it did not significantly alter the overall biological patterns observed. The threshold selection reflects standard practices for analyzing publicly available GEO datasets and can be refined in future studies based on specific contexts or validation datasets. The modest module-trait CCs (CC: 0.15-0.2) observed, though statistically significant, are typical in large, high-dimensional datasets where many factors influence gene expression. These modest correlations still hold substantial biological relevance, especially in complex traits like exercise-induced muscle adaptations, where even small shifts in gene expression can reflect important underlying processes. Notably, these correlations remained stable across different datasets and subsampling, reinforcing the reliability of the identified gene modules. It is important to clarify that the ROC analyses performed in this study were based on transcriptomic data from the mouse model, rather than clinical diagnostic tests.

Our study had several limitations. First, it relied on transcriptomic data from public databases, and while cross-validation was performed, batch effects and platform inconsistencies may still affect generalizability. While the animal model used showed promising findings, clinical validation in human cohorts is essential to confirm their relevance. Additionally, we did not evaluate glucose tolerance, insulin resistance, or HOMA-type indices, which limits our understanding of how exercise influences lipid and glucose homeostasis in MASLD, and the mechanisms by which the identified exercise-responsive genes regulate these processes remain unclear. Furthermore, while we propose these genes as potential biomarkers of exercise response, their expression in human serum or muscle biopsies in MASLD has not been fully explored. The next steps involve validating these biomarkers in human cohorts, conducting proteomic analyses, and performing functional assays to establish causality. Overall, these efforts will bridge the gap between transcriptomic discoveries and clinical relevance, offering new avenues for MASLD treatment. Future investigations should focus on human validation, proteomic studies and longitudinal trials to confirm the therapeutic potential of these biomarkers.

CONCLUSION

Our study identified key exercise-responsive biomarkers in skeletal muscle that are mechanistically linked to MASLD, including LAMA4, PECAM1, PXDN and THBS4. By integrating bioinformatics analysis with experimental validation, we have uncovered new insights into the molecular adaptations in skeletal muscle induced by exercise and their role in improving liver health. Our results position exercise-responsive candidates as both therapeutic targets and diagnostic biomarkers for MASLD. Future research should prioritize validating these biomarkers in human cohorts, conducting proteomic analyses, and performing functional assays to establish causality, which will ultimately facilitate the development of more effective treatment strategies for MASLD.