Multi-omic network analysis identified betacellulin as a novel target of omega-3 fatty acid attenuation of western diet-induced nonalcoholic steatohepatitis

Abstract

Clinical and preclinical studies established that supplementing diets with ω3 polyunsaturated fatty acids (PUFA) can reduce hepatic dysfunction in nonalcoholic steatohepatitis (NASH) but molecular underpinnings of this action were elusive. Herein, we used multi-omic network analysis that unveiled critical molecular pathways involved in ω3 PUFA effects in a preclinical mouse model of western diet induced NASH. Since NASH is a precursor of liver cancer, we also performed meta-analysis of human liver cancer transcriptomes that uncovered betacellulin as a key EGFR-binding protein upregulated in liver cancer and downregulated by ω3 PUFAs in animals and humans with NASH. We then confirmed that betacellulin acts by promoting proliferation of quiescent hepatic stellate cells, inducing transforming growth factor-β2 and increasing collagen production. When used in combination with TLR2/4 agonists, betacellulin upregulated integrins in macrophages thereby potentiating inflammation and fibrosis. Taken together, our results suggest that suppression of betacellulin is one of the key mechanisms associated with anti-inflammatory and anti-fibrotic effects of ω3 PUFA on NASH.

Keywords: betacellulin; docosahexaenoic acid; multi-omic network; nonalcoholic steatohepatitis; ω3 PUFA.

Figure 1. NASH mouse model outlines all expression/omics data, DHA/EPA treatment effects on outcome

1. All omics data collected and shown in Fig 1A bar graph and its associated table was acquired using samples (N = 8/treatment group) from a preclinical NASH prevention study previously described (Depner et al, ,; Lytle et al, 2015). The comparison of the western diet + olive oil (WD + O) group vs. reference diet (chow, RD) group showing the differentially expressed genes or parameters with a P‐value < 0.05 and FDR < 10% and those that have treatment effect reversal uniquely by DHA (blue), or EPA (orange) or similar with significance in both EPA and DHA is EPA&DHA (green).

2. Heatmap of differentially expressed genes in WD + O vs. RD fed mice and organized by treatment effect category: DHA, EPA, or EPA & DHA. Data shown are the geometric mean expression for each gene per each treatment group. Row max is displayed as red, row min is displayed as blue.

3. Heatmap of parameters (lipids and metabolites) in WD + O vs. RD fed mice and organized by treatment effect category (DHA, EPA, or DHA &EPA) that are significant at P‐value < 0.05 and FDR < 0.1. Data shown are a geometric mean of each parameter measurement for each treatment group. Row max: red, row min: blue.

4. Shown are selected representatives (Cd36, Notch2, Cx3cr and Egfr) from each category according to treatment effect (top left: DHA, top right: EPA, bottom left: EPA & DHA, bottom right: no treatment effect). Shown here is an example of one gene from each category (Data are mean ± SD, N = 8 mice/treatment group; One‐way ANOVA, with multiple comparisons test with WD + O, ns (not significant), *P < 0.05, **P < 0.001, ***P < 0.005, ****P < 0.0001).

5. Scatterplot of fold change differences between WD + O and DHA or EPA treated mice with number of genes regulated by DHA (Blue) or EPA (Orange) displayed (Pearson's Chi‐squared test, ****P < 0.0001).

6. Top enriched biological process changed by DHA treatment. Top + Blue bars: Induction by DHA, Bottom + Green bars; Repression by DHA. Bars show −log₁₀(P‐value) transformed values for visualization.

Figure EV1. Extent of DHA reversal effects are significantly higher than those of EPA

1. A–C

   The Log₂FC of EPA & DHA gene expression (A) and lipids and metabolites (B, C) show the extent of DHA reversal effects as significantly higher than EPA though similar in profile. The up and down regulation are shown in separate plots for clarity (A, B) (paired, two‐sided t‐test, ns [not significant], ****P < 0.0001).

2. C

   Scatterplot of fold change differences between WD + O and EPA (x‐axis) or DHA (y‐axis) treated mice with number of lipids and metabolites regulated similarly by DHA & EPA displayed (Pearson's Chi‐squared test, **P < 0.004).

3. D

   Heatmap of differentially expressed genes in individual mice (WD + O vs. RD fed mice and organized by treatment effect category: DHA, EPA, or EPA & DHA). Data shown are the log transformed quantile normalized expression for each gene per prevention group. Row max is displayed as red, row min is displayed as blue.

Figure EV2. Interrogation of network model from the multi‐omic NASH data

1. An outline with the steps involved in deriving, analyzing and interrogation of network model from the multi‐omic NASH data (see Materials and Methods; GE: gene expression; P: Phenotype/Anthropometric/Biochemical data; LM: Lipids/Metabolites; PUC: Proportion of unexpected correlation). The distribution of treatment effects from the NASH preventive model among the network parameters is shown in the right panel.

2. Outline of cell–cell interactions used to calculate the interaction BiBCs (see Materials and Methods). The nodes represent all genes part of cell types as a cluster in the network and edges are the interaction strength among the nodes. The node size indicates the number of genes represented by the cell type; color chart is proportional to the treatment effects in each cell type.

3. Box plot of maximum cell–cell interaction BiBC for the genes belonging to each treatment effect category (DHA [blue], EPA&DHA [green], EPA [orange] and no category [gray]). From the network cell–cell BiBC analysis, genes regulated by DHA, DHA&EPA have higher BiBC that indicates a higher contribution to cell–cell communication than genes regulated by EPA.

4. Bar plots for abundance of top BiBC lipids, shown are the PGs (cardiolipin precursors) and SM in NASH preventive study (Data are mean ± SD, N = 8 mice/treatment group; Ordinary One‐way ANOVA, multiple comparisons test with WD + O, ***P < 0.005, ****P < 0.0001).

Figure 2. Multi‐omic network (NW) reconstructed to model NASH in vivo using the data from preventive model (see Fig EV2; Materials and Methods)

1. The cytoscape visualization of the network has nodes (circles, rectangles) representing genes, lipids, metabolites, plasma biochemical, and anthropometric data (Depner et al, 2013a), and edges representing correlation in a color ranging light red to light blue depending on correlation (1 to −1). The nodes are colored based on their treatment effect category membership, with DHA (blue), EPA&DHA (green), EPA (orange) and no category (gray). Network clusters are based on infomap modules additionally characterized by gene and lipids functional enrichment.

2. Bar plot of number of NW genes from each treatment category (DHA [blue], EPA [orange], or EPA&DHA [green]) with assignment to a given cell type. Subplot: figure shows t‐SNE plot with all cell type clusters from a reanalyzed NASH mouse single cell RNA‐seq dataset (Xiong et al, 2019) used in our study to assign cell type information. NAM‐ Nash associated macrophages, KC‐Kupffer like cells, DC‐ dendritic cells, HSC‐ hepatic stellate cells and E‐ endothelial cells.

3. Violin plot of average cell–cell interaction BiBC for the genes belonging to each treatment effect category (DHA [blue], EPA&DHA [green], EPA [orange] and no category [gray]) shown with n = 141, 718, 1,299, 4,007 genes respectively; solid lines indicate median, dashed line—quartiles; P values ** < 0.01; **** < 0.0001.

4. The scatterplot shows the network cell–cell BiBC verses node degree for the top hepatic lipids in the NASH network. The figure insert is the structural representation of the lipid TG58:11 (TG 16:0/20:5/22:6) with DHA and EPA as two of its acyl chains from the NASH prevention study. These top hepatic lipids are shown with treatment effect category EPA&DHA (green).

5. Bar plots for abundance of top BiBC lipids, shown are the triglycerides with DHA and EPA as acyl chains in NASH preventive study (Data are mean ± SD, N = 8 mice/treatment group; Ordinary One‐way ANOVA, multiple comparisons test with WD + O, ns [not significant], *P < 0.05, **P < 0.001, ***P < 0.005, ****P < 0.0001).

Figure 3. Combining NASH network with meta‐analysis of liver cancer identifies betacellulin as a key pathogenic regulator downregulated by DHA

1. The heatmap of genes from human liver cancer meta‐analysis ordered by corresponding genes from the mouse NASH on the effect of EPA and DHA on prevention and treatment of NASH in mice (Depner et al, ; Lytle et al, 2017).

2. 3D scatter plot showing gene set enrichment analysis (GSEA) for the DHA effects with Rank score on the x‐axis, GSEA −Log₁₀(P‐value) on the y‐axis, and cell–cell interaction BiBC on the z‐axis. Relevant pathways are labeled in the figure with the (BTC)‐ERBB pathway (red) ranked highly by all metrics.

3. Scatterplot of network BiBCs between gene expression and anthropometric parameters plotted against cell–cell BiBCs. Members of the Btc‐Erbb pathway genes are overlayed and indicate the importance of Btc‐Erbb signaling in the NASH multi‐omic network model with preventive effects. Each circle is a node in the network, filled circles (preventive effects) with red outer circle are part of the Btc‐Erbb pathway.

4. Bar plot of Btc gene expression in mouse livers from the prevention and treatment studies (N = 5 or 6/group) experiment. Data are mean ± SD, Ordinary One‐way ANOVA, multiple comparisons test of reference diet (chow, RD), western diet + DHA (WD + DHA) with western diet + olive oil (WD + O), *P < 0.05, **P < 0.001, ***P < 0.005, ****P < 0.0001).

5. BTC gene expression from meta‐analysis of human normal and liver cancer datasets indicates a significant increase in liver cancer samples (Effect size FDR = 0.0011). The size of each individual dot represents the number of patients associated with the dataset with higher the number of patients darker the dot color, Normal (green) and Cancer (red; 7 transcriptomic datasets with a total of 544 tumor, 260 non‐tumor, and 32 healthy patient liver samples).

6. BTC gene expression in human datasets of NASH. Left: Control and NASH liver (GSE193080) with normal 3 patients and Nash with 32 patient samples (Welch's t‐test, One‐tailed, **P < 0.001). Middle: GSE48452, Control, N = 14, Obese Healthy, N = 27, Steatosis, N = 14, Nash, N = 18 (Mann–Whitney test, Two‐tailed, *P < 0.05). Right: normalized Btc gene expression (GSE96971), NASH patients before and after 1 year of treatment with EPA and DHA, N = 9 (Paired t‐test, One‐tailed, *P < 0.05).

Figure EV3. Liver cancer meta‐analysis and network interrogation points to significance of BTC

1. Outline for human liver cancer meta‐analysis for the orthologous genes representing the NASH mouse model. Right panel shows the heatmap of individual human datasets for selected genes in the meta‐analysis (with treatment effects in NASH model).

2. The bar graph of gene set enrichment analysis using GSEA (KEGG pathways) for the human cancer meta‐analysis genes with DHA effects in mouse models. Data are displayed as −log₁₀(P‐value).

3. The gene expression for BTC‐EGFR‐ERBB pathway and cell cycle related genes in the NASH mouse preventive and treatment models. The color scale is indicated from high expression of the genes in red to low in blue.

4. Distribution of Btc BiBC values calculated between DHA treatment reversed genes and metabolomic/lipidomic data (x‐axis) and DHA controlled genes and anthropometric data (y‐axis) in 5,000 random networks. Dark regions represent areas where a calculated Btc BiBC value is more likely to be found due to random chance. The probability of finding an actual Btc BiBC value equal to or higher than those seen in the random networks (43/5,000) is 0.009.

5. Violin plot of the BiBC of Btc between DHA reversed genes and anthropometric data in 5,000 random networks compared to its actual BiBC value in the reconstructed network (one sample Wilcoxon test P < 1 * 10⁻¹⁵).

Figure EV4. BTC and other Ligand‐Receptor interaction network with DHA/EPA effects in NASH preventive mouse model

1. A

   NASH Ligand–Receptor network with DHA/EPA effects in NASH preventive mouse model. The individual cells expressing the genes inferred from single cell RNA sequence reanalysis are overlayed and cell specific ligand–receptor interactions are shown. Nodes are shown as diamonds (ligands), or hexagon (receptors) and edges are the correlation between the nodes. Node color is according to the treatment effect. The thickness of the node is the cell–cell interaction BiBC, higher the better and is shown thicker. Each cell type and their respective genes are colored differently as labeled. The genes belonging to other cells in addition to specific cells are shown as part of ‘Not assigned’ group.

2. B

   The normalized Btc expression form others experiments (Left panel; Control and Nash fed with high fat high sugar mouse model; GSE197884; Mann–Whitney test, One‐tailed, N = 3). (Right panel) The normalized Btc expression form other datasets (GSE222576; Control and CCl4 treated liver fibrosis mouse model; Ordinary one‐way ANOVA, *P < 0.05, **P < 0.005).

Figure 4. Hepatic stellate cell (LX2) proliferation could promote fibrosis via BTC‐TGFβ‐2, reversed by DHA

1. LX2 cell proliferation assay after treatment with BTC, 0–25 ng/ml (left panel, N = 6 individual experiments). Fold change in DNA concentration after treatment with BTC 0–10 ng/ml (right panel, N = 3 experiments). Green bar in each plot is untreated control. Data are displayed as mean ± SD with each point being an individual experiment; Ordinary One‐way ANOVA, multiple comparisons test with control vehicle, ns (not significant), *P < 0.05, **P < 0.001, ***P < 0.005, ****P < 0.0001).

2. Collagen staining (Pico Sirius Red, PSR) in LX2 cells indicating increased fibrosis (collagen production) when stimulated with BTC (20 ng/ml; N = 5 experiments) or a vehicle control. A representative image of stained cells is shown (left panel). A bar plot of staining intensity normalized by total protein (μg) per well (Data are displayed as mean ± SEM, N = 5 experiments; unpaired, two‐sided t‐test **P = 0.0058).

3. Gene enrichment (biological process) of LX2 cells (see Materials and Methods) while induced by BTC (BTC vs. Vehicle; one‐sided t‐test) identifies pathways reversed by DHA treatment in vivo. Data are displayed as −log₁₀(P‐value).

4. Bar plot of TGFB2 gene expression in LX2 cells treated with vehicle or BTC (20 ng/ml; Data are displayed as mean ± SD, N = 5 experiments) in vitro (paired, two‐sided t‐test *P < 0.05, FDR < 0.1).

5. Tgfb2 gene expression in vivo. DHA reversed the gene expression significantly in the in vivo experimental model both in Preventive & Treatment models (Data are displayed as mean ± SD, N = 8 mice/group [preventative study] or N = 5 or 6 mice per group [treatment study]; Ordinary One‐way ANOVA, multiple comparisons test of reference diet [chow, RD], western diet + DHA [WD + DHA], western diet + EPA [WD + EPA] each with western diet + olive oil [WD + O], ns [not significant], *P < 0.05, **P < 0.001, ***P < 0.005, ****P < 0.0001).

6. TGFB2 gene expression from meta‐analysis of human normal and liver cancer data indicates a highly significant increase in expression in the liver cancer samples (Effect size FDR = 0.006). The size of each individual dot represents the number of patients associated with the dataset with higher the number of patients darker the dot color, Normal (green) and Cancer (red; 7 transcriptomic datasets with a total of 544 tumor, 260 non‐tumor, and 32 healthy patient liver samples).

7. Gene enrichment analysis of genes upregulated in TGFβ‐2 treated cells (GSE45382) that were reversed by DHA in vivo, identifies top enriched pathways. Data are displayed as −log₁₀(P‐value).

8. Gene enrichment analysis of downregulated genes in TGFβ‐2 treated cells (GSE45382) that were reversed by DHA in vivo, identifies highly enriched gene ontologies. Data are displayed as −log₁₀(P‐value).

Figure 5. BTC promotes TLR‐dependent inflammation and integrin production by macrophages

1. A

   Scatterplot of differential expression with BTC‐LPS‐PGN (BTC and TLR2/4 ligands) treated cells to control THP‐1 cells in vitro (BTC‐LPS‐PGN/control in THP‐1 cells; P‐value < 0.05; Details about concentrations and duration of the TLR agonists treatment are detailed in Materials and Methods) against in vivo differential expression with WDO to WD + DHA in NASH preventive model. Filled circles are the genes in DHA treatment category.

2. B

   Gene enrichment is shown for the genes upregulated from BTC‐LPS‐PGN treatment in THP‐1 cells that were reversed by DHA treatment in vivo (Gene ontology‐biological process). Data are displayed as −log₁₀(P‐value).

3. C

   Gene enrichment is shown for the genes downregulated from BTC‐LPS‐PGN treatment in THP‐1 cells that were reversed by DHA treatment in vivo (Gene ontology‐biological process and cellular components). Data are displayed as −log₁₀(P‐value).

4. D

   The heatmap for summary metric of all major molecular pathways affected by BTC treatment in THP‐1 cells shown to be prevented by DHA in both in vivo NASH preventive and treatment models (see Materials and Methods).

5. E–G

   The individual pathways enriched in BTC‐LPS‐PGN treated THP‐1 cells but reversed by DHA in vivo in the NASH preventive and treatment models displayed as summary metric (Data are displayed as mean ± SD, N = 5 experiments).

6. H

   Individual gene expression for IL‐1b and CSF1 from THP‐1 cells in vitro treated with BTC and or TLR2/4 ligands (N = 5 experiments, paired, one‐sided t‐test, ns (not significant), *P < 0.05, **P < 0.001; see Materials and Methods) and in vivo NASH preventive model (Data are mean ± SD, N = 8 mice/treatment group; Ordinary One‐way ANOVA, with Dunnett's multiple comparisons test with WD + O, ns (not significant), *P < 0.05, **P < 0.001, ***P < 0.005, ****P < 0.0001).

7. I

   The bar plots for summary metric are shown for integrin signaling and collagen pathway in THP1 cells (Data are displayed as mean ± SD, N = 5 experiments).

8. J

   The bar plot for ITGA1 gene expression is shown in THP‐1 cells in vitro treated with BTC and or TLR2/4 ligands (N = 5 experiments, paired, one‐sided t‐test, ns (not significant), *P < 0.05, **P < 0.001; see Materials and Methods) and in vivo NASH preventive model (Data are mean ± SD, N = 8 mice/treatment group; Ordinary One‐way ANOVA, with Dunnett's multiple comparisons test with WD + O, *P < 0.05, ****P < 0.0001).

Figure EV5. Potential mechanism of BTC regulation by the ω3 PUFA

1. The representation of network model and network interrogation to identify genes (nodes) top ranked according to Bipartite betweenness centrality (BiBC) and degree which potentially mediate effects of ω3 PUFA lipids on BTC‐dependent gene expression profile. The legend for nodes is shown.

2. The scheme of analysis to identify upstream regulator of Btc. The three modules of analysis with experimental data, network interrogation and identification of transcription regulator for Btc in Cholangiocytes. Identified transcription factors are labeled. Predicted binding site and motif for Foxo3 mediated transcription regulation in mice and human for Btc gene expression in the liver cholangiocytes is shown.

3. Foxo3 gene expression in vivo. DHA reversed the gene expression significantly in the in vivo experimental model both in Preventive & Treatment models (Data are displayed as mean ± SD, N = 8 mice/group (preventative study) or N = 5 or 6 mice per group (treatment study); Ordinary One‐way ANOVA, multiple comparisons test of reference diet (chow, RD), western diet + DHA (WD + DHA), western diet + EPA (WD + EPA) each with western diet + olive oil (WD + O), ns (not significant), *P < 0.05, ***P < 0.005).