A prognostic molecular signature of hepatic steatosis is spatially heterogeneous and dynamic in human liver

Abstract

Hepatic steatosis is a central phenotype in multi-system metabolic dysfunction and is increasing in parallel with the obesity pandemic. We use a translational approach integrating clinical phenotyping and outcomes, circulating proteomics, and tissue transcriptomics to identify dynamic, functional biomarkers of hepatic steatosis. Using multi-modality imaging and broad proteomic profiling, we identify proteins implicated in the progression of hepatic steatosis that are largely encoded by genes enriched at the transcriptional level in the human liver. These transcripts are differentially expressed across areas of steatosis in spatial transcriptomics, and several are dynamic during stages of steatosis. Circulating multi-protein signatures of steatosis strongly associate with fatty liver disease and multi-system metabolic outcomes. Using a humanized "liver-on-a-chip" model, we induce hepatic steatosis, confirming cell-specific expression of prioritized targets. These results underscore the utility of this approach to identify a prognostic, functional, dynamic "liquid biopsy" of human liver, relevant to biomarker discovery and mechanistic research applications.

Keywords: diabetes; liver-on-a-chip; metabolic dysfunction-associated steatotic liver disease; non-alcoholic fatty liver disease; proteomics.

Graphical abstract

Figure 1

Proteins related to hepatic steatosis are primarily expressed in the liver and identify pathways of metabolism (A) Volcano plot of proteins associated with hepatic steatosis after adjustment for age, gender, race, and BMI. For visualization, proteins with an FDR < 5% in CARDIA derivation subsample are visualized with the beta coefficient and p values presented coming from models using the CARDIA validation subsample. (B) Heatmap of the top 25 positively associated and top 25 negatively associated proteins with hepatic steatosis in the CARDIA validation sample. MASLD is defined as CT liver attenuation <40 HU. (C and D) Tissue expression analysis at the transcriptional (C) and tissue protein (D) level for proteins related to hepatic steatosis in CARDIA using the full SomaScan 7k platform as the background demonstrated enrichment of proteins expressed in the liver. (E) KEGG and Reactome pathway analysis. (F) Hub gene analysis of significant proteins associated with liver attenuation showing the hub genes (≥5 connections; circles) and all proteins with high confidence connections to the hub genes (rectangles). p values are from 2-sided tests.

Figure 2

Development of a proteomic score of hepatic steatosis and its relation with clinical outcomes (A) A protein score of liver attenuation by CT (less attenuation ∼ more steatosis) demonstrated moderate correlation with the parent variable in both CARDIA derivation and validation samples. (B) Replication of the association between a protein score of liver attenuation and MRI-based measure of hepatic steatosis (proton density fat fraction: higher ∼ more steatosis, opposite directionality as with CT-based liver attenuation) in UK Biobank. (C) The protein score is related to controlled attenuation parameter (higher value ∼ more steatosis) in CCHC. (D) Forest plot of associations with clinical outcomes in UK Biobank along with C-index comparisons of models with and without the protein score (see Table S15). p values reported are for comparisons of C-indices with two-sided tests. Error bars represent 95% confidence intervals.

Figure 3

Transcriptional architecture of the hepatic steatosis proteome (A) Uniform manifold approximation projection (UMAP) of integrated single-nuclear and single-cell RNA sequencing in steatotic and healthy liver, colored by a composite expression score (derived from gene expression of implicated target proteins, n = 19, single-nuclear RNA-seq n = 5, single-cell RNA-seq n = 14; see STAR Methods). (B and C) UMAP of spatial transcriptomic data (Visium, 5 samples from 4 individuals) in the liver colored by liver pathology diagnosis (healthy vs. fatty, B) and the composite expression score (C). (D) Violin plot comparing composite expression score across Visium spots by fatty vs. healthy state (Wilcoxon rank-sum test). (E) Representative images of healthy and steatotic hematoxylin-eosin-stained liver tissue overlaid with Visium spots, colored by composite expression score, demonstrating increased activity of implicated targets in steatotic regions. Images are from livercellatlas.org (Guilliams et al.³⁴). All sections presented in our parent manuscript are shown in Figure S4. (F) Composite expression across liver zonation groups. (G) Differential expression of implicated targets between healthy and fatty liver (Visium) versus circulating proteomic regression coefficient. A more positive proteomic coefficient indicates less liver fat, and a more positive log₂ fold change indicates greater expression of a given transcript in healthy (non-steatotic) liver. Highlighted in purple are targets that were considered as differentially expressed using spatial data (see text). For targets with multiple aptamers (e.g., IGFBP2 has multiple SomaScan aptamers), we present the mean regression coefficient for that target. (H) Gene expression of significantly different implicated targets between healthy and steatotic regions (Visium), liver zonation (Visium), and across cell types (single-cell RNA sequencing).

