Main

NAFLD is the most common chronic liver disease. Its increasing prevalence is driven by high-calorie diets and sedentary lifestyles. Central obesity, insulin resistance and type 2 diabetes (T2D) are strong independent risk factors of NAFLD1,2. NAFLD exists on a histological continuum encompassing stages ranging from isolated steatosis (NAFLD) to NASH, characterized by lobular inflammatory infiltrates, hepatocyte ballooning and cell death, to fibrosis and ultimately cirrhosis3. NASH development results from complex interactions of metabolic and stress pathways in hepatocytes, initiated by chronic excessive lipid accumulation, with inflammatory processes driven by various immune cell populations, collectively inducing the histological picture of an active steatohepatitis1. Several lesions can be present, but lobular inflammation and ballooning are the most relevant histological markers of NASH and their combination is referred to as NASH activity, clearly distinguishing the activity of the steatohepatitis from the features of steatosis4,5.

Ballooned hepatocytes are thought to be stressed and damaged cells that lose their rectangular shape and swell due to cytoskeleton degeneration, possibly responding inadequately to pro-apoptotic and danger signals6.

Inflammatory infiltrates within the liver lobules are a hallmark of active steatohepatitis, and specialized immune populations, both resident and infiltrating, are linked with NASH7. Although certain circulating and hepatic immune populations have been associated with NASH8,9,10, a systematic and in-depth analysis of the cellular immune system in NASH is missing. Despite current strategies to treat NASH targeting to reduce both lobular inflammation and ballooning11,12, the molecular mechanisms underlying these components of NASH are poorly understood. Previous studies mainly aimed to identify molecular pathways correlated with NASH versus no NASH12,13, without distinguishing steatosis from disease activity. Functional signatures associated with NASH and the activity that distinguishes NASH from steatosis as such remain unexplored. Transcriptional signatures of NASH and its activity can be identified by comparing liver transcriptomics from patients with histologically proven lobular inflammation and ballooning versus patients with simple steatosis, as well as by longitudinally assessing regression of NASH on LSI or bariatric surgery (BS), which can lead to NASH resolution14,15.

Using systems biology and experimental approaches, we set out to identify gene sets associated with NASH presence and activity at baseline and whose expression is normalized in patients displaying reduced ballooning and lobular inflammation on LSI. We identify a NASH transcriptomic signature strongly enriched in genes controlling immune inflammatory processes, antigen presentation and cytotoxic cells. We further show that NASH activity is associated with altered blood immune cell populations, including cDC subsets and cytotoxic CD8 T cells. Moreover, in an obesity-driven mouse model of NASH exhibiting profound liver inflammation and hepatic damage, we also find increased hepatic expression of genes from the NASH transcriptomic signature, and altered hepatic populations of cDC and CD8 T cells. These results show that distinct immune cell populations play an important role in NASH activity and therefore constitute targets for NASH therapy.

Results

NASH associates with a hepatic immune-related gene module

Using transcriptome data16, we first searched for groups of genes with hepatic expression patterns linked with NASH in a large cohort of obese patients with or without histologically proven NASH (Supplementary Table 1). To this end, weighted gene co-expression network analysis (WGCNA), a method allowing identification of clusters (modules) of co-expressed genes with similar expression patterns in different experimental conditions, was used17. WGCNA of the approximately 30% most variable hepatic transcripts across patients identified nine co-expressed gene modules (Fig. 1a,b and Supplementary Table 2). To assess whether these gene modules linked with NASH, we first analysed alterations in the modules’ expression in response to improvement in NASH activity in patients with NASH at baseline and in whom a liver biopsy was available 1 year after Roux-en-Y gastric bypass (RYGB) (n = 21, Supplementary Table 3). In agreement with the reported effects of RYGB on NAFLD15, steatosis, lobular inflammation and ballooning were significantly reduced in these patients at 1 year after RYGB (Supplementary Table 3). Interestingly, RYGB resulted in a significant alteration in transcriptional activity of four gene modules at follow-up compared to baseline (Fig. 1c). To decipher which gene modules were associated with common mechanisms of NASH regression, rather than with transcriptional alterations in the liver due to body weight loss, we tested these modules in LSI responders (patients with NASH at baseline with decreased lobular inflammation and/or ballooning at 1 year follow-up, for details see Methods) and LSI non-responders in terms of NASH activity reduction, but who showed similar body weight loss at 1 year follow-up (Supplementary Table 4). Among the four gene modules affected by RYGB, only module ‘blue’ displayed significant downregulation at follow-up versus baseline in LSI responders, whereas its expression pattern hardly changed in non-responders (Fig. 1c,d). As overall gene expression levels in module ‘blue’ were similarly decreased in patients with RYGB and LSI responders, but not in LSI non-responders (Fig. 1c,d), this module is associated with reduction of NASH activity, independent of body weight changes. Indeed, the change in module ‘blue’ expression was not correlated with percentage body weight change in either the LSI or BS groups (Supplementary Fig. 1b).

Fig. 1: Identification of hepatic transcriptomic signature of NASH.
figure 1

