Introduction

Hepatic steatosis (HS) is characterized by the accumulation of triglycerides (TG) in the liver (at least 5% of liver weight) in the absence of secondary contributing factors (Nassir et al., 2015). Metabolic dysfunction-associated fatty liver disease (MAFLD) encompasses a spectrum of conditions, ranging from simple HS to nonalcoholic steatohepatitis (NASH), which is characterized by HS combined with various degrees of liver inflammation and liver injury, with or without fibrosis. NASH can progress to advanced liver disease including cirrhosis and hepatocellular carcinoma (Machado and Cortez-Pinto, 2023). The global prevalence of MAFLD is estimated to be 20∼25% (Younossi et al., 2023), and it is rapidly increasing in parallel with the epidemic of obesity and type 2 diabetes mellitus (Kotsiliti, 2023; Le et al., 2023). Given the absence of approved pharmacological therapy for MAFLD by the US Food and Drug Treatment (FDA) (Llovet et al., 2023; Nassir et al., 2015), the discovery and development of effective therapies are of great importance in MAFLD treatment.

The progression of MAFLD is associated with various factors, including overnutrition, genetic determinants, and associated comorbidities (Sveinbjornsson et al., 2022; Yki-Jarvinen et al., 2021). Strong evidence from meta-analysis, prospective cohort studies, and observational studies indicate that excessive sucrose, a disaccharide that comprises equal amounts of glucose and fructose and is prevalent in all kinds of sugar-sweetened beverages, maybe an initial factor for the development of MAFLD (Geidl-Flueck et al., 2021; Geurtsen et al., 2021; Karlsen et al., 2022; Malik and Hu, 2022; Zhang et al., 2021). For instance, a previous study reported that rats subjected to a 12-week high-sucrose diet (25% w/v sucrose), developed a metabolic syndrome phenotype characterized by glucose intolerance, HS, and insulin resistance (IR) (Sousa et al., 2018). Recent research further reveals that the consumption of sucrose at a concentration of 10% w/v can disrupt energy balance and induce HS within 10-20 weeks (Jang et al., 2020; Stephenson et al., 2022). In this context, c-Jun-N-terminal Kinase (JNK), a mitogen-activated protein kinase family member, emerges as a critical player due to its known function in response to intra- and extracellular stress, including hypercaloric diet consumption (Nikolic et al., 2020; Sabio et al., 2009; Softic et al., 2020). Indeed, accumulating evidence implicated JNK in the pathogenesis of diet-induced obesity, IR, development of metabolic syndrome and MAFLD (Czaja, 2010; Sabio and Davis, 2010; Seki et al., 2012; Singh et al., 2009). JNK has three isoforms: JNK1 and JNK2, which are ubiquitously expressed, and JNK3 is specifically expressed in the brain, testis and heart (Chang and Karin, 2001). Previously, Singh et al. demonstrated that increased hepatic JNK signalling occurred simultaneously with steatosis and lipid accumulation in a steatohepatitis model induced by a methionine- and choline-deficient diet (MCD) (Singh et al., 2009). Notably, MCD diet-fed JNK1 null mice exhibited significantly reduced levels of key hallmarks of MAFLD, including HS, TG accumulation, inflammation, lipid peroxidation, hepatocyte injury, and apoptosis (Singh et al., 2009). Similarly, high-fat diet-fed mice with hepatic ablation of both JNK1 and JNK2 are protected from HS and IR (Vernia et al., 2014). Altogether, these findings suggest that the suppression of the JNK pathway may prevent the development of HS and related metabolic syndrome (Kodama and Brenner, 2009; Schattenberg et al., 2006). However, whether JNK mediates the metabolic disruption induced by sucrose overconsumption remains to be elucidated. Moreover, our previous study demonstrated that JNK-IN-5A, a selective JNK2 and JNK3 inhibitor (Angell et al., 2007; Cerbone et al., 2012), effectively decreased fat accumulation and steatosis-related protein expression in an in vitro steatosis model (Kim et al., 2023). However, the inhibitory effect of JNK-IN-5A to JNK and JNK-related pathways in in vivo is unknown. It is established that JNK exerts its metabolic regulation as a result of multi-tissue communication, including those between liver and adipose tissues (Azzu et al., 2020), and skeletal muscle (Nikolic et al., 2020; Vernia et al., 2014). A systematic investigation into tissue-specific and tissue-crosstalk metabolic changes in response to diet and therapy interventions is, therefore, paramount by using genome-scale metabolic models (GEMs). GEMs are the collection of biochemical reactions and the associated enzymes and transporters between the cellular compartments and such models can be used in the analysis of omics data for gaining insights and interpretation of the omics data (Mardinoglu et al., 2013; Mardinoglu et al., 2018; Mardinoglu and Nielsen, 2012).

In this study, we first performed a 7-day oral toxicology study and studied the tolerability of JNK-IN-5A in rats. Next, we explored the potential therapeutic effect of JNK-IN-5A under high sucrose consumption in in vivo by performing systems biology analysis. To elucidate the systems-wide impact of sucrose overconsumption and JNK-IN-5A treatment, we investigated the intra- and inter-tissue response through integrative analysis of transcriptomics profiles from the metabolically active tissues, including liver and three other extrahepatic tissues, including visceral white adipose tissue (vWAT), skeletal muscle (SkM) and brain tissues. We further predicted the genome-wide metabolic activity changes induced by liquid sucrose ingestion and JNK-IN-5A treatment by performing genome-scale metabolic modelling.

Results

The acute effect of the JNK-IN-5A in rats

We previously showed that the JNK-IN-5A can effectively decrease the TG accumulation and steatosis-related protein levels in an in vitro model simulating the activated de novo lipogenesis in MAFLD patients (Kim et al., 2023). Here, we first performed a 7-day oral toxicology study and studied the tolerability of JNK-IN-5A in rats. Four groups of three males and three females were given vehicle and doses of 30, 100 and 300 mg/kg of JNK-IN-5A, once a day for 7 days (Fig. 1, Fig. 2a). We assessed clinical adverse effects, body weight, and clinical pathological parameters during the study.

Schematic representation of experimental setup.