Figure 4

Transcriptional heterogeneity of spatial targets in human liver across steatosis stages (A) Bulk transcript log₂ fold change in human liver (over control samples without histologic steatosis) for those genes (among 33 significant on spatial studies) that were significantly differentially expressed in at least one comparison (stage 1 vs. control; stage 2 vs. control; stage 3 vs. control). Of the 33 genes passed forward for assessment in bulk transcriptomics, 12 were not significantly expressed in any of the stages of steatosis (by FDR adjusted p < 0.05) and were not included in visualization. The “liver attenuation beta” represents the regression coefficient against liver attenuation in the CARDIA derivation sample. A positive coefficient (red) indicates a greater protein level is related to higher attenuation (lower steatosis); a negative coefficient (blue) indicates a greater protein level is related to lower attenuation (higher steatosis). For proteins with multiple aptamers (e.g., IGFBP2 has multiple SomaScan aptamers), we present the mean regression coefficient for that protein. This analysis excluded individuals with stage F4 fibrosis, given differences in hepatic physiology at this stage of decompensation. (B) Violin plots of example gene expression (in log₂ counts per million) for genes that displayed a “concordant” directionality between the proteome and the bulk transcriptome (top and middle) and “discordant” directionality between the proteome and the bulk transcriptome. Differential gene expression analysis was performed using limma-voom with the protein-coding genes using an FDR of 5% (Benjamini-Hochberg). ∗p < 0.05, ∗∗p < 0.01, ∗∗∗p < 0.001, ∗∗∗∗p < 0.0001.

Figure 5

Transcriptional and proteomic architecture of MASLD on a humanized liver-on-a-chip (A and B) The structure and experimental design of MASLD induction on the liver-on-a-chip (LOC). (C) Successful MASLD model generation on a representative LOC. On the left, lipid droplet accumulation was visualized after 5-day treatment period of FAs. Representative bright-field and fluorescent confocal images of the LOC cells (scale bar, 100 μm). DAPI (nuclear) and LipidSpot (lipid droplet) stains are shown. Red arrows represent lipid droplet accumulation. (D) mRNA expression of canonical genes implicated in steatosis demonstrates expression patterns consistent with MASLD in both hepatocytes and NPCs. A total of 6 chips were included (3 FAs and 3 control). Results were analyzed by an unpaired t test and expressed as mean ± standard error of 3 independent experiments. Each data point represents the average of 3 technical replicates. Control is in blue and FA treated is in red. Relative expression is shown as fold change (delta-delta CT) relative to control, normalized to β-actin expression. (E) mRNA expression of top genes on the LOC that were prioritized by proteomic and transcriptomic studies (see text). Breaks in y axis are presented given disparate expression of some genes (e.g., HMGCS1 had low expression in non-hepatocytes, while HSPA1A was expressed at low levels in hepatocytes). ME1, CTSZ, and DEFB1 were not expressed in hepatocytes; CDA, SHBG, IL1RAP, and IGFBP2 were not expressed in NPCs. See text for details. (F) Secretory protein expression of top targets (both cell types: HMGCS1, SERPINE1, and PYGL; NPC specific: DEFB1 and CTSZ; hepatocyte specific: ACY1 and AKR1C4) on the LOC was consistent with the mRNA expression. A total of 8 chips were included (4 FAs and 4 control). Results were analyzed by an unpaired t test and expressed as mean ± standard error of 4 independent experiments. Each data point represents the average of 2 technical replicates. Control (Ctrl) is in green and FA treated (FA) is in pink. Abbreviations: ne, not expressed (raw Ct > 40); nd, not detected; ns, non-significant. ∗p < 0.05; ∗∗p < 0.01; ∗∗∗p < 0.001; ∗∗∗∗p < 0.0001.