Male C57BL/6J mice were purchased from Shanghai Research Center for Model Organisms, and all experimental procedures and animal care were following the guidelines of the Ethics Committee of Animal Experiments of East China Normal University. The mice were maintained under specific-pathogen-free conditions at 22°C as the standard housing temperature and under a 12-h light/dark cycle with free access to food and water. The NAFLD group of mice was fed with high-sucrose and high-fat (HSHF) diet with 40 kcal% fat, 20 kcal% fructose, and 2% cholesterol (D09100310, Research Diets) for 12 weeks, and the NASH group of mice was fed with HSHF diet for 28 weeks, while the normal group was fed with control chow diet (Medicience Ltd., MD17111). The mice were anesthetized with sodium pentobarbital (Nembutal, 80 mg/kg, ip) and sacrificed. Samples, such as the liver and quadriceps muscles, were harvested for further analysis.

Formalin-fixed, paraffin-embedded mouse liver sections were stained. Liver fibrosis was assessed by Sirius red (Polysciences, catalog 24901) staining as per the manufacturer’s instructions. Frozen sections from NAFLD and NASH mice were stained with 0.5% oil red O staining reagent for 20 min and were counterstained with hematoxylin for 5 min, and H&E staining was performed for the quantitative analysis of lipid accumulation in the liver. The different histological features of NAFLD and NASH were assessed using the NASH Clinical Research Network Scoring System as described previously (22).

Quadriceps muscles from NAFLD, CON1, NASH, and CON2 groups of mice were collected and frozen for subsequent RNA extraction. We have randomly selected 3 quadriceps muscle samples from each group of mice for RNA sequencing and 6 samples for qPCR confirmation. Total RNA was extracted from the liver or quadriceps muscle tissues using RNAiso Plus (Takara, 9108), following the manufacturer’s instructions, and the concentration was determined via NanoDrop Microvolume Spectrophotometers and Fluorometer (Thermo Fisher Scientific).

The RNA quality was checked using Bio-analyzer instrument (Agilent, USA), quantified with ND-2000 (NanoDrop Technologies), and then subjected to the construct cDNA library (CloudSeq Biotech Inc., Shanghai, China). Briefly, Ploy A mRNA was enriched using Oligo dT magnetic beads from 2 µg total RNA and broken up to 200 bp for each replicate. The double-strand cDNA was synthesized and purified for end repair, poly A addition, and adapter ligation. Then, the products enriched with PCR amplification generated cDNA libraries that were subsequently sequenced on Illumina HiSeq™ 3000 platform. The RNA-seq raw reads generated by Illumina sequencer were performed by trimming of adaptors and removal of low-quality reads using Skewer (v0.2.2). FastQC (v0.11.5) was used to check the quality of the pre-treated data which were mapped to GRCm38 using STAR (2.5.3a). The transcripts were assembled using StringTie (v1.3.1c), and the differential gene transcript expression between the early- and late-phase groups was analyzed with DESeq2 (v1.16.1). The differential threshold value was P-value ≤0.05 and fold-change ≥1.5.

For the qPCR analysis, 1 μg of total RNA was reverse-transcribed as cDNA by PrimeScript™RT Master Mix (Takara, RR036A). The qPCR analysis was performed on quantitative real-time PCR system (Roche, LightCycler 480) with SYBR Green Master Mix (Thermo Fisher Scientific, 4309155). Internal control using Gapdh gene for data analysis and cycle threshold (Ct) values was calculated using the 2^-ΔΔCt method. The primer sequences used in this study are shown in Supplementary Table S1 .

Protein extraction was performed with RIPA buffer containing 50 mM Tris (pH 7.4), 150 mM NaCl, 1% Triton X-100, 1% sodium deoxycholate, 0.1% SDS, sodium orthovanadate, sodium fluoride, EDTA, leupeptin, supplemented with 1 mM PMSF, 10 mM DTT, and 10 μM protein kinase inhibitor. The protein concentrations were quantified using a BCA protein quantification kit (Beyotime), and the samples from the SDS-PAGE gels were transferred to a membrane of nitrocellulose by electrophoretic transfer. The membranes were blocked with 5% skimmed milk for 1 h at room temperature and then incubated with anti-MSTN (Santa Cruz, sc-393335) and anti-GAPDH (Aksomics Biotechnology, KC-5G4) antibodies. The membranes were washed with Tris-buffered saline with Tween three times and incubated with IRDye secondary antibodies (926-68071 and 926-322210) from LI-COR Biosciences (NE, USA). The protein bands were detected by Odyssey^® CLX imaging system (LI-COR Biosciences, NE, USA).

Database for Annotation, Visualization, and Integrated Discovery Functional Annotation Tool (23, 24) was used to perform Gene Ontology (GO) terms, including biological process (BP), cell component, and molecular function, as well as Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis. The top altered GO terms and KEGG pathways are shown in Supplementary Table S2 .