a, WGCNA performed with 11,784 transcripts in the liver in patients with or without histologically proven NASH (n = 155 patients, see Supplementary Table 1). Clustering of co-expressed genes in nine gene modules. b, Number of transcripts in each gene module. c, Overall transcriptional regulation of gene modules on RYGB (n = 21 patients), LSI in responders (n = 10 patients) and LSI in non-responders (n = 10 patients) at 1 year follow-up compared to baseline (see Supplementary Tables 3 and 4). P values are calculated by mean-rank gene set test using the geneSetTest function as described in detail in Methods. d, Volcano plots of average log2(fold changes) versus P values (paired moderated t-test using limma package) of all transcripts (grey dots) and transcripts from gene module ‘blue’ (red dots) in patients with RYGB (n = 21 patients), LSI responders (n = 10 patients) and LSI non-responders at 1 year follow-up compared to baseline. e, Top hallmark pathways enriched in gene module ‘blue’ (n = 786 transcripts), calculated using GSEA software as described in detail in Methods. f, Venn diagrams with transcripts in gene module ‘blue’ downregulated (P < 0.05 by moderated paired t-test using limma package) in three groups of patients at follow-up versus baseline; immune-related genes are shown. g, Spearman correlations between NASH activity index and hepatic expression levels of genes in gene module ‘blue’ in the 155 patients at baseline; top genes with maximal positive correlation coefficients are shown. h, CXCL9, CXCL10 and LYZ expression (by microarray) in patients with NAFL (n = 22 patients) and NASH (n = 106 patients) at baseline. Data are presented as median with first and third quartiles as the box edges. *P < 0.05, **P < 0.01 by unpaired two-sided Mann–Whitney U-test.

Full size image

Module ‘blue’ included 786 transcripts (Supplementary Table 2) and was significantly enriched for inflammation-related pathways, such as complement, tumour-necrosis factor-α (TNF-α) and interleukin (IL)-6 signalling, as well as Kristen rat sarcoma virus (KRAS) signalling, coagulation and apoptosis (Fig. 1e). Moreover, fasting plasma insulin and C-reactive protein (CRP) levels were also associated with module ‘blue’ expression, highlighting the close link between insulin resistance, systemic inflammation and NASH (Supplementary Fig. 1c). Among the 786 transcripts in module ‘blue’, 507 were significantly downregulated at 1 year follow-up, including 195 transcripts significantly downregulated in both patients with RYGB and LSI responders, but not in LSI non-responders, with other transcripts being preferentially decreased in patients with RYGB (93 genes) or LSI responders (187 genes) (Fig. 1f). Most of these genes were related to immune system functions, including pro-inflammatory cytokines and chemokines, cytotoxic cells, infiltration of immune cells into tissues, and major histocompatibility complex (MHC) I and II antigen presentation (Fig. 1f). Interestingly, among key interferon-α (IFN-α)-responsive genes (IFIT1, OAS1, MIX1 and ISG15), only MX1 showed a modest, non-significant association with NAFLD severity at baseline and IFIT1 was significantly increased and ISG15 decreased only after RYGB (Supplementary Fig. 1d,e).

To determine whether module ‘blue’ is associated with NASH and its severity in terms of activity at baseline, we correlated the expression levels of its transcripts with the NASH activity index (sum of lobular inflammation and ballooning score) in the 155 patients. Importantly, many genes from module ‘blue’, including immune-related genes, associated positively with activity index (Fig. 1g), suggesting links with NASH disease activity. Moreover, some of these genes were also decreased in patients with RYGB or LSI responders at 1 year follow-up (Fig. 1f). In addition, hepatic expression levels of some of these immune-related genes, for instance chemokine (C-X-C motif) ligand (CXCL) 9 and 10 and lysozyme (LYZ), were significantly higher in patients with NASH versus those with NAFL at baseline (Fig. 1h). Taken together, these results indicate that NASH regression associates with the response of a specific gene module, containing multiple co-regulated genes involved in inflammation, antigen presentation, cytotoxic response and activation of T cells.

Blood immune cell signatures of NASH

As the gene module ‘blue’ identified the immune system as a key player in NASH regression on LSI, we next assessed whether blood immune cells correlate with the presence of NASH and with its activity by performing a deep immunophenotyping analysis in patients without (n = 17) or with NASH (n = 21), stratified for T2D, a major risk factor for NASH (Supplementary Table 5). All patients were obese, with the highest body mass index (BMI) in the no-NASH T2D group. Patients with T2D were also older than patients without T2D. Although lobular inflammation was absent, the level of steatosis was already higher in the no-NASH T2D group compared to no-NASH no-T2D group and a few patients featured some ballooning (Supplementary Table 5), indicating that early stages of hepatocyte damage are already present in these patients.

Flow cytometry analysis of blood immune populations was performed in these 38 patients and correlations assessed between the 39 measured immune cell populations and clinical parameters (Fig. 2a and Supplementary Fig. 2). Unsupervised hierarchical clustering grouped immune cell populations based on the similarity of their correlation patterns with clinical parameters for NASH and glucose metabolism, and yielded three main clusters: Cluster 1 associated with NASH but showed weak or absent correlations with T2D; Cluster 2 positively associated with NASH and T2D; and Cluster 3 negatively correlated with most parameters specific for NASH and T2D (Fig. 2a). Within Cluster 1, NASH activity positively correlated with natural killer (NK) cells, atypical CD16++ monocytes and HLA-DR+CD123−CD11c+CD141− cells, similar to classical dendritic cells type 2 (cDC2). Likewise, TH1 lymphocytes and NKT cells were positively associated with lobular inflammation, ballooning and NASH activity index. Immune cell populations within Cluster 2 were positively associated with lobular inflammation, ballooning, and glucose or HbA1c levels, thus linking NASH activity and T2D (Fig. 2a). Cluster 2 included pro-inflammatory CD16+ monocytes and IL-10+ CD4 T lymphocytes, the latter probably reflecting a compensatory anti-inflammatory response to increased disease activity, as previously described in mice18. Interestingly, within Cluster 2 we also found activated and cytotoxic CD8 T lymphocytes correlating with NASH activity (Fig. 2a). CD8 T cells appear to be functionally linked to hepatic damage and inflammation in mouse models of NASH10,19,20, but have not been well studied in human NASH. Cluster 3 immune cell populations were mostly negatively associated with NASH and glucose parameters. Among them are TH2 lymphocytes, including IL-5+ cells, as well as regulatory T cells (Treg). While cDC2 were positively associated with NASH (Cluster 1), both HLA-DR+CD123−CD11c+CD141+ cDC1 and HLA-DR+CD123+ plasmacytoid DC (pDC) were inversely correlated with NASH and glucose levels (Cluster 3). This indicates that pDC, cDC1 and cDC2 may play opposing roles in NASH. These results show that both NASH and T2D are associated with pronounced changes in blood immune cell populations, which may contribute to the high incidence of NASH in T2D.