Study groups included the healthy control group received tap water, the sucrose group received 10% sucrose water, the sucrose+JNK_D1 group received 10% sucrose and 30 mg/kg/day JNK-IN-5A, and the sucrose+JNK_D2 group received 10% sucrose and 60 mg/kg/day JNK-IN-5A (N = 44, n = 11 per group).

Sucrose consumption and JNK-IN-5A treatment exhibit distinct effects on the mRNA expression of genes encoding JNK isoforms in the major metabolic tissues.

(a) The acute effect of JNK-IN-5A in rats (n = 3/Sex/Group). Statistically significant changed clinical chemistry parameters were presented as mean ± standard deviation (SD), see also Data S1. ALB: albumin; ALP: alkaline phosphatase; BUN: blood urea nitrogen; PHOS: phosphate; GLU: glucose; K+: potassium; Na+: sodium; TP: total Protein; GLOB: globulin. Statistical significance was assessed by one-way ANOVA testing followed by Dunn’s multiple comparisons post-test, p-value < 0.05 is considered as statistical significance. (b) Boxplot showing the levels of plasma triglycerides in control, sucrose, sucrose+JNK_D1, and sucrose+JNK_D2 rats. Statistical significance was assessed by one-way ANOVA testing followed by Dunn’s multiple comparisons post-test, p-value < 0.05 is considered as statistical significance. * p < 0.05, ** p < 0.01, *** p < 0.001, **** p < 0.0001. (c) Boxplot showing the relative mRNA expression of Mapk8, Mapk9, and Mapk10 in the studied tissues. The count-based abundance of genes was transformed using the vst function and the adjusted p-value (p.adj) is derived from DESeq2 (Love et al., 2014). * p.adj < 0.05, ** p.adj < 0.01, *** p.adj < 0.001, **** p.adj < 0.0001. (d) Correlation between the relative mRNA expression of Mapk8 and (e) Mapk9 in skeletal muscle tissue (SkM) and plasma triglycerides level.

A few animals including those given vehicles displayed noisy respiration and reflux. These signs were only observed in connection with dosing, suggesting to be related to the high viscosity of the formulation and foul taste, and are not considered as systemic side effects of the JNK-IN-5A administration. We observed none of the administrated doses led to significant changes in body weight and haematological parameters (Data S1). Clinical chemistry parameters revealed lower plasma levels of glucose, blood urea nitrogen and total proteins in female rats given JNK-IN-5A at all three doses (Fig. 2a, Data S1). In addition, the two higher doses of JNK-IN-5A decreased the levels of sodium in female rats and elevated the level of phosphate in both sexes (Fig. 2a, Data S1). Overall, the administration of JNK-IN-5A is well tolerated with regard to adverse clinical signs, body weight and clinical haematology parameters.

Treatment of liquid sucrose-induced fatty liver with JNK-IN-5A

We next investigated the effect of JNK-IN-5A treatment in a rat model of MAFLD induced by liquid sucrose ingestion (Fig. 1b). The rats feeding a standard chow diet received either tap water (control group) or sucrose-containing water (10% w/v sucrose, sucrose group) for 4 weeks, after which rats on sucrose watering were further divided into subgroups: sucrose only group and JNK-IN-5A-treated groups at a regular dosage of 30 mg/kg/day (sucrose+JNK_D1 group) or a higher dosage of 60 mg/kg/day (sucrose+JNK_D2 group), for 3-weeks. A total of 44 rats were used in this study (n = 11 rats per group).

We observed that liquid sucrose consumption resulted in significantly increased plasma levels of TG (p < 0.0001) and lactate dehydrogenase (LDH, p < 0.0001), which is an inflammation marker, compared to the healthy controls (Fig. 2b, Data S2). On sucrose water, JNK-IN-5A-treated rats exhibited significantly lower levels of TG (JNK_D1, p < 0.0001; JNK_D2, p < 0.001), aspartate aminotransferase (AST; JNK_D1, p = 0.09; JNK_D2, p = 0.005), and LDH (JNK_D1, n.s.; JNK_D2, p = 0.015) but had no significant change on body weight, glucose, high-density lipoprotein (HDL) and low-density lipoprotein (LDL) as compared to their corresponding controls (Fig. 2b, Data S2). We also performed neurological tests to study the potential effect of JNK-IN-5A on the brain. We found that there is no significant difference between the groups in terms of learning and memory abilities assessed by the passive avoidance test (Data S3). Furthermore, the locomotor activity test (Data S4) and elevated plus maze test (Data S5) collectively indicate the JNK-IN-5A-treated rats exhibit no evidence of anxiety and behaviour disorders.

To investigate the gene-level inhibition of JNK upon JNK-IN-5A treatment, we further examined the mRNA expression of JNK1 (encoded by the gene Mapk8), JNK2 (encoded by the gene Mapk9) and JNK3 (encoded by the gene Mapk10) in the study. As anticipated, Mapk8 and Mapk9 have high expression in all these tissues (liver, vWAT, SkM, and brain) whereas Mapk10 is highly expressed specifically in the brain tissue (Fig. 2c, Fig. S1a). Our analysis indicates that the treatment of JNK_D1 effectively alleviated the JNK1 gene expression in the liver (padj = 0.03) and exerted more profound inhibitory effects on both JNK1 (padj = 0.0009) and JNK2 gene expression (padj = 2.2e-07) in SkM (Fig. 2c, Data S6). Interestingly, we found that plasma TG level is significantly correlated with the JNK2 gene expression in SkM (Fig. 2d&e). In brain tissue, we found JNK3 gene expression exhibited a trend of dose-dependently inhibition in the JNK-IN-5A-treated groups. Particularly, the downregulation of JNK3 in the sucrose+JNK_D2 group approaches the borderline of significance (padj = 0.07), compared to those, which drank sweetened water group only (sucrose group) (Data S6).

JNK inhibition rewires metabolic perturbation in the liver and vWAT

To gain in-depth insight into system-level transcriptional response to sucrose consumption and JNK-IN-5A treatment, we next studied the global gene expression profiles in both liver and extrahepatic tissues (vWAT, SkM, and brain). A total of 19,780 protein-coding genes in the rat genome were analyzed (Martin et al., 2023). Principal component analysis (PCA) revealed that liquid sucrose ingestion resulted in the largest transcriptional changes in the liver and vWAT (Fig. 3a), which is consistent with previous results (Stephenson et al., 2022). Of note, we observed that JNK-IN-5A treatment, especially JNK_D1, corrected to a large extent of the gene expression profile in the liver and vWAT of treated-rats towards the profile of healthy controls on a principal components space (Fig. 3a). In SkM, the transcriptional difference is primally driven by JNK_D1 treatment and showed a clear separation (Fig. 3a).

System-wide transcriptomics profiling revealed tissue-specific metabolic rewiring by sucrose consumption and JNK-IN-5A treatment

(a) Principal component analysis (PCA) of liver, vWAT, SkM, and brain transcriptome data showing global gene expression profiles in control (n = 9), sucrose (n = 11), sucrose+JNK_D1 and sucrose+JNK_D2 (n = 10/per group). Each data point represents a sample in the respective coloured group. (b) Bar graphs represent the percentage of differential expressed genes (or regulated genes) associated with sucrose (sucrose vs. control), sucrose+JNK_D1 (sucrose+JNK_D1 vs. sucrose), or sucrose+JNK_D2 (sucrose+JNK_D2 vs. sucrose) in all protein-coding genes (N = 19,780) for each tissue. (c) Heatmap of all genes differentially upregulated (red) and downregulated (blue) in three pairwise comparisons for sucrose (sucrose vs. control), sucrose+JNK_D1 (sucrose+JNK_D1 vs. sucrose), and sucrose+JNK_D2 (sucrose+JNK_D2 vs. sucrose) (adjusted p-value, p.adj < 0.01) for each tissue. Column annotation represents tissues. (d) Venn diagram showing the overlap between significantly regulated genes response to sucrose consumption and JNK_D1 treatment, see also Fig. S1.

In parallel with the results presented above, differential gene expression analysis revealed that the liver had the largest number of differentially expressed genes (DEGs, DESeq2-Benjamini-Hochberg adjusted p < 0.01) associated with liquid sucrose ingestion, comprising a total of 2,771 DEGs (14%), followed closely by vWAT with 2,724 DEGs (13.8%), SkM with 257 DEGs (1.3%), and the brain with a modest 30 DEGs (0.2%) (Fig. 3b&c, Fig. S1b-e, Data S6). Moreover, we observed that JNK_D1 treatment demonstrated a remarkable effect in the reversal of gene expression patterns related to sucrose consumption. Specifically, in the liver, JNK_D1 treatment substantially suppressed the expression of 60.7% (1,143 out of 1,884, hypergeometric p ≍ 0) sucrose-induced DEGs and elevated the expression of 22% of sucrose-repressed genes (195 out of 887, hypergeometric p = 6.9e-129) (Fig. 3d, Fig. S1b). Likewise, we observed a reversal of 42.2% and 31.4% (698 of 1,655 and 336 out of 1,069, respectively) DEGs in vWAT (Fig. 3d, Fig. S1c). It should be noted that, in both liver and vWAT, the effect of JNK_D2 treatment is relatively smaller than JNK_D1 in terms of the corrections of DEGs (Fig. S2a&b). Pathway enrichment analysis of JNK_D1 reversed genes identified significantly enriched KEGG pathways related to carbon metabolism, fatty acid metabolism, protein processing in ER, and oxidative phosphorylation (OXPHOS) (Benjamini-Hochberg adjusted p < 0.01) (Fig. S2c&d). Interestingly, some of the genes exhibited opposing regulation in the liver and vWAT in response to the interventions of sucrose only and the combination of sucrose intervention and JNK-IN-5A treatment. This is exemplified in genes related to fatty acid degradation, elongation and biosynthesis (up in the liver and down in vWAT in sucrose only group, while down in the liver and up in vWAT in sucrose plus JNK-IN-5A-treated groups), including Cpt2, Acads, Acadvl, Acat1, Acot5, Eci3, Hsd17b12, and Oxsm (Fig. S3, Data S6); genes related to fatty acid (FA) uptake, including fatty acid binding proteins Fabp1, Fabp4 and transporters Slc27a1 and Slc27a5 (Fig. S4, Data S6); and genes related to redox homeostasis, including Prdx4, a secreted enzyme that scavenges redox oxidative stress (ROS), Selenos, and Ndufa6 (Fig. S5, Data S6). Taken together, these results suggest that JNK-IN-5A treatment reshapes sucrose-induced transcriptional changes in both the liver and vWAT, reinforcing their metabolic synergy in FA uptake and metabolism (Stephenson et al., 2022).

JNK inhibition in SkM correlated with the regulation of energy metabolism in the liver and vWAT

It has been reported that JNK acts as a critical mediator of muscle phenotype and adaptive remodelling in animals and humans (Lessard et al., 2018). Here, we found that in contrast to the massive transcriptional regulation observed in the liver and vWAT, liquid sucrose intake has minimal impact on SkM in terms of both genome-wide transcriptional variance (Fig. 3a) and the number of DEGs (Fig. 3a, Fig. 4a&b). While in JNK_D1-treated SkM, a substantial quantitative difference in gene expression emerged, encompassing 4,159 (20.9%) DEGs, of which 97% up-regulated (1,552 out of 1,588) and 99% down-regulated (2,557 out of 2,571) genes being exclusively dependent on the JNK_D1 treatment (Fig. 4a&b). Functional analysis revealed that many of the upregulated genes are involved in immune-related pathways, including hematopoietic cell lineage (HSCs), osteoclast differentiation, chemokine signalling pathway, phagosome, and Fc gamma R-mediated phagocytosis, and down-regulated genes are involved in OXPHOS, thermogenesis, TCA cycle, carbon metabolism, cardiac muscle contraction, and insulin signalling pathway (Fig. 4c&d, Fig. 5). Among the genes regulated by both sucrose-feeding and JNK_D1 treatment in SkM were those involved in the insulin signalling pathway (e.g., Ppp1r3f, Phkg1, Phkb) (Fig. 5). As this intra-tissue effect was observed in JNK_D1 not in the JNK_D2-treated group (Fig. 3b, Fig. S5a), we subsequently investigated the inter-tissue regulation based on the significantly transcriptional changes (as indicated by the log2-Fold change, log2FC) in response to JNK_D1 treatment and found that log2FC of shared DEGs between JNK_D1-treated SkM and liver, vWAT or brain exhibited strong agreements (Fig. 4e, Fig S5b). Specifically, 86% and 69% of shared DEGs between SkM and vWAT, and the brain, exhibited consistent changes in response to JNK_D1 treatment (the same log2FC sign, indicating either upregulated or downregulated in both tissues being compared), with a Spearman correlation of 0.78 and 0.54, respectively. These shared DEGs were found to be involved in several critical metabolic processes, including cholesterol metabolism, bile acid biosynthesis and secretion (such as Lipc, Apob, Slc7a1, Angptl3, Slc10a1), and regulation of energy homeostasis (e.g., Lep), suggesting that SkM had a major metabolic adaptation to JNK_D1 treatment and exerts its metabolic effect, at least partially, through coordinated interactions with other tissues.