The protein–protein interaction (PPI) network was built by using the Search Tool for the Retrieval of Interacting Genes (http://string.embl.de/) (25). The uploaded genes were clustered into networks to detect significant functional modules. Cytoscape (26) was used to subsequently visualize the PPI network and to calculate the properties of the network based on the interaction pair information.

All the data are represented as mean ± SEM. The data analysis was performed using GraphPad Prism7. The data normality test was examined with Shapiro–Wilk normality test. Two-tailed t-test analysis was performed to identify differences between two groups. P-value <0.05 was considered statistically significant.

It has been reported that the development of NAFLD and NASH is accompanied with skeletal muscle dysfunctions (27). To understand the detailed molecular changes of skeletal muscles in NAFLD and NASH, we established NAFLD and NASH mice models by feeding mice with HSHF diet (40 kcal% fat, 20 kcal% fructose, and 2% cholesterol) for 12 weeks (as the NAFLD group) or 28 weeks (as the NASH group) as previously reported (28, 29). Their control group of mice was fed with chow diet (Medicience Ltd., MD17111) for 12 weeks (CON1) or 28 weeks (CON2) separately ( Figure 1A ). A detailed analysis on liver morphology revealed increased lipid accumulation in the livers of both the NAFLD and NASH groups of mice, in accordance with increased triglyceride contents, and a strong fibrotic area in the livers of the NASH group as shown by hematoxylin–eosin (H&E), oil red O, and Sirius red staining ( Figures 1B–E ).

The QU muscle is one of the largest skeletal muscles mixed with oxidative and glycolytic muscle fibers. Histological analysis and oil red O staining showed increased lipid deposition, which was confirmed by increased triglyceride contents in the QU muscle in the NAFLD and NASH groups of mice ( Figures 1F–I ). To comprehensively understand the molecular signature of the skeletal muscle in NAFLD and NASH, compared to their age-matched controls, we further took QU muscle samples for RNA sequencing ( Figure 2 ).

Via RNA sequencing, we firstly identified 140 (54 upregulated and 86 downregulated) differentially expressed genes in the NAFLD group of mice as shown by a volcano plot ( Figure 3A ). The GO analysis of BP revealed that cytokine production and secretion, protein secretion, and regulation of hormone level were increased ( Figure 3B ), suggesting that a crosstalk between the muscle and the metabolic organs happened in NAFLD, consistent with a previous finding that blocking cytokines, such as IL6, in NAFLD mice could alleviate insulin resistance in the skeletal muscle (30, 31). In addition, the KEGG pathway analysis demonstrated that the PPAR signaling pathway, fatty acid metabolism, and AMPK signaling pathway were upregulated ( Figure 3D ), which were consistent with previous reports claiming that these pathways play critical roles in lipid and glucose metabolism on NAFLD (32). Besides these, muscle morphogenesis, muscle development, muscle contraction, and hypertrophic pathways were reduced in the QU muscle from the NAFLD group compared with the control group, as shown by GO and KEGG analysis ( Figures 3C, E ), suggesting the decline of muscle functionality. Interestingly, these declined genes are also shared by the cardiac system. Indeed NAFLD patients are reported to be associated with inflammation-related atrial and ventricular myopathy, which is the predictor of coronary heart disease (33, 34). Notably, we validated the changes of the top altered gene using qPCR and found that Fasn and Scd1 were increased and that Tpm3, Myl2, Myl3, Myh7, Tnnt1, and Tnni1 were reduced in the NAFLD group compared with the controls ( Figures 3F, G and Supplementary Table S2 ). Overall, these data suggested that NAFLD is associated with increased lipid metabolism and inflammation, while it suppressed muscle functionality in the skeletal muscle.

Similarly, via RNA sequencing, data analysis identified 141 differentially expressed genes, including 53 upregulated and 88 downregulated genes, in NASH group of mice as shown in a volcano plot ( Figure 4A ). The GO enrichment analysis of biological process suggested that the regulation of lipid metabolic process and carbohydrate catabolic process was enriched to be upregulated ( Figure 4B ). Consistently, the KEGG pathway analysis revealed that the PPAR signaling pathway, fatty acid metabolism, and AMPK signaling pathway were dramatically enriched in upregulated genes in the skeletal muscle of the NASH group of mice ( Figure 4D ). Besides this, the GO analysis revealed that muscle tissue morphology, muscle development, and muscle contraction were downregulated in the NASH group muscle, similar to NAFLD ( Figure 4C ), which was consistent with the KEGG analysis showing that muscle contraction, adrenergic signaling, and hypertrophic genes were reduced ( Figure 4E ). Besides this, we confirmed the top altered genes by qPCR and found that Fasn, Scd1, Adipoq, Lep, and Srebf1 were largely increased, while Tpm3, Myl2, Myl3, Myh7, Tnnt1, and Tnni1 were sharply reduced in the NASH group in comparison with the control group ( Figures 4F, G and Supplementary Table S2 ). Overall, these results suggested an increase in lipid and carbohydrate metabolism and a reduction of muscle development and functions in the skeletal muscle in NASH.

We next set out to explore the common altered gene programs in both the NAFLD and NASH groups. Overlapping analysis identified 84 genes as displayed by a Venn diagram ( Figure 5A ). Among them, 19 upregulated and 56 downregulated genes were analyzed. In accordance with previous data, the GO and KEGG analysis showed that 19 upregulated genes were associated with fatty acid and lipid metabolic processes ( Figure 5B ) as well as the PPAR signaling pathway, fatty acid metabolism, and AMPK signaling pathway ( Figure 5D ). Besides this, 56 downregulated genes were related to muscle contraction and development as well as to muscle contraction, adrenergic signaling, and hypertrophic gene programs ( Figures 5C, E ). Furthermore, we examined protein–protein interactions for upregulated and downregulated genes with a PPI network ( Figure 5F ) and found that lipid synthetic genes like Fasn, Scd1, and Scd2 showed dramatic increases, while muscle fiber markers including Tnni1, Tnnt1, Myl3, and Myh7 showed sharp decreases. Overall, these results suggested increased lipid metabolism, majorly lipid synthetic genes, and reduced muscle developmental and functional genes that were commonly regulated genes in the skeletal muscle of both NALFD and NASH mice.

Next, we investigated the heterogeneity in gene expression patterns in the skeletal muscle of NAFLD and NASH mice. For 280 DEGs uniquely changed in the skeletal muscle of the NAFLD group of mice ( Figure 6A ), the GO and KEGG analysis reveal that NF-kB transcription factor activity and apoptotic processes were uniquely upregulated ( Figures 6B, D ), suggesting inflammatory and apoptotic events that happen in the early stage of liver diseases as reported previously (35); meanwhile, chronic aerobic exercise could ameliorate it (36). Besides this, we found that response to toxic substance, regulation of growth, angiogenesis, axis specification as well as HIF-1 signaling pathway, mineral absorption, PI3K-Akt signaling pathway were enriched in NAFLD uniquely downregulated genes ( Figures 6C, E ), suggesting that the skeletal muscle may exhibit a disturbance in response to a toxic substance and muscle growth. Furthermore, the PPI network diagram revealed that NF-kB and transforming growth factor-β signaling pathway, including Tab2, Map3k7, Traf3, Tgfβ3, and Tgfβ2, were highly induced ( Figure 6F ).

Moreover, 271 DEGs were uniquely changed in the skeletal muscle in the NASH group of mice ( Figure 7A ). Notably, lipid metabolic process as well as AMPK signaling pathway, PPAR signaling pathway, and insulin signaling pathway were uniquely upregulated in NASH mice as shown by the GO and KEGG analysis ( Figures 7B, D ), suggesting enforced lipid deposition and enhanced lipid and glucose metabolic dysfunctions in the skeletal muscle of NASH than the NAFLD group, while lipoprotein metabolic process as well as tight junction, MAPK signaling pathway, and focal adhesion were uniquely downregulated in the muscle from the NASH group ( Figures 7C, E ). In addition, the PPI network analysis revealed increased levels of Adipoq, Srebf1, Lep, Fabp4, and Cebpa, which were associated with increased lipid deposition, in the NASH group of mice ( Figure 7F ).

Lastly, to understand the possible crosstalk between muscle and other metabolic organs, including liver, we analyzed commonly altered myokines in skeletal muscle from the NAFLD and NASH groups of mice compared to their controls. Notably, 4 upregulated myokines (Mstn, Stc2, Gdf11, and Ostn) and 8 downregulated myokines (Serpina3k, Serpina1a, Kng1, Gbp6, Apoa1, Alb, Fndc5, and Irisin) were shown in volcano plots ( Figures 8A, B and Supplementary Table S3 ). The changes of majority of these myokines were confirmed by qPCR, while Ostn and Gbp6 in the NAFLD and NASH groups showed marginal changes without significance ( Figures 8C, D ). Besides this, we validated the protein levels of MSTN, a classic atrophic marker gene that promotes muscle wasting. Notably, the protein expression of MSTN was induced in the quadriceps muscle of both the NAFLD and NASH groups compared with their controls ( Figures 8E, F ).

Notably, among these genes, the increased levels of Mstn, Gdf11, and Stc2 have been reported to be involved in muscle wasting or muscle dysfunction (37–39). Besides this, Fndn5/Irisin and Apoa1, which have been shown to promote muscle hypertrophy and the systematic improvement of lipid and glucose metabolism (40–42), were largely declined in the muscle of NAFLD and NASH mice, overall suggesting the possible contribution of myokines on the systematic and pathological performance during NAFLD and NASH. In addition, other genes, including Serpina3k and Serpina1a, have not been well understood as myokines and might be involved in muscle functionality and whole-body metabolism during NAFLD and NASH.

In addition, we investigated the uniquely changed myokines in the NASH group and found that Lep, Ly96, Hp, Adipoq, and Serpina3n were uniquely upregulated, while Chad, Cilp2, Wif1, Fetub, Comp, and Mup14 were uniquely downregulated in the NASH group compared with the CON2 groups ( Supplementary Table S3 ), which were potentially affected by NASH progression.