Fig. 2: Correlations between blood immune cell populations, disease activity in NASH and genes in module ‘blue’.
figure 2

a, Hierarchical clustering of correlation coefficients in 38 patients (see Supplementary Table 5) between proportions of blood immune cell populations and histological liver parameters (Spearman’s correlation), T2D-associated parameters (Pearson correlation) and systemic inflammation markers (Pearson correlation). Asterisks indicate P < 0.05 for the given correlation. b, Pearson correlations between selected immune cell populations in blood from a subset of 29 patients (see Supplementary Table 5) and hepatic expression levels of genes from module ‘blue’. NAS, NAFLD activity score; AI, activity index; EM, effector memory; CM, central memory; EMRA, effector memory CD45RA+.

Full size image

Blood immune cell populations correlate with liver gene signature

Analysis of hepatic transcriptome and blood immune cell populations suggest that interplay between antigen-presenting cells (probably cDC), CD4 T cells and cytotoxic lymphocytes is linked with the presence and activity of NASH. We thus probed for correlations between blood immune cell populations and expression levels of genes from the NASH transcriptomic signature (module ‘blue’) in the 28 patients with both liver microarray and blood flow cytometry data. CD16+ monocytes correlated with elevated expression of immune and stress genes from module ‘blue’ (Fig. 2b). cDC1 and cDC2 oppositely linked with NASH activity (Fig. 2a), and cDC1, but not cDC2, displayed negative associations with hepatic expression of genes involved in immune regulation and antigen presentation (Fig. 2b and data not shown). Interestingly, elevated IL-10+ CD4 T cells were linked with increased hepatic expression of genes related to inflammation, chemotaxis and cytotoxic responses (Fig. 2b), again indicative of a possible compensatory increase of immunoregulatory IL-10 expression in CD4 T cells in NASH. Finally, we tested correlations between genes from module ‘blue’ and cytotoxic CD8 T cells expressing perforin (strongly associated with NASH activity and T2D) (Fig. 2a,b). This immune cell population showed associations with hepatic genes related to cytotoxic and IFN-γ responses (GZMA, CD226, IRF1), T helper differentiation (ITK) and TNF-α signalling (TNFAIP2) (Fig. 2b). Taken together, these results show a remarkable similarity and associations between hepatic immune pathways predicted by gene module ‘blue’ and the repertoire of NASH-associated immune populations in the blood. These data suggest a cross-talk between immune cells and hepatocytes, probably via cytokines, and direct infiltration of immune cells into the liver.

Diet-induced NASH drives increased liver cDC2 and CD8 T cells

Next, we assessed whether lifestyle-induced NASH leads to hepatic accumulation of the immune cell populations identified in the clinical study, such as cDC and cytotoxic cells. Since disease-inducing protocols are unethical and flow cytometric analysis of hepatic immune populations in humans is challenging due to the limited tissue access, we studied liver immune cells in a diet-induced model of murine NASH. Most published NASH models used to analyse hepatic immune cells do not develop obesity2 and hepatic immune populations have not yet been characterized in two published mouse models of obesity-associated NASH21,22. The diet-induced NASH model recapitulated the main known dietary drivers of this disease in humans23,24 and displayed hallmarks of well-established NASH, such as steatosis, hepatocyte damage and lobular inflammatory infiltrates (Supplementary Fig. 3 and Supplementary Information). Flow cytometric analysis of mouse liver showed unaltered neutrophil population and decreased proportions of B lymphocytes and CLEC4F+ Kupffer cells, whereas inflammatory macrophages and monocytes were increased in livers of mice with NASH (Supplementary Fig. 4a–f). Hepatic NK, NKT and CD4 T cells were not affected on NASH induction, whereas hepatic FOXP3+ Treg cells tended to be elevated (Supplementary Fig. 4g–j).

Systems biology analysis in humans predicted an important role of hepatic antigen-presenting cDC and cytotoxic cells in NASH regression on LSI (Fig. 1f). Interestingly, mouse CD172a+ cDC, phenotypically similar to human cDC2, were significantly increased in mouse livers on NASH induction (Fig. 3a,b). In contrast, NASH induction in mice led to a decrease of hepatic XCR1+ cDC (analogous to human cDC1) (Fig. 3a,c), driving a marked decrease in the cDC1/cDC2 ratio (Fig. 3d). This shift in cDC populations was specific to the liver, as no changes were found in these cDC subsets in the spleen (data not shown). Moreover, the proportion of CD8 T cells was significantly increased in livers of mice with NASH (Fig. 3e,f), suggesting that hepatic accumulation of CD8 T cells may be linked with antigen presentation by CD172a+ cDC2. These results in a lifestyle-induced NASH model corroborate the hepatic NASH gene expression and blood immune signatures in humans and indicate that hepatic accumulation of CD8 T cells and disturbance in hepatic cDC populations are the main immune hallmarks of NASH activity and progression.