JNK inhibition in SkM correlated with the regulation of energy metabolism in the liver and adipose tissues.

(a-b) Venn diagram showing the overlap between significantly regulated genes response to sucrose consumption and INK_D1 treatment in SkM, see also Fig. S1&Fig. S5. (c-d) Functional annotation of up-(red) and down-regulated (blue) genes in JNK_D1-treated SkM with a Benjamini Hochberg adjusted p-value < 0.05. (e) Correlation of log2(Fold change) for shared DEGs between the SkM and Liver, vWAT, or Brain after JNK_D1 treatment.

JNK inhibition regulates insulin signalling-related genes in SkM

(a) A representative diagram of metabolic pathways associated with insulin resistance and insulin signalling was found to be largely inhibited by JNK_D1 treatment. (b) Correlation between the relative mRNA expression of Phkg1, Phkb, and Ppp1r3f (those three genes are also significantly upregulated in sucrose-feeding only group) in SkM and plasma triglycerides level. (c) Relative changes of genes involved in insulin resistance and insulin signalling pathways in response to sucrose consumption and JNK-IN-5A treatment in the tissues. The adjusted p-value (p.adj) is derived from DESeq2 (Love et al., 2014). * p.adj < 0.05, ** p.adj < 0.01, *** p.adj < 0.001, **** p.adj < 0.0001.

Inter-tissue network analysis revealed the JNK-IN-5A inhibition-associated molecular mechanism

To better understand the regulation of lipid metabolism and energy homeostasis at a system-wide level in an unbiased manner, we next performed gene co-expression network (GCN) analysis to study intra- and inter-tissue coordination of metabolic pathways as well as biological processes that are associated with sucrose-feeding and/or JNK_D1 inhibition. GCN, a robust systems biology approach (Barabasi and Albert, 1999), allows us to study the functionality of the gene module and its constituent genes. These genes are more likely to be co-expressed if they are either regulated by the same transcription factors or share similar topological properties in protein-protein interaction networks (Califano et al., 2012; Dobrin et al., 2009; Zhang et al., 2016). By investigating the importance and topology of these functional modules, it is possible to understand the biological mechanisms underpinning the disease and/or therapy intervention (Choobdar et al., 2019; Yang et al., 2021; Yang et al., 2014).

We first constructed tissue-specific GCN based on the whole transcriptomics data and selected highly connected genes from the network (the top 10% positively correlated genes that fulfilled FDR < 0.05), followed by module detection using the Leiden clustering algorithm (Fig. S7a; Data S7). We subsequently evaluated the enrichment of DEGs in the modules to identify those associated with sucrose consumption, JNK_D1 treatment, or both (so-called perturbated module). As expected, GCN analysis revealed distinct gene expression patterns across different tissues (Fig. 6a). For instance, in module 1 in the liver (Liver.M1), vWAT.M4, SkM.M4, and Brain.M4, we found a significant enrichment of genes that exhibited elevated expression in response to sucrose consumption while simultaneously being suppressed by JNK_D1. Conversely, Liver.M0 and vWAT.M3 exhibited significant enrichment of genes repressed by sucrose intake but elevated after JNK_D1 treatment. These distinct expression patterns collectively reflect the combined metabolic effect of sucrose-feeding and JNK_D1 interventions. Additionally, GCN analysis identified modules associated with either sucrose-feeding or JNK_D1 intervention at the tissue level, for example, vWAT.M0 (sucrose), vWAT.M1 and Brain.M5 (JNK_D1) (Fig. 6a).

Inter-tissue network analysis identifies JNK-IN-5A inhibition-associated molecular mechanism.

(a) Significantly enriched modules (hypergeometric test, p < 0.05) by the differentially expressed genes related to sucrose, sucrose+JNK_D1 and sccrose+JNK_D1 in each tissue. (b) Dot-plot showing the significant correlation among inter-tissue module pairs. The size and colour of the connected line are proportional to the correlation coefficient and statistical significance indicated by the adjusted p-value. The module pairs with Benjamini Hochberg adjusted p-value (p.adj) < 0.05 are presented. (c) Spearman coefficient correlation of the first component from PCA analyses on gene expression of module pairs. (d) Significantly enriched MSigDB hallmark gene sets by each module.

To pinpoint gene modules involved in inter-tissue communication, we conducted a correlation analysis among inter-tissue modules using their module eigengene (see details in Method) and assessed the enrichment of genes from these modules in MSigDB hallmark gene sets. Our analysis revealed a significant correlation between SkM.M0 and Liver.M0 (R = 0.67, padj = 0.006) as well as vWAT.M3 (R = 0.58, padj = 0.03) (Fig. 6b&c). This correlation aligns with the functional annotation of genes within these modules, including those associated with mTORC1 signalling, PI3K/AKT/mTOR signalling, DNA Repair, and unfolded protein response (Fig. 6d). Interestingly, we observed that Mapk9 is among the module members of SkM.M0, indicating a potential metabolic link coordinated by JNK inhibition, connecting SkM with other tissues, including the liver and vWAT (Fig. 6b, Data S7). SkM.M1 and Liver.M2 (R = −0.66, padj = 0.006), SkM.M3 and vWAT.M0 (R = 0.60, padj = 0.03), and SkM.M4 and vWAT.M3 (R = −0.66, padj = 0.006) were also found to be significantly correlated (R = -0.57, padj = 0.03). Note that genes in Liver.M1, vWAT.M3, SkM.M0, SkM.M4, and Brain.M4 are commonly involved in adipogenesis, fatty acid metabolism, OXPHOS, and reactive oxygen species pathway (Fig. 6d). These observations suggest the potential crosstalk or regulatory relationship between these tissue modules after JNK_D1 treatment.

Metabolic modelling identified the metabolic reprogramming linked to JNK-IN-5A inhibition

To predict how tissue-level metabolic activity responds to sucrose feeding and JNK-IN-5A treatment, we performed an integrative analysis of transcriptomics data by generating tissue-specific GEMs and applying Compass (Wagner et al., 2021). Compass is a flux balance analysis-based (Orth et al., 2010) algorithms to model the metabolic state of a system (e.g., single cell or tissue), taking into account the observed expression levels of enzyme-coding genes in each sample from RNA-seq data. In this study, we employed Recon3D, a generic GEMs of human metabolism, comprising 10,600 metabolic reactions associated with 5,835 metabolites and 2,248 metabolic genes (Brunk et al., 2018). Consistent with our previous analysis, we observed a shift in metabolic activity in the liver and vWAT after JNK_D1 treatment (Fig. 7, Data S8). Specifically, Compass predicted that 2,840 and 4,170 metabolic reactions showed significantly altered activity (Wilcoxon test, adjusted p-value < 0.05) following the sucrose consumption in the liver and vWAT, respectively. 89.2% and 88.3% of which were correspondingly reversed after JNK_D1 treatment (Data S8). The common altered reactions are partitions into the subsystems named citric acid cycle, fatty acid oxidation, fatty acid synthesis, and tryptophan metabolism in Recon3D (Fig. 7c). Interestingly, we found 8 metabolic reactions significantly altered in sucrose+JNK_D2-treated rats, including 2 reactions involved in glycolysis/gluconeogenesis, 2 in nucleotide interconversion, 1 in pentose phosphate pathway, 1 in fructose and mannose metabolism, 1 in propanoate metabolism (Cohen’s d < 0, padj < 0.05), and 1 in drug metabolism (Cohen’s d > 0, padj < 0.05) (Data S8). Hence our systems biology analysis revealed the molecular mechanisms altered due the JNK-IN-5A treatment.

Metabolic modelling highlights metabolic reprogramming linked to JNK-IN-5A inhibition

(a) The numbers of differentially altered metabolic reactions in response to sucrose, sucrose+JNK_D1, and sucrose+JNK_D2, respectively, see also Data S13. (b) Distribution of Cohen’s d statistic for differential metabolic reactions in response to sucrose, sucrose+JNK_D1 and sucrose+JNK_D2 in each tissue, see also Data S13. (c) Principal component analysis of the Compass score of metabolic reaction potential activity scores showing the metabolic heterogeneity in response to sucrose ingestion and JNK-IN-5A treatment at the tissue levels, for principal components 1 and 2. (d) Heatmap showing the tissue-level differential activity of metabolic reactions in response to sucrose intake and JNK-IN-5A treatment. Reactions are partitioned by Recon3 pathways (Brunk et al., 2018) and coloured by the sign of their Cohen’s d statistic derived from different contrasts: sucrose vs control, sucrose+JNK_D1 vs sucrose, or sucrose+JNK_D2 vs sucrose, respectively. Abbreviations: vWAT, visceral white adipose tissue; SkM, skeletal muscle.

Discussion

The aetiology of MAFLD is multifactorial and remains incompletely understood, impeding the progress in developing effective therapies for MAFLD. Our investigation into the metabolic changes from the major metabolically active tissues enabled us to systematically assess the potential therapeutic effects of JNK-IN-5A under high sucrose consumption. Here, we demonstrated that JNK-IN-5A attenuated the accumulation of circulating TG and inflammation exaggerated by liquid sucrose consumption. Importantly, the rats treated with JNK-IN-5A did not exhibit indications of anxiety or behavioural disorders. These results have important implications for the development of JNK3 inhibitors, taking into account the high expression of JNK3 in the brain (Chang and Karin, 2001). Additionally, our data indicated that JNK-IN-5A, a selective inhibitor of JNK2 and JNK3 (Angell et al., 2007; Cerbone et al., 2012), also exerted an inhibitory effect on the gene expression of JNK1 in in vivo, especially in skeletal muscle (Fig. 2).

By studying the hepatic and extrahepatic tissues using genome-wide transcriptome data and integrative analyses with GEMs, we found that the intervention with sucrose-sweetened water triggered extensive transcriptional changes in the liver and adipose tissues, including those involved in fatty acid and oxidative metabolism. The metabolic axis between the liver and adipose tissues is a critical regulatory network that governs energy balance and metabolic homeostasis in the body (Azzu et al., 2020). It has been reported that upon chronic overnutrition, adipose lipogenesis is significantly downregulated, which leads to ectopic lipid accumulation and insulin resistance (Jeon et al., 2023). In line with this, we found that in the sucrose-sweetened water group, the mRNA expression of many genes related to fatty acid metabolism and uptake is downregulated in adipose tissue. Interestingly, those genes showed the opposite regulation between liver and adipose tissue, indicating their complementary role in these metabolic processes in response to the diet intervention. These findings also pinpoint the critical role of extrahepatic metabolism in the development of MAFLD.

JNK plays essential roles in a range of signalling pathways that have implications for metabolic processes such as insulin sensitivity, thermogenesis, adipogenesis, and lipid metabolism (Nikolic et al., 2020; Solinas and Becattini, 2017). Our data demonstrated that JNK-IN-5A, especially when administered at a regular dosage (JNK_D1, 30mg/kg/day), effectively opposite a substantial portion of the transcriptional changes induced by sucrose consumption, particularly those related to carbon, fatty acid metabolism, OXPHOS, and thermogenesis in the liver and adipose tissues. This suggests that JNK-IN-5A might hold promise as a potential strategy for restoring metabolic perturbation in tissues challenged by sucrose-induced stress. However, our study also revealed that JNK_D1-treated skeletal muscle exhibited profound transcriptional changes related to HSCs, osteoclast differentiation, and immune-related pathways. Importantly, we found that these changes were not isolated to SkM but exhibited correlations with changes in liver, vWAT, and even brain tissue, evidenced by the correlation analysis and co-expression network analysis across tissues. Previously, Xiao et al. (Xiao et al., 2019) demonstrated that the JNK pathway regulates HSCs expansion and self-renewal in JNK-IN-8-treated CD34+CD45RA-CB cells. Lessard et al. (Lessard et al., 2018) demonstrated that mice with muscle-specific JNK1/2 knockout, have a muscle phenotype that mimics that of endurance-trained athletes. Together with our data, these results suggest a multifaceted role of JNK signalling in the regulation of HSCs and metabolism in skeletal muscle. Further investigation is warranted to understand the precise mechanisms underlying these changes and their implications for skeletal muscle function and overall physiological homeostasis.

In summary, our study provides a comprehensive understanding of the therapeutic potential of JNK-IN-5A treatment in a rat model of sucrose-induced MAFLD. Through an analysis of genome-wide transcriptome data in the major metabolic tissues, we shed light on the JNK-IN-5A-mediated metabolic effects in both hepatic and extrahepatic tissues. This implies a potential therapeutic strategy for the treatment of diet-induced MAFLD and other metabolic complications, potentially by regulating cellular energy homeostasis and mitochondrial metabolism. Additional studies are needed to determine whether JNK inhibition can exert its therapeutic benefits in long-term dietary interventions, resembling the different stages of MAFLD, including the NASH model.

Materials and Methods

Experimental models and subject details

1. Oral Toxicity Study in Rats

Twelve male and twelve female rats of the Wistar strain (Envigo, Venray, Netherlands) were used to assess the 7-day oral (gavage) toxicology, in terms of clinical signs, body weight and clinical pathological parameters, of JNK-IN-5A in rats. Animals were group-housed under standard conditions and had free access to water and 2916C Teklad global 16 % protein rodent irradiated diet (Envigo, Venray, Netherlands). Dose levels of 0 (vehicle control), 30, 100, or 300 mg/kg of JNK-IN-5A were administered once daily via oral gavage for 7 days. The vehicle was 2 w/w% HPMC (Colorcon), 2 w/w% PS80 (Merck) in 10 mM PBS, pH 7 and dose volume was 5 mL/kg for JNK-IN-5A. Body weights were recorded the week before and within 24 hours before the start of dosing and on Day 3 and on necropsy Day 8. All animals were bled at 2 hours after dosing on Day 1 and Day 7. Blood samples were placed on ice after sample collection and centrifuged at approx. 3000xG, 5 min., 4°C. 7.5 μL plasma was transferred into a separate low-binding protein tube and the remaining plasma was transferred into a second low-binding tube. The tubes were frozen upright at approximately - 20°C. Blood and plasma were analysed using an Exigo analyser and Vetscan equipment. The VetScan analyser was calibrated with calibration samples with known plasma levels of each parameter before the start of analysis, occasionally during the analysis period of study samples and after completed analyses. All the animal procedures (including housing, health monitoring, dosing, etc.) and ethical revisions were performed according to the 2010/63/EU Directive on the protection of animals used for biomedical research.

2. JNK-IN-5A Treatment Study in Rats

All experiments were performed in accordance with institutional ethical guidelines, and the present in vivo study was approved by the Ethics Committee for Animal Research of Atatürk University (with the reference number 2023/04 on 30.03.2023). Rats were obtained from the Atatürk University Experimental Research Center. All groups were housed on a 12-hour light and dark cycle. The temperature and humidity values of the environment were kept at optimal values. To prevent neophobia, rats were divided into groups one week before sucrose treatment by following the double-blind method to adapt to the laboratory environment and each other.

A total of 4 groups consisting of 11 subjects each were formed, including the healthy Control group consisting of rats received ad libitum tap water, the Sucrose group consisting of rats received ad libitum 10% Sucrose (BioShop, Canada, Lot No: 2E77269) (equivalent to the concentration of fructose in typical soda (Jang et al., 2020)), and groups treated by JNK-IN-5A at two concentration: sucrose+JNK_D1 group with rats received ad libitum 10% sucrose and 30 mg/kg/day JNK-IN-5A gavage and sucrose+JNK_D2 group with rats received ad libitum 10% sucrose and 60 mg/kg/day JNK-IN-5A gavage. The rats had ad libitum access to tap water (Control) or 10% sucrose water (sucrose, sucrose+JNK_D1, and sucrose+JNK_D2) and a standard chow diet throughout the experimental period. The amount of water consumed by all groups was recorded daily, and drinking water was renewed every day. After feeding with sucrose for 4 weeks, drug treatment was given at the following 3 weeks. Behavioural tests were taken before the experiment started and after 3 weeks of drug treatment. Locomotor activity tests, elevated plus maze tests and passive avoidance tests were used as behavioural tests. Behavioural tests were repeated at the end of the 21-day treatment regimen. While the last behavioural tests were being taken, drug treatment continued so that the drug level in the blood would not change. After the experimental animals were sacrificed with a high dose (200 mg/kg thiopental sodium) anaesthetic, the relevant tissues were taken for further analysis.

2.1 Monitoring behaviour

All tests were performed twice in total, initially and after the 21-day treatment regimen. According to the results obtained from the initial measurements, rats with outliers for the experiment were determined. After the treatment measurement, it was determined whether there was any change in the behaviour of the rats. The Locomotor Activity Test is a system used to test differences in motor activity and anxiety-related behaviours of rats. In the system scanned with infrared rays, each of the ambulatory, stereotypical, rearings, resting percentage and travelled distance movements of the rats are digitized and converted to computer output. These outputs are evaluated and provide information that needs to be interpreted about the animal’s adaptation to the new environment and the continuation of the sense of curiosity or anxiety. Stereotypical movement shows rats grooming and recognizing their environment. Ambulatory movement is all the movements that rats make on the ground without standing up. Distance shows the distance travelled by the rats in the cage. The rearing movement indicates the time the rats stand up in the cage. The resting percentage is an indicator of anxiety. Ambulatory movements are recorded with horizontal sensors that surround the cage 360 degrees at the base, and stereotypic and rearing movements are recorded with vertical sensors located 10 cm above the floor. The entire cage continues to be measured continuously with 250 ms infrared. Unlike the open field test, activity measurement is performed in this study and the computer automatically records all the recordings instead of an independent observer (Lu et al., 2019). The experiment is performed in a dark environment. Rats are placed in cages for 10-minute periods and the system automatically records with infrared rays. After each rat, the cage is cleaned with alcohol to remove odor and urine (Ferah Okkay et al., 2022). The Elevated Plus Maze Test is used to measure anxiety-like behaviour in rodents.

The elevated plus maze apparatus is a system consisting of two open arms (without walls; 50 cm × 10 cm) and two closed arms (with black walls, 50 cm × 10 cm). The open and closed arms are connected by a central platform. The labyrinth is 55 cm above the ground. At the start of each test, rats are placed in the central zone and allowed to explore the platform for 5 minutes. The rats are recorded with a camera located at the top of the maze and the test is done in a quiet, dim room. The number of times the rats entered the open and closed arms and the time spent in these arms were recorded. The device was cleaned with an alcohol solution between tests (Chellian et al., 2020; Walf and Frye, 2007). Results are given as arm entries and time spent (durations) in these areas. The more entries into open arms and the more duration spent there, the less anxiety the rat has. The Passive Avoidance Test is a system consisting of 2 rooms (light/dark), the base of which is covered with rods that can apply electroshock, in which memory, behaviour and learning abilities are investigated. This system is used to test the learning ability of the animal by fear-driven conditional avoidance. In this test, the animal is expected to learn to avoid entering the electroshock chamber (dark room). Rats are nocturnal animals. By nature, they sleep during the day and are active at night. Although rats are not afraid of light, they generally live in dark compartments because they are uncomfortable with the bright environment. In the passive avoidance system, lighting the bright compartment with LED lights and high lumen light creates an uncomfortable environment for the animals, and they want to move to the closed compartment. Here, the low-intensity electroshock they are exposed to causes a deterrent stimulus and allows the rats to understand and learn how the system works. The experimental procedure is built on these understanding and learning abilities (Gimenez De Bejar et al., 2017; Hosseini et al., 2013; Wu et al., 2020). The experiment is carried out on two consecutive days. Learning is measured on the first day and recall on the second day. The time that each rat took to enter the dark compartment was considered as initial latency. The next day (retention phase), the rats were reintroduced into the bright chamber and the step-through latency was noted down as memory retention. Shorter latencies indicated poorer cognition. At this stage, rats are expected to remember the shock in the dark compartment. For this reason, it was accepted that the memory of the rats that stayed longer in the bright compartment was stronger. In retention measurements, it is understood that the longer the rats stay in the bright compartment, the stronger the memory and learning.

2.2 RNA extraction and library preparation for transcriptome sequencing

Total RNA is extracted from the snap-frozen rat tissues with Zymo Research Quick RNA Miniprep kit (California USA) by following the protocol provided on the manufacturer’s website with slight modifications made to improve the quality and quantity of RNA from each tissue. Before the cell lysis, the tissues are subjected to short (1 min or longer for muscle tissues) homogenization with the beads in PowerLyzer 24 by keeping samples on ice to prevent overheating. The homogenized tissues are further enzymatically in a digestion buffer for at least 30 minutes or more for muscle tissue at room temperature with continuous inversion and agitation of the tubes. After the lysis of the cells, the RNA is cleaned and collected with columns. The quality and quantity of the isolated samples were examined spectrophotometrically with NanoDrop (Thermofisher, USA). After the evaluation of the sample A260, A260/A280, and A260/A230 values, the RNA integrity number (RIN) value of total RNA samples was measured with TapeStation (Agilent Tech, USA). The amount of total RNA was measured more precisely fluorometrically with the RNA Broad Range kit (Thermofisher, USA) by using Qubit (Invitrogen, USA). RNA sequencing library was prepared with Illumina Stranded Total RNA Prep and Ligation with Ribo-Zero Plus kit was used by following the standard protocol provided by the manufacturer. The libraries were then pair-end (2×100bp) sequenced on the NovaSeq6000 system yielding, on average, 25 million fragment reads per sample. Raw sequencing data (.bcl) was converted to FASTQ with the Dragen Bio-IT platform (v3.9.5). The quality of RNA-seq data was assessed by FastQC (v0.11.9).

2.3 RNA-seq data processing, differential expression, and functional enrichment analysis

Tissue bulk RNA-seq data were aligned and quantified using a standard protocol of Kallisto (v0.46.2) (Bray et al., 2016) against the Rattus norvegicus genome (mRatBN7.2) downloaded from Ensembl (Martin et al., 2023) official website (https://www.ensembl.org/index.html). The output of Kallisto, both estimated counts and TPM (transcript per kilobase million)-based transcript-level expressions were then transformed into gene-level expressions using the Bioconductor package tximport (v1.22.0) with the tx2gene option set to connect transcripts to genes (Soneson et al., 2015). Protein-coding genes were considered for the above step and downstream analyses. Differential analysis was performed using the DESeq2 Bioconductor package (v1.34.0) (Love et al., 2014), following a standard protocol for all three pairwise intervention comparisons for determining the effect of sucrose consumption only (sucrose vs control), the effect of 30 mg/kg/day JNK-IN-5A treatment (sucrose+JNK_D1 vs sucrose), and effect of 60 mg/kg/day JNK-IN-5A treatment (sucrose+JNK_D2 vs sucrose) in all four studied tissues. Significantly expressed genes (DEGs) were identified with a significance threshold of an adjusted p-value < 0.01. To examine differences in global gene expression, we conducted principal component analysis (PCA) based on variance-stabilizing transformation (VST)-normalized count data in DESeq2. PC1 and PC2 were plotted for visualization.

Gene Co-expression Network Analysis

For the generation of the tissue-specific co-expression networks, we first filtered out lowly-expressed genes based on their TPM-based expression level (TPM <1) and constructed co-expression networks by generating gene pairs Spearman correlation ranks within the tissue, which was performed using the “spearmanr” function from SciPy (Virtanen et al., 2020) in Python 3.8. Next, considering the network with negative correlation has relatively low correlation scores, we retained the top 10% positively correlated genes that fulfilled FDR < 0.05 on the network (Arif et al., 2021; Yang et al., 2021) and performed module detection analysis using Leiden algorithm (Traag et al., 2019), implemented in the leidenalg (v0.7.0) package with “ModularityVertexPartition’’ to find the optimal partition. Modules with less than 30 genes were discarded to be able to get significant functional analysis results in the downstream analysis. To identify gene modules involved in tissue crosstalk, we correlated the module eigengenes of inter-tissue modules, and assessed the significance of these correlations after Benjamini & Hochberg (BH) correction, with BH-adjusted p < 0.05 was considered as significance. The module eigengene is defined as the first principal component of the expression matrix of a given module.

Gene Set Functional Enrichment Analysis

We used Bioconductor package enrichR (v3.2) to annotate gene sets (e.g., differentially expressed genes or co-expression genes) against the curated genes sets obtained from the KEGG pathway database (Kanehisa, 2002), or against hallmark gene sets from MSigDB (Liberzon et al., 2015) using R package msigdbr (v.7.5.1). Benjamini-Hochberg (BH)–adjusted P values (p.adj) < 0.05 are considered as statistically significant and are provided in the relevant figures/datasets.

2.4 Compass-based Metabolic Activity Analysis

We inferred tissue-level metabolic state from RNA-seq profiles using Compass (Wagner et al., 2021). Compass is a flux balance analysis-based algorithm (Orth et al., 2010) that enables modeling the metabolic state of a system (e.g., single-cell or organ) from RNA-seq data based on a reconstructed genome-scale metabolic model of human metabolism (also known as Recon3D) (Brunk et al., 2018). We downloaded Recon3D from BiGG Models (http://bigg.ucsd.edu/models/Recon3D/) and created model metadata files required for the input of downstream analysis, including gene, reaction, and metabolite metadata, using in-house scripts. We used human orthologs of rat gene expression data (TPM-based) for computation of the Compass scores matrix using default parameters and the human species as detailed in (Wagner et al., 2021). Downstream analysis was done using the compassR R package (v.1.0.0) following the protocol (https://github.com/YosefLab/compassR).

Statistical Analyses

One-way analysis of variance (ANOVA) test was used to evaluate the significance between the control and treatment groups in behaviour tests using SPSS software (version 20). For the statistical analysis of biometric and plasma parameters, data were shown as mean ± standard deviation (SD), unless stated otherwise. The assumption of normality was determined using the Shapiro–Wilk test. Differences between multiple groups were tested by ANOVA or a non-parametric a non-parametric Kruskal-Wallis test if the data were not normally distributed, followed by Dunn’s multiple comparisons post-test. All results were considered statistically significant at p < 0.05, unless stated otherwise. A hypergeometric test was used for testing the overlapping between two gene sets. The Jaccard index was used to evaluate the similarity and diversity of two gene sets.

Acknowledgements

The authors would like to acknowledge financial support from ScandiEdge Therapeutics and the Knut and Alice Wallenberg Foundation (No. 72110). A.M. and H.Y. acknowledge support from the PoLiMeR Innovative Training Network (Marie Skłodowska-Curie Grant Agreement No. 812616) which has received funding from the European Union’s Horizon 2020 research and innovation program. The computations were performed on resources provided by SNIC through the Uppsala Multidisciplinary Center for Advanced Computational Science (UPPMAX) under Project sllstore2017024 and sctatlas.

Author contributions

A.M. conceived and designed the analysis; M.K., C.B., I.B., O.O.T., C.B., N.Y., S.Y., S.I., J.S., A.H., H.T. collected data for the study; H.Y., C.Z., W.K., M.N.SH. performed the data analysis, visualization, and interpretation; H.Y. drat the manuscript; All authors read and approved the final manuscript.

Declaration of interests

The authors declare no conflict of interest.

(a) Gene expression pattern of JNK isoforms (b-e) UpSet plot showing the number of shared differentially expressed genes (DEGs) (adjusted-p < 0.01) associated with sucrose consumption and JNK-IN-5A treatment in the liver, visceral white adipose tissue (vWAT), skeletal muscle tissue (SkM) and brain tissues in the experimental rats. The histogram bar graph (right) shows the total number of DEGs in each comparison corresponding. The histogram bar graph (top) shows the number of DEGs shared among pairwise comparisons. Connected lines indicate what comparisons are sharing DEGs.

(a) Venn diagram showing the overlap between significantly regulated genes response to sucrose consumption and JNK_D2 treatment in liver and (b) in vWAT, related to Fig. 3. (c) Top-10 enriched KEGG pathways enriched by reversed genes after JNK_D1 treatment in liver and (d) in vWAT.

Relative changes of genes involved in fatty acid degradation, elongation, biosynthesis, and metabolism in response to sucrose consumption and JNK-IN-5A treatment in the tissues, related to Fig. 2.

Relative changes of genes involved in fatty acid uptake in response to sucrose consumption and JNK-in-5A treatment in the tissues, related to Fig. 2.

Relative changes of genes involved in redox homeostasis in response to sucrose consumption and JNK-in-5A treatment in the tissues, related to Fig. 2.

(a) Venn diagram showing the overlap between significantly regulated genes response to sucrose consumption and JNK_D2 treatment in SkM, see also Fig. S1&Fig. S5. (b) Correlation of log2(Fold change) for shared DEGs between the SkM and Liver, vWAT, or Brain after JNK_D2 treatment. (c) Venn diagram showing the overlap between significantly regulated genes response to sucrose consumption and JNK_D1 treatment as well as (d) JNK_D2 treatment in the Brain, respectively.

(a) The number of gene modules are identified from the tissue-specific co-expression networks. The text on the module indicates the module name and the number of genes of individual modules in each tissue. The edges between the clusters were aggregations of the inter-cluster edges. (b) PCA of Liver.M0, vWAT.M3, and SkM.M0 gene expression data. (c) Module neighbors of Mapk9 in SkM.M0.