Fig. 3: A diet-induced NASH model alters cDC and CD8 T cells and inflammation in the liver.
figure 3

Male C57BL/6J mice were fed conventional diet (CD) or NASH-diet (ND) for 24 weeks (see Supplementary Information). a, Representative flow cytometry plots of cDC in the liver: proportions of XCR1+ and CD172a+ of total cDC are shown (n = 8 mice CD; n = 6 mice ND). b, CD172a+ cDC2 cells as a proportion of CD45+ cells (n = 8 mice CD; n = 6 mice ND). c, XCR1+ cDC1 cells as a proportion of CD45+ cells (n = 8 mice CD; n = 6 mice ND). d, Ratio of cDC1/cDC2 cells (n = 8 mice CD; n = 6 mice ND). e, Representative flow cytometry plots of TCR-β+ T cells in the liver: proportions of CD4+ and CD8+ of total TCR-β+ T cells are shown. f, The proportion of CD8+ T cells of CD45+ immune cells in the liver (n = 9 mice CD; n = 20 mice ND). g, Quantitative PCR analysis of inflammatory gene expression in mouse livers (n = 9 mice CD; n = 20 mice ND). Data are shown as mean ± s.e.m. Statistical significance of differences between groups was analysed by unpaired two-sided t-test. mRNA, messenger RNA.

Full size image

Finally, we studied hepatic expression of inflammatory cytokines and chemokines in the diet-induced NASH model. NASH induction increased expression of Tnf, macrophage colony-stimulating factor (Csf)1, Il23a, Il33, C-C motif chemokine ligand (Ccl)2 and Cxcl10 (Fig. 3g), with the latter four genes present in the gene expression signature of NASH responsive to intervention (Fig. 1f). IFN-α-responsive genes Ifit1 and Oas1a were also induced in NASH-diet- compared to conventional-diet-fed mice (Supplementary Fig. 3g). These data indicate that this diet-induced NASH model in mice affects the same hepatic immune populations and inflammatory pathways that are modified on regression of NASH activity induced by LSI in humans.

Activated CD8 T cells are elevated in NASH and T2D in humans

Given that both gene expression analysis in human livers and flow cytometric analysis in mouse livers revealed an association of cytotoxic CD8 T cells with obesity-related NASH, we next focused on these cells to better understand their potential role in NASH in humans. Indeed, both NASH and T2D enhanced the proportions of IFN-γ+ and TNF-α+ circulating CD8 T cells (Fig. 4a). Cytotoxic functions of CD8 T cells were also increased both in NASH and T2D as evidenced by increased expression of granzyme A and B, and perforin (Fig. 4b,c). Because T2D is a major risk factor for NASH, its effects towards CD8 T cells may link these diseases. Indeed, in the absence of NASH (no-NASH T2D group), T2D was associated with elevated expression of IFN-γ, TNF-α and cytotoxic molecules in CD8 T cells and some of these patients already displayed steatosis and ballooning. These results show that both patients with NASH and T2D have a blood signature of increased populations of activated and cytotoxic CD8 T cells, which may link these diseases.

Fig. 4: NASH and T2D alter activity of cytotoxic CD8 T cells.
figure 4

a, Proportion of IFN-γ+ and TNF-α+ CD8 T lymphocytes in blood from patients with/without NASH and/or T2D (see Supplementary Table 5). Groups: n = 7 patients with no NASH and no T2D; n = 6 with NASH and no T2D; n = 7 patients with no NASH and T2D; n = 7 patients with NASH and T2D. Data are shown as median with first and third quartiles. Statistical significance of differences between groups was analysed by unpaired two-way ANOVA for effects of NASH and T2D followed by Tukey’s post-hoc test. b,c, Representative flow cytometric plots (b) and proportions of perforin, granzyme A and B expression in blood CD8 T lymphocytes from patients (c). Groups: n = 7 patients with no NASH and no T2D; n = 6 patients with NASH and no T2D; n = 7 patients with no NASH and T2D; n = 7 patients with NASH and T2D. Data are shown as median with first and third quartiles. Statistical significance of differences between groups was analysed by moderated paired t-test or two-way ANOVA (for effects of NASH and T2D) followed by Tukey’s post-hoc test.

Full size image

Hepatic CD8 T cells correlate with NASH in humans

Finally, we investigated whether the increase of circulating cytotoxic CD8 T lymphocytes in obese patients with NASH is accompanied by an accumulation of these cells in the liver. CD8+ cells with lymphoid morphology were detected on sections from all groups of patients (Fig. 5a). While both NASH and T2D were associated with significantly increased levels of CD8 T lymphocytes in the liver, no further increase was observed in patients with NASH and T2D (Fig. 5b). CD8 T lymphocytes localized within inflammatory foci in the parenchymal compartment in close proximity to steatotic and ballooned hepatocytes (Fig. 5c). Liver CD8 T lymphocyte numbers showed relatively little association with steatosis, but significantly correlated with lobular inflammation, ballooning and activity index (Fig. 5d). No significant correlation was found between liver CD8 T lymphocytes and glucose levels or HbA1c (data not shown), suggesting that the T2D-associated accumulation of CD8 T lymphocytes is not driven by dysfunctional glucose metabolism. Notably, liver CD8 T lymphocytes significantly correlated with blood CD8 T lymphocytes and with the subset of perforin-expressing CD8 T cells (Supplementary Fig. 5a). Hepatic CD8 T cells also significantly correlated with expression of perforin (PRF1) and granzyme A (GZMA) in the liver, two cytotoxic genes from the module ‘blue’, but not granzyme B (GZMB) and granulysin (GNLY) (Supplementary Fig. 5b). Finally, we found that hepatic CD8 T cells positively correlated with expression levels of multiple genes from module ‘blue’, including genes related to T cells (CD2, CD226), cytotoxic responses (KLRC4-KLRK1) and MHCII-mediated antigen presentation (HLA-DQB1) (Fig. 5e). Interestingly, the strongest correlation with hepatic CD8 T cells was observed for PTPN22 (protein tyrosine phosphatase, non-receptor type 22), which controls T-cell receptor (TCR) responsiveness and is associated with inflammatory and autoimmune diseases in humans25. Taken together, these results indicate that disease activity (lobular inflammation and ballooning) in NASH is associated with the accumulation of CD8 T lymphocytes in the liver. Moreover, hepatic CD8 T cells are linked to an increased expression of genes from the signature of active NASH. Thus, hepatic gene expression and immune signatures of NASH activity regression reveals CD8 T cells as a potential target to control hepatic inflammation, cytotoxic responses in the liver and NASH activity.

Fig. 5: Hepatic CD8 T lymphocytes correlate with lobular inflammation, ballooning and transcriptomic signature of NASH.
figure 5

a, Representative immunostaining for CD8 (red) with haematoxylin counterstaining on liver biopsies from patients with/without NASH and/or T2D. b, Quantification of CD8+ cells per mm2. No NASH, no T2D n = 10; NASH, no T2D, n = 10; no NASH, T2D, n = 7; NASH, T2D, n = 9. Data are shown as median with first and third quartiles. c, Localization of CD8 T lymphocytes (red) near immune infiltrates, steatosis and ballooned hepatocytes (indicated by arrows) in the liver from a patient with NASH. Scale bar, 50 μm. d, Correlations between hepatic CD8 T lymphocyte number and histological features in the liver (n = 36). e, Pearson correlations and −log10(P values) between hepatic CD8 T lymphocyte and expression levels of genes from module ‘blue’ (n = 29, Supplementary Table 5). Statistical significance of differences between groups was analysed by unpaired two-way ANOVA (for effects of NASH and T2D) followed by Tukey’s post-hoc test.

Full size image

Discussion

The global epidemic of NASH is an important health problem, as effective pharmacological approaches to treat NASH are not yet available. LSI and BS are the best available strategies to reduce NASH activity in some patients14,15,26. Here, we identify a hepatic transcriptomic signature of NASH in humans that is distinct from NAFL and responsive to RYGB and LSI, suggesting a shared mechanism of NASH regression by these different weight-loss methods. Interestingly, weight loss alone was apparently not sufficient to improve NASH in all patients as weight loss was similar in LSI responders and non-responders. The reversible transcriptomic signature of NASH is highly enriched by immune genes such as cytokines and chemokines and their receptors, genes involved in activation of lymphocytes by antigen-presenting cells and genes linked to cytotoxic cells. Some of these genes have already been described as associated with NASH. For example, CD44, a molecule mediating leucocyte recruitment into the liver, plays a role in methionine-choline deficient diet (MCD)-induced NASH in preclinical models27. Although the MCD model has obvious limitations in representing human NASH within the framework of the metabolic syndrome, it reproduces NASH in its hepatic phenotype and is therefore an interesting model to study intrahepatic changes in NASH. Furthermore, the same study also found elevated CD44 in human NASH biopsies. The pro-inflammatory cytokine CXCL10 also drives NASH in the MCD NASH model and is increased in blood from patients with NASH28. Here we show that downregulation of hepatic pro-inflammatory pathways accompanies regression of NASH activity on intervention in patients. Moreover, these pathways increase on lifestyle induction of NASH in a mouse model. Thus, the reversible transcriptomic signature of NASH reveals new biological targets for specific therapy of NASH focusing on its activity. While the present study focuses on NASH improvement driven by weight loss or lifestyle modification, effective pharmacological treatments could target this immune signature and act in synergy with other metabolism-focused pathways. Because the transcriptomic analysis here was performed on total RNA from liver biopsies, it is conceivable that single-cell RNA analysis may reveal more subtle differences in hepatocyte versus non-parenchymal cell responses in NASH.

Increasing evidence demonstrates a close interaction between the immune system and metabolism in the development of NASH29. Our data are in line with recent publications on associations of different immune cells with NAFLD in humans8,9,30,31. However, our study includes a global analysis of circulating immune cell populations in relation to the presence of NASH and its severity. Integrative analysis revealed a species-conserved phenotype of NASH implicating antigen-presenting cDC subtypes and cytotoxic CD8 T cells. Whereas cDC1 are thought to present antigens to CD8+ T cells and cDC2 present antigens to CD4+ T cells32,33, we surprisingly found that blood cDC1 correlate inversely with NASH activity, whereas blood cDC2 and CD8+ T cells both positively correlate with NASH. Similarly, CD172a+ cDC2 replace XCR1+ cDC1 in mouse livers on NASH induction. Moreover, elevated hepatic expression of multiple genes involved in antigen presentation in NASH is reversible on intervention. In the MCD model, total cDC have been shown to protect against NASH34. However, whether these effects are mediated by cDC1 and/or cDC2 subpopulations is unclear. As cDC subtypes control pro-inflammatory and tolerogenic immune responses depending on the tissue environment35, a shift between cDC1 and cDC2 in mice and humans might be associated with elevated hepatic inflammation and hepatocyte damage in NASH. Interestingly, depletion of CD141+ cDC1 has been shown to occur in human liver on inflammation36, suggesting an inflammation-mediated mechanism of suppression of cDC1 in relation to NASH activity. Among the hepatic myeloid populations identified in mice (such as, Kupffer cells, monocyte-derived macrophages and DC), the present gating strategy was unable to exhaustively phenotype these clearly important immune populations. Further studies are thus necessary to better understand the evolution of the hepatic immune milieu and confirm our findings in human NASH biopsies.

Previous studies demonstrated associations of CD8 T cells with liver damage in mouse models of NASH10,19,20. We show that, in patients, circulating and hepatic cytotoxic CD8 T lymphocytes are significantly linked to histological hallmarks of NASH, such as lobular inflammation and ballooning. Moreover, CD8 T lymphocytes correlate with hepatic markers of inflammation and antigen presentation in NASH. Polymorphisms in one of these genes, HLA-DQB1, a MHCII haplotype, is associated with NAFLD37,38. Another gene, PTPN22, is an important regulator of TCR signalling and linked with several autoimmune diseases25. Furthermore, in the lifestyle-induced NASH model, elevated hepatic CD8 T cells associated with CD172a+ cDC2, suggesting a possible role of antigen presentation and TCR activation of CD8 T cells in NASH. CD8 T lymphocytes appear to contribute to insulin resistance in a mouse model of diet-induced obesity38. Similarly, hepatic CD8 T lymphocytes and type I IFN signalling promote glucose intolerance and insulin resistance in mice39,40. Interestingly, the reversible transcriptomic signature of NASH that we identified to be enriched by cytokines and cytotoxic molecules in circulating CD8 T cells is linked to NASH activity and T2D in humans. Curiously, dysregulated glucose metabolism in patients with T2D does not explain the increased numbers of CD8 T cells in the liver, suggesting that other factors contribute to the T2D-mediated accumulation of hepatic CD8 T cells. Moreover, in addition to elevated hepatic CD8 T cells, some patients with T2D, without an unequivocal NASH diagnosis, already display ballooning. Together, these results suggest a role of CD8 T cells in the interplay between T2D and NASH. We found that the reversible transcriptomic signature of NASH is enriched by apoptosis pathway and NF-κB target genes (data not shown). These hepatic NF-κB target genes are downregulated on NASH activity reduction accompanied by decreased ballooning. Ballooning is a poorly understood but essential hallmark of NASH, representing a specific form of ‘undead cells’ with features of initiated but not resolved apoptosis6. Because the apoptotic machinery is controlled by the NF-κB signalling pathway41, downregulation of the pathway on LSI or BS might resolve apoptosis and thereby eliminate ballooned cells from the liver. Moreover, associations of the blood immune signature of NASH, particularly cytotoxic CD8 T cells, with ballooning suggest a cross-talk between the immune system and stressed hepatocytes in NASH. Although our results identify a pathway contributing to NASH resolution in a manner independent of body weight, it is clear that metabolic control through direct action on the liver or via extra-hepatic organs also contributes to NAFLD. Indeed, adipose tissue may directly contribute to NAFLD progression through systemic cytokine and immune signalling and drive hepatic injury in the context of obesity39,42,43.

Our study provides insights into the molecular mechanisms implicated in NASH and its regression on intervention. In patients and a mouse model, we comprehensively and consistently show surprisingly pronounced associations of immune pathways and cell populations with NASH activity. These pathways, as well as the exact molecular targets revealed by this study, underline an important role of innate and acquired immunity in the development and severity of NASH, which may be targeted for treatment.

Methods

Patients

All patients were consecutively recruited at the Liver Clinic and Obesity Clinic of the Antwerp University Hospital and underwent hepatologic and metabolic work-ups. Blood analysis included blood cell count and white blood cell formula. Exclusion criteria were alcohol consumption >2 units per day for women and >3 units per day for men, liver diseases other than NAFLD, age <18 years and liver cirrhosis. For the baseline gene expression analysis in the liver, patients (n = 155, including no liver disease n = 27, simple steatosis n = 22, NASH n = 106) were selected from a cohort recruited since 2006 (ref. 16). Selected obese patients with NASH with paired biopsies at 1 year follow-up (LSI n = 20, RYGB n = 21) were included for gene expression analysis (Supplementary Fig. 1). For immunophenotyping analysis, patients (n = 38) were consecutively recruited between 2014 and 2016. Gene expression analysis was performed on a subset of these patients (n = 29) who were included in the baseline gene expression analysis (Supplementary Fig. 1). The study protocol is part of the Hepadip protocol (Belgian registration number B30020071389) and was approved by the Ethical Committee of the Antwerp University Hospital (file 6/25/125). Written informed consent was obtained from all patients.

Selection of patients

From patients visiting the Obesity Clinic at the Antwerp University Hospital, who were recruited from 2006 to 2016, 155 patients with hepatic RNA microarray and clinical data available at baseline were enrolled16,44. Among the 155 patients, 41 obese non-diabetic patients with NASH at baseline were selected for comparative baseline and 1 year follow-up hepatic RNA microarray and clinical data analysis. Of these 41 patients with NASH at baseline, 31 patients displayed an improvement in disease activity, with decreases of lobular inflammation and/or ballooning, 1 year after RYGB (n = 21) or LSI (n = 10), whereas 10 patients did not improve 1 year after LSI. We defined the LSI patients with NASH at baseline who improved ballooning and/or lobular inflammation at 1 year follow-up as responders to intervention and those who did not as non-responders to intervention. In a second independent cohort, immune cell populations were analysed in consecutive patients with four clearly distinct phenotypes enrolled based on their metabolic (T2D) and histological (NASH) phenotype (total n = 38): 17 patients without NASH in which some degree of simple steatosis was allowed (without (n = 7) and with (n = 10) T2D) and 21 patients with unequivocal NASH based on histological parameters (without (n = 11) and with (n = 10) T2D). The group without T2D and NASH consisted of a slightly lower number of patients, as liver biopsies were only performed on clinical indication of the potential presence of NASH. High-quality RNA appropriate for microarray analysis was obtained from liver biopsies in 29 out of 38 patients (13 patients without NASH and 16 patients with NASH).

Clinical assessment and biological measurements

Fasting blood glycaemia was analysed in the morning. A 2-h oral glucose tolerance test (75 g of glucose), including insulin quantification, was performed. Homeostatic Model Assessment for Insulin Resistance (HOMA-IR) was calculated45. BMI, glucose and percentage of HbA1c were measured46. Liver biopsies were performed for suspected NAFLD as indicated by elevated serum transaminase levels or a steatotic liver on ultrasound as described previously44.

Histological assessment of the liver biopsies

All liver biopsies were stained (haematoxylin-eosin, Sirius red, reticulin and Perls’ iron stain) and scored by two expert pathologists blinded to all clinical information. Histological features of NAFLD (steatosis, ballooning, lobular inflammation and fibrosis) were assessed using the NASH Clinical Research Network Scoring System criteria47. The NAFLD activity score was calculated as the sum of steatosis, lobular inflammation and ballooning scores. An activity index was also calculated as the sum of ballooning (range 0–2) and lobular inflammation (range 0–3) in line with recent observations of the distinct roles of steatosis versus activity of disease5. NASH was defined by the simultaneous presence of steatosis ≥1 AND ballooning ≥1 AND lobular inflammation ≥147. As outlined, the sum of ballooning and lobular inflammation was calculated and indicated as the activity score.

Mouse model of diet-induced NASH

Wild-type male C57BL/6J mice (8 weeks of age) were purchased from Charles River Laboratories (France). Mice were maintained in a pathogen-free environment (12:12 h light/dark cycle, 21 °C–24 °C) and were given ad libitum access to water and food. Alternatively, Foxp3-YFP reporter mice48 maintained on a C57BL/6J background were used to monitor hepatic Treg populations. Littermate animals were randomized by body weight before the start of the diet. No power calculations were performed to determine sample size. Mice were fed either a control diet (standard rodent chow, 5% kcal fat) or a ‘NASH’ diet (45% kcal fat, 40% kcal carbohydrate, 15% kcal protein with 1% (by weight) cholesterol; SAFE diets, France) for 24 weeks. All experiments were performed following approval by the Ethics Committee for Animal Experimentation from Nord-Pas de Calais Region (APAFIS#5746-2016040109244171 and APAFIS#7160-2017040313471173).

Histology and immunohistochemistry

Human liver biopsy sections were incubated for 1 h at 22 °C with rabbit monoclonal anti-CD8a (SP16) (Thermo Fisher Scientific), followed by staining with a mouse anti-rabbit HRP-conjugated antibody, and revealed with Vector ImmPRESS HRP Reagent Kit and Vector NovaRED Substrate Kit. Sections were counterstained by haematoxylin. Mouse liver samples were fixed with 4% paraformaldehyde, embedded in paraffin and stained with haematoxylin-eosin. Images were obtained by conventional microscopy using a Nikon Ti-U microscope.

Flow cytometry

EDTA blood samples from patients were collected and peripheral blood mononuclear cells (PBMC) isolated using Percoll. For cytokine staining, cells were incubated with 20 ng ml–1 phorbol 12-myristate 13-acetate (PMA), 1 μg ml–1 ionomycin and 1 μg ml–1 brefeldin A in Roswell Park Memorial Institute (RPMI) 1640 medium for 4 h at 37 °C, 5% CO2. Then 1 × 106 PBMC were pre-incubated with Fc Block to minimize non-specific binding and they were labelled with conjugated antibodies CD3 (UCHT1, PE-CF594), CD4 (OKT4, APC-CY7), CD45RA (HI100, V500), CD56 (HCD56, BV605), CD8 (HIT8a, PE-CY7), CCR10 (6588-5, PE), CD11b (ICRF44, AF700), CD11c (3.9, PB), CD123 (6H6, PE-CY7), CD127 (A7R34, BV605), CD14 (HCD14, PerCP-CY5.5), CD141 (B-A35, FITC), CD16 (3G8, V500), CD161 (HP-3G10, PB), CD172a (15-414, APC), CXCR3 (G025H7, APC), CD19 (HIB19, PE-CF594), CCR6 (R6H1, PE), CD197/CCR7 (3D12, FITC), HLA-DR (G46-6, AF700 or APC-CY7), Perforin (dG9, PE), Granulysin (DH2, AF647), Granzyme A (CB9, AF700), Granzyme B (GB11, PB), TNF (MAb11, AF700), IFN-γ (4S.B3, BV421), IL-10 (JES3-9D7, eFluor660), IL-17A (BL168, AF700), IL-22 (22URTI, eFluor660), IL-5 (TRFK5, APC). Intracellular staining was performed using the Cytofix/Cytoperm kit (BD Biosciences) according to the manufacturer’s instructions.

Immune cells were isolated from mouse liver, digested 45 min with collagenase D, using centrifugation at 1,000g with 30% Percoll. Cells were treated with Zombie UV to discriminate live and dead cells, and then incubated with Fc Block and labelled with conjugated antibodies: CD45 (BUV737, clone 104), CD11b (BUV395, clone M1/70), CCR2 (BV421, clone SA203G11), Ly6C (BV785, clone HK1.4), F4/80 (BV711, clone BM8), NK1.1 (AF700, clone PK136), CD4 (BV605, clone RM4-5), CD8a (BV510, clone 53-6.7), Ly6G (PE-Cy7, clone 1A8), IA/IE (BV650, clone M5/114.15.2), CD11c (APC-Cy7, clone N418), CD19 (PE-CF594, clone 1D3), TCRb (APC, clone H57-597). For myeloid cell staining, an additional panel was used with the following antibodies: CD45 (PE-CF594, clone 30-F11), CD172a (BUV395, clone P84), CCR2 (BV421, clone SA203G11), B220 (BV510, RA3-6B2), CD26 (BV605, clone H194-112), XCR1 (BV650, clone ZET), F4/80 (BV711, clone BM8), Ly6C (BV786, clone HK1.4), MHCII (FITC, clone M5/114.15.2), CD64 (PE, clone x54-5/7.1), CLEC4F (PE-Cy7, clone 370901), CD19 (APC, clone eBio1D3), CD3 (AF700, clone 500A2), NK1.1 (AF700, clone PK136), CD11c (APC-Cy7, clone N418). The CLEC4F antibody was coupled to PE-Cy7 using the Abcam Antibody Coupling Kit (reference ab102903) and used at a final dilution of 1:100.

Flow cytometry analysis was performed on a BD LSRFortessa X-20 (Becton Dickinson). Results were acquired with the Diva software (Becton Dickinson) and analysed using FlowJo software (Tree Star). Additional details are provided in the Reporting Summary.

RNA extraction

For human samples, the entire biopsy was homogenized for RNA extraction, purification and processing as described previously44. Total RNA was isolated from mouse liver using Trizol reagent and used for PCR with reverse transcription and real-time PCR. RNA from human liver biopsies was isolated using RNeasy Mini Kit (Qiagen) and used for microarray analysis. Quantification and RNA integrity number (RIN) was tested using the Agilent 2100 Bioanalyzer System. Only RNA samples with RIN ≥ 6 were used for microarray analysis.

PCR with reverse transcription and real-time PCR

Total RNA (500 ng) isolated from mouse liver was treated with DNase I and used to generate complementary DNA with a high-capacity cDNA reverse transcription kit. Gene expression was measured by SybrGreen-based real-time PCR. Results were normalized to the normalization factor calculated as average expression of housekeeping genes Ppia and Tbp, and the ΔΔCt method was used for all real-time PCR analyses. Primers sequences are provided in Supplementary Table 6.

Microarray analysis

Transcriptome analysis was performed with the Affymetrix GeneChip HuGene 2.0 ST arrays16,49. All liquid handling procedures were performed on a GeneChip Fluidic Station 450. GeneChips were scanned with a GeneChip Scanner 3000-7G (Affymetrix) using Command Console v.4.1.2. Quality controls were performed using the Affymetrix expression console.

Microarray data processing and WGCNA

Microarray data were normalized by the robust multi-average method50 using oligo/Bioconductor and corrected for batch effects using SVA/Bioconductor R packages51. In total, 38,598 annotated transcripts were selected for analysis, and 11,784 transcripts with maximal variability across all patients at baseline (n = 155) based on median absolute deviation were selected for WGCNA and tested using the WGCNA R package17. Biweight midcorrelations and weighted adjacency matrix were calculated using the power threshold of 11, which was selected based on the scale-free topology fit model. Gene modules were identified using the ‘hybrid’ method and parameters deepSplit = 4 and mergeCutHeight = 0.15. For creation of volcano plots, log10(P values) and log2(fold changes) were calculated using the limma R package52. Selected sets of transcripts of gene modules were tested for differential expression with geneSetTest from the limma R package53. The Gene Expression Omnibus repository accession number for the microarray data is GSE106737. Gene set enrichment analysis of gene modules was performed using GSEA software (http://software.broadinstitute.org/gsea/).

Statistical analysis

No statistical method was used to predetermine the sample size. For animal studies, mice were randomized by body weight before dietary challenge and no blinding was performed for subsequent analysis. For comparison in patients at baseline and 1 year after BS or LSI, the paired moderated t-test or the paired Mann–Whitney U-test was used. For comparison of the four groups of patients, data were analysed by two-way analysis of variance (ANOVA), using NASH and T2D as factors, followed by Tukey’s post-hoc test for multiple comparisons. For histological quantification of hepatic CD8+ cells, the observer was blinded from the clinical parameters of each patient. Parametric Pearson correlation (continuous data) or non-parametric Spearman’s rank-order correlation (categorical data) and hierarchical clustering were performed using R. Gene networks were visualized using Cytoscape. Values of P < 0.05 were considered significant. P values were adjusted using the Benjamini–Hochberg procedure. Statistical analyses were performed with Prism 6 (GraphPad Software, Inc.) or R v.3.4.4.

Reporting Summary

Further information on research design is available in the Nature Research Reporting Summary linked to this article.