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.

The paper explained Problem

Metabolic-associated fatty liver disease (MAFLD, also known as non-alcoholic fatty liver disease, NAFLD) and non-alcoholic steatohepatitis (NASH) are major risks factor for cirrhosis, hepatocellular carcinoma, and liver failure. Lack of specific FDA-approved treatments makes this a prominent health issue affecting roughly one third of the population, particularly in the Americas and South-East Asia. Although omega-3 fatty acids are known for their positive effects on the liver in NASH, a comprehensive analysis to identify the key molecular factors involved is still necessary.

Results

In this work, we used top-down system biology approach with causal inference analysis to reconstruct a multi-omic network (transcriptome, metabolome, lipidome, single cell RNA sequencing) and to reveal main molecular players responsible for beneficial effects of ω3 fatty acids on NASH and potentially on liver cancer. In agreement with our previous findings, out of two omega-3 fatty acids tested, docosahexaenoic acid (DHA) was more potent than eicosapentaenoic acid (EPA) in its effects on multi-omic readouts. In our search for NASH mechanisms that also play a role in liver cancer, we found that a key aspect of omega-3 fatty acid action involves inhibiting the expression and function of betacellulin (BTC), a less studied member of the epidermal growth factor family. We verified our network predictions in vitro and found that BTC stimulates TGFβ-2-driven collagen production in hepatic stellate cells and enhances microbial signals in the induction of integrins by macrophages. These processes work together to promote liver fibrosis and inflammation. Conversely, omega-3 fatty acids, especially DHA, can interrupt and even reverse these processes. Additionally, we identified the transcription factor FOXO3 as the most likely upstream regulator of the effects of omega-3 fatty acids on BTC.

Impact

Our research has revealed a novel mechanism through which omega-3 fatty acids regulate liver health, mitigating harmful processes during NASH. These newfound mechanistic insights have the potential to facilitate the development of innovative therapies that target the BTC pathway for NASH treatment and liver cancer prevention. Additionally, the gene expression signature triggered by BTC holds promise as a potential biomarker for guiding clinical trials involving ω3 PUFA, potentially advancing personalized medicine for liver disease management.

Introduction

Metabolic diseases associated with obesity have increased to epidemic proportions in recent years and are one of the leading causes of morbidity and mortality (Konerman et al, 2018; Dufour et al, 2022). Metabolic-associated fatty liver disease (MAFLD) or nonalcoholic steatohepatitis (NASH), type 2 diabetes and cardiovascular diseases are all associated with obesity and a sedentary lifestyle (Marjot et al, 2020; Lazarus et al, 2022). About 80 million adults and 13 million children are obese in the US alone. Among these, 60 percent of patients with body mass index > 30 have evidence of liver steatosis (Jump et al, 2015; Dufour et al, 2022). NASH is a progressive form of nonalcoholic fatty liver disease (NAFLD) and is a major risk factor for cirrhosis, hepatocellular carcinoma (HCC) and liver failure. While treatments to manage the co-morbidities associated with NASH, i.e., obesity and type 2 diabetes are available, NASH has no specific FDA-approved treatment. Thus, lifestyle modifications and dietary interventions are the current options available to NASH patients. Most if not all drugs which targeted individual molecules or specific pathways have failed to significantly improve the NASH patient (Neuschwander-Tetri, 2020; Pfister et al, 2021; Ampuero et al, 2022; Dufour et al, 2022). This strategy that has been effective in the treatment of other diseases might not be adequate for NAFLD/NASH therapy because it does not address entirely the complexity of this disease and may miss the master regulators involved in disease onset and progression.

Omega-3 polyunsaturated fatty acids (ω3 PUFA) are known to be consistently lower in livers of NASH patients when compared to healthy patients or patients with benign steatosis (Burke et al, 1999; Fridén et al, 2021). This prompted us to hypothesize that dietary supplementation with ω3 PUFA would restore liver functions. Indeed, this strategy was very successful in a preclinical mouse model, not only reducing liver steatosis but also in attenuating western diet-induced hepatic fibrosis (Depner et al, 2013a,b; Jump et al, 2015). Moreover, ω3 fatty acid treatment of children and adults with NAFLD have demonstrated these dietary lipids reduce hepatic steatosis and hepatic injury (Iannelli et al, 2013; Spooner & Jump, 2019). While it is well established that ω3 PUFA have the capacity to regulate hepatic mechanisms controlling fatty acid synthesis and oxidation, as well as inflammation, it is less clear if these pathways form the extent of ω3 fatty acid regulation of hepatic function (Jump et al, 2018). Despite several studies demonstrating the therapeutic effects of ω3 PUFA in NAFLD/NASH models the mechanism of action has been elusive. Nevertheless, it is important to note that many studies demonstrated diverse effects of ω3 in the liver ranging from immune-modulatory activities (Gutiérrez et al, 2019) and improvement of oxidative stress (Yang et al, 2019) to structural effects related to incorporation of phospholipids into the mitochondrial membrane, altogether with potentially positive impact on liver function. In addition to the effects of ω3 PUFA on hepatic cells, recent studies suggest that ω3 PUFA can alter gut microbiota, a known player in the pathogenesis of NAFLD/NASH (Watson et al, 2018). Furthermore, preclinical and clinical studies demonstrated that docosahexaenoic acid (DHA, 22:6, ω3) might be a more efficient agent than eicosapentaenoic acid (EPA, 20:5, ω3) in preventing and treating NAFLD (Depner et al, 2013a,b; Spooner & Jump, 2019).

Thus, while there are many studies describing different molecular and cellular effects of ω3 fatty acids (Hodson et al, 2017; Zöhrer et al, 2017; Musa-Veloso et al, 2018; Okada et al, 2018; Tobin et al, 2018; Šmíd et al, 2022) and some evidence that they may be an effective therapy for NAFLD/NASH, the key mechanisms of how they improve liver health is unknown. As such, we used a comprehensive unbiased systems approach to answer these questions. For this, we evaluated the liver transcriptome, metabolome and lipidome and assessed causal inferences via multi-omic network analysis to identify prospective mechanism operating in the diseased liver that were restored by EPA and/or DHA. Since NASH is one of the precursors of liver cancer, we also performed a meta-analysis of human liver cancer to evaluate which aspects of NASH pathogenesis leading to cancer are reversed by ω3 fatty acids. Together, our studies pointed to betacellulin (BTC), one of several epidermal growth factor receptor (EGFR) agonists, as a master regulatory molecule that was downregulated by ω3 PUFA in the NASH liver. We further validated the impact of BTC in cell culture experiments that established TGFβ-2 and integrins as the main downstream molecular targets of BTC in human hepatic stellate cells and macrophages, respectively. Suppression of these pathways specifically by DHA leads to attenuated fibrosis in NASH. Thus, our study disclosed an entirely novel mechanism for ω3 fatty acid control of hepatic function and its beneficial action against detrimental molecular events in the liver leading to NASH.

Results DHA reverses the effects of WD more effectively than EPA

We first performed a comprehensive analysis of molecular changes contributing to prevention of NAFLD/NASH by two ω3 PUFA, namely docosahexaenoic (DHA, 22:6, ω3) and eicosapentaenoic acid (EPA, 20:5, ω3). For this, we evaluated histological markers of NASH (Dataset EV1a), transcriptomic, metabolomic, and lipidomic changes caused by DHA and EPA in the whole tissue liver samples from Ldlr^−/− mice fed a western diet (WD) with the addition (or not) of DHA or EPA (Depner et al, 2013a; Fig 1A–C, Dataset EV1b). To focus our analysis on disease features corrected by ω3 PUFA, we first established which changes induced by WD were reversed by DHA or EPA treatment (see Materials and Methods for details). We then established four categories of parameters: (i) regulated by DHA only (e.g., Cd36, Fig 1D top left); (ii) regulated by EPA only (e.g., Notch2, Fig 1D top right); (iii) regulated similarly by DHA and EPA (e.g., Cx3Cr1, Fig 1D bottom left); (iv) not regulated by either DHA or EPA (e.g., Egfr, Fig 1D bottom right). Although there was a large overlap between the effects of each ω3 PUFA, overall, DHA showed stronger effects than EPA in restoring alterations caused by WD. Specifically, while both EPA and DHA showed similar effects on 19% of the genes affected by disease, DHA alone reversed more genes (11%) than EPA alone (3%), (Fig 1B and C). In line with this result, focusing on genes regulated by both fatty acids (DHA and EPA) we observed more pronounced changes by DHA than EPA (P < 0.0001; Figs 1E and EV1A and B). The genes upregulated by DHA were enriched for several functional categories: mitochondrial organization, translation, and energy derivation by oxidation of organic compounds were among the most prominent pathways affected. Regulation of cytokine production was the top enriched pathway among the downregulated genes (Fig 1F). Thus, the first step of transcriptome analysis indicated that DHA had a stronger effect than EPA in reversing damage inflicted by WD on the liver potentially by restoring mitochondrial function and inhibiting inflammation.

In the second largest omics data set, represented by lipidomes, we did not see differences in the number of lipids regulated by DHA or EPA and ω3 PUFA restored most lipids impaired by WD (Fig 1A and C). However, we detected a significantly stronger effect on lipids regulated by DHA vs. by EPA (Fig EV1C). Analysis of anthropometric features and plasma biochemicals also showed a more pronounced effect by DHA than EPA, as we previously reported (Depner et al, 2013a; Fig 1A). Overall, DHA demonstrated stronger effects than EPA in reversing WD-induced changes in both gene expression (Fig EV1D) and lipid concentrations.

Mapping effects of ω3 PUFA onto multi-omic network model of NASH and its cellular components

To investigate how changes in different omics data relate to each other and which of the effects are contributing to the effects of ω3 PUFA, we have reconstructed a multi-omic network model of NAFLD/NASH and mapped ω3 PUFA effects into this model (Fig EV2A). After filtering out data features that did not pass statistical (Dong et al, 2015) and causality (Yambartsev et al, 2016) criteria thresholds (see Materials and Methods), the multi-omic network consisted of 6,743 nodes connected by 80,811 edges. Specifically, the network included 6,346 gene transcripts, 5 anthropometric nodes, 11 plasma biochemicals, 357 lipids, and 24 metabolites (Fig 2A). We next identified clusters (sub-networks) and using functional enrichment analyses found that different subnetworks were enriched in different pathways, including mitochondrial organization, myeloid leukocyte activation, cell and mitochondrial membrane fluidity, remodeling, signaling, and energy metabolism. It also included processes such as macromolecule catabolic and fatty acid metabolic processes (Fig 2A, Dataset EV1b).

Multiple cell types in the liver contribute to NASH pathogenesis (Ramachandran et al, 2019; Xiong et al, 2019; Seidman et al, 2020) but which cells are responding to DHA treatment has not been comprehensively studied. Therefore, we mapped genes regulated by ω3 PUFA to cell type information using a previously published single cell RNA-seq dataset generated from diet-induced NASH mouse livers (Xiong et al, 2019). Among different cell types, most of the ω3 PUFA-regulated genes were assigned to one of the two major hepatic macrophage subpopulations including NASH-associated macrophages (NAM, 1,455 genes) and Kupffer-like cells (KC, 585 genes), named according to the previous study (Xiong et al, 2019). These were followed by hepatocytes, cholangiocytes, and hepatic stellate cells (522, 365, and 192 genes respectively; Fig 2B). In line with our initial observation (Fig 1A), we found that DHA reversed expression of a larger number of genes than EPA, irrespectively of the cell type (Fig 2B).

As a next step, we wanted to ensure that our model (i.e., multi-omic network enhanced with information about cell type assignment) is consistent with known effects of ω3 PUFA on NASH. For this, we analyzed connectivity related topological network properties, known as bipartite betweenness centrality (BiBC) that was shown by us (Morgun et al, 2015; Li et al, 2022) and others (Lam et al, 2021) as a metric reflective of causal relationships between nodes in a co-variation network. High ranked BiBC nodes represent parameters mediating impact of one part (e.g., module/cluster) of a network on another (Morgun et al, 2015; Li et al, 2022). Thus, we calculated cell–cell interaction BiBC based on the expression of genes assigned to different cells and affected or not by ω3 PUFA in NASH multi-omic network (Fig EV2B). The results showed that DHA regulated genes had higher BiBCs than those regulated only by EPA or unaffected by either fatty acid (Fig 2C; Fig EV2C). This result is in line with our previous observations of DHA being more potent than EPA in improving NASH (Lytle et al, 2017).

This result supported the general validity of the transcriptomic part of our network, thus, we asked which of the lipids/metabolites may have major contribution to cell–cell interactions in NASH. For this analysis, in addition to BiBC we accounted for node degree as this topological property has been the most prevalent network parameter in computational biology (Sorrells & Johnson, 2015; Choobdar et al, 2019) reflecting importance of the node in controlling direct and indirect neighbors and therefore corresponding biological function. Specifically, nodes of high degree (also called hubs) generally represent master regulators of a part of the large network (clusters or subnetworks) in which these hubs are situated (Sorrells & Johnson, 2015; Choobdar et al, 2019).

We selected nodes from lipids/metabolite part of the network that were top ranked by BiBC and had high degree (top 10% for both parameters) and were also regulated by ω3 PUFA. We identified five lipids which represented < 2% of total lipids detected in our data. Importantly, two out of the five lipids (Figs 2D and EV2D) were triglycerides containing EPA and DHA (TG 58:11 and TG 60:10) as acyl chains (Fig 2E). Only 7 other lipids out of 357 had similar chemical composition (i.e., TG containing EPA/DHA) were detected by our lipidomic assay, thus making this finding being random virtually impossible (P < 0.0001). Thus, this analysis identified two lipids, which are drastically depleted in the liver during and have protective effect reversing NASH (Burke et al, 1999; Jump et al, 2015; Fridén et al, 2021). Although this result might be expected, it is nevertheless important as a “positive control” providing additional confidence in our network.

The other three lipids were phosphatidyl glycerol (PG) 36:3, PG 36:2 (cardiolipin precursors) and sphingomyelin (SM) d42:2, which are both main membrane lipids (Fig EV2D). Since top inferences from this network validated key previously known molecular aspects of ω3 PUFA action in NASH liver we can rely on this model for inference of new yet undiscovered aspects of this process.

Combining the NASH network with the meta-analysis of liver cancer identifies betacellulin as a key pathogenic regulator of NASH inhibited by DHA

NASH can progress to liver cirrhosis and cancer in humans (Samuel & Shulman, 2018; Anstee et al, 2019; Pfister et al, 2021). Furthermore, this was also shown in the mouse model used in this study (Chen et al, 2021). Although ω3 PUFA have been studied as a potential prevention strategy for colorectal cancer progression and other cancers (Van Blarigan et al, 2018; Dierge et al, 2021), there is still little understanding which cancer pathways can be inhibited by these fatty acids. Thus, as a next step, we mapped the molecular model of transcriptome alterations by DHA in liver tissue (Fig 2A) to transcriptomic alterations in human liver cancer. For this, we first performed a meta-analysis of human liver cancers (7 transcriptomic datasets with a total of 544 tumor [Hepatocellular carcinoma and Cholangiocarcinoma], 260 non-tumor, and 32 healthy patient liver samples) and established which genes were expressed concordantly by cancer and WD-induced alterations in our preclinical mouse model (Fig EV3A, Dataset EV2). Among 2,080 concordantly expressed genes between NASH and cancer, 22% (456 genes), 10% (221 genes) and 3% (56 genes) were reversed by DHA and EPA, DHA alone and EPA alone, respectively (Fig 3A left panel, Dataset EV3). While the current mouse study was designed to assess the effects of ω3 PUFA on preventing NASH (Prevention study), in another study (Lytle et al, 2017) we evaluated liver transcriptomes of DHA-treated mice with already established NASH (Treatment study). In this case, we observed that treatment effects of DHA cover ~ 54% (361 out of 677genes) of its preventive effects relevant for human liver cancer (Fig 3A right panel). These results suggest that DHA can potentially be used as cancer preventive strategy in patients with already established NAFLD/NASH.

Next, we asked which of the molecular pathways regulated in cancer and reversed by DHA in mice may mediate beneficial effects of DHA. For this, we combined gene enrichment analysis with BiBC to focus on genes with the largest impact on cell–cell interactions (see details in Materials and Methods). The top enriched pathway was Oxidative Phosphorylation with the well-known pathways such as TGFβ and p53 signaling being enriched to a lesser extent (Fig EV3B). One pathway, however, that stood out as highly enriched was the ERBB signaling pathway; it was also top ranked in mediating DHA-driven cell–cell interactions (i.e., BiBC; Figs 3B and EV3C). ERBBs are known homologs of EGFR, which are activated through binding to EGF and related members of the EGF family of growth factors. These include EGF-like ligands or cytokines that are comprised of at least 10 proteins including betacellulin, transforming growth factor-alpha (TGF-α), amphiregulin, HB-EGF, epiregulin, and neuregulins and the various other heregulins (Olayioye et al, 2000; Wieduwilt & Moasser, 2008; Chen et al, 2016).

Using our multi-omic network (including DHA affected and not affected genes in NASH), we ranked genes from the ERBB pathway based on their potential capacity to mediate effects of DHA and identified betacellulin (BTC), an alternative ligand of EGFR (Fig EV4A), as a top gene among the important genes (Grb2, Gsk3 Gsk3β/α, and Cbl) in the pathway (Fig 3C). Interestingly, EGFR itself, although not regulated by DHA was the second best potential regulatory gene for this pathway in NASH. To ensure statistical robustness of BTC potential causal role reflected by its high BiBC rankings, we used two additional approaches (see Materials and Methods) that demonstrated very low probability of Btc being randomly highly ranked (with P < 10⁻¹⁵ and probability density P = 0.009, respectively; Fig EV3D and E).

We next confirmed that Btc gene expression was increased in other models of NASH, liver cancer and fibrosis including a chemically induced disease (Green et al, 2022; Hammad et al, 2023; Fig EV4B), and it was downregulated by both EPA and DHA in our prevention model, albeit DHA had a stronger effect than EPA (Fig 3D, left panel).

In agreement with mouse models, Btc gene expression was increased in human hepatic carcinoma meta-analysis (Fig 3E) and in human NASH in two studies (Ahrens et al, 2013; Fujiwara et al, 2022; Fig 3F, Dataset EV4). Importantly, EPA + DHA treatment of patients with NASH for one year resulted in a significant decrease of Btc gene expression (Fig 3F).

BTC promotes NASH fibrosis via activating hepatic stellate cells to produce TGFβ-2

Given that BTC was predicted as a new target, central for effects of DHA on NASH, we next verified which cell types express BTC and its receptor (EGFR). For this, we integrated the network with available interaction information from ligand–receptor database (Abugessaisa et al, 2021) and single cell RNA-seq data (Xiong et al, 2019; Fig EV4A). We observed that cholangiocytes were the primary population of cells expressing BTC (Fig EV4A). Although its source is restricted to cholangiocytes, BTC's role as a secreted growth factor indicates it can act on many different types of neighboring cells (e.g., hepatocytes, macrophages, hepatic stellate cells and others) that express EGFR/ERBBs (Fig EV4A and B).

Among different liver cell types which can respond to BTC and proliferate during NASH progression, hepatic stellate cells (HSCs), frequently called mesenchymal cells in humans (Ramachandran et al, 2019; Carter & Friedman, 2022) produce collagens and represent a major contributor to hepatic fibrosis. Indeed, the number of mesenchymal cells counted in humans with cirrhosis is markedly higher than those of healthy livers (Appendix Fig S1C). Importantly, DHA prevents and reverses fibrosis in a NASH mouse model (Depner et al, 2013a; Lytle et al, 2015) and decreases expression of two out of three collagen encoding genes (COL1A1, COL1A2, COL4A1) that are upregulated in liver cancer and NASH (Appendix Fig S1D). Therefore, we next tested effects of BTC on human hepatic stellate cells using the LX2 cell line (Xu et al, 2005). LX2 cells were grown and pretreated with and without BTC using EGF as a growth factor positive control. LX2 growth was significantly increased by BTC (Fig 4A) and to a similar extent as was observed for EGF (Appendix Fig S1E). Moreover, we observed increased collagen staining (Sirius red) in cells treated with BTC (Fig 4B).

We next performed a RNASeq transcriptomic analysis of LX2 cells treated with BTC and compared it to genes regulated by DHA in vivo. We found 63 genes upregulated by BTC and downregulated by DHA and 16 genes downregulated by BTC and upregulated by DHA. Among the enriched pathways for the set of genes induced by BTC and repressed by DHA were transcripts involved in cell growth including ERBB signaling (Fig 4C, Appendix Fig S1F–H, Dataset EV5). Strikingly, TGFB2 was the only gene found in common across several enriched categories. Moreover, expression of TGFB2, but not TGFB1 was increased by BTC in vitro (Fig 4D, Appendix Fig S1G), repressed by DHA in both prevention and treatment mouse studies (Fig 4E), and increased in human liver cancer meta-analysis (Fig 4F). A recent study demonstrated that TGFβ-2, but not TGFβ-1 has a critical non-redundant role in promoting lung and liver fibrosis (Sun et al, 2021). Therefore, we hypothesized that downstream effects of reduction of BTC by DHA can be explained by reduction of TGFβ-2. For this, using publicly available in vitro data on TGFβ-2 effects and our in vivo data, we evaluated which genes regulated by DHA were regulated in the opposite direction by TGFβ-2 (Dataset EV6). We found 62 genes upregulated by TGFβ-2 and downregulated by DHA and another 62 genes downregulated by TGFβ-2 and upregulated by DHA. DHA downregulated/TGFβ-2-upregulated genes were highly enriched for production of collagen trimers and extracellular matrix organization (Fig 4G) while genes downregulated by TGFβ-2 and reversed by DHA were enriched for mitochondrial inner membrane and other mitochondrially related functions (Fig 4H). Altogether these results suggest that one of a key mechanism of fibrosis reduction by DHA is achieved through an inhibition of BTC and a consequent reduction of HSCs proliferation and TGFβ-2-induced collagen production (Sun et al, 2021).

TLR-dependent inflammatory processes in NASH are exacerbated by BTC

Our results from stellate cells support that inhibition of BTC by DHA would explain reduction of fibrosis and improvement of mitochondrial function. However, the reduction in inflammatory pathways and macrophage gene expression, another major effect of DHA we observed (Figs 1 and 2), could not be explained by effects of BTC on stellate cells. Furthermore, EGFR deficiency specifically in macrophages has been shown to attenuate liver cancer in a mouse model (Lanaya et al, 2014). Finally, we have previously reported that systemic levels of TLR2 ligands are decreased by DHA (Lytle et al, 2015; Appendix Fig S2A). These results, along with a well-known fact that a genetic deficiency of microbiota sensors (TLR2 and TLR4) attenuates NASH (Spruss et al, 2009; Miura et al, 2013; Wu et al, 2020), indicate that a reduction in TLR2 and/or TLR4 ligand levels by DHA might interact with the reduction of BTC in preventing the disease. While different receptors for BTC are widely distributed across different liver cell types (Appendix Fig S1B and S2B), TLR2 and TLR4 are predominantly expressed by macrophages in the liver (Appendix Fig S2C).

Taken together, our next question was: which processes regulated by DHA in the liver can be explained by the effects of BTC and TLR2/4 agonists on macrophages? We also asked whether BTC modulates TLR2/4-dependent immune stimulation in macrophages. To answer these questions, we differentiated the human monocyte cell line (THP-1) to a macrophage-like phenotype and stimulated with BTC, LPS (TLR4 ligand) and PGN (TLR2 ligand) and compared global gene expression in these cells to cells stimulated with TLR2/4 ligands only or unstimulated control cells (Appendix Fig S2D, Dataset EV7).

To investigate a potential interaction effect, we tested a range of concentrations of BTC and TLR-agonists evaluating expression of IL-8, and CCL2 (MCP-1) expression, well-known targets of BTC (Lanaya et al, 2014; Shi et al, 2014) and TLR-agonists (Seki et al, 2007; Spruss et al, 2009; Miura et al, 2013; Wu et al, 2020) and chose the lowest doses that induces their expression (Appendix Fig S2E and F). We next performed transcriptomic analysis of THP-1 cells treated with BTC-LPS-PGN and compared it to genes regulated by DHA in vivo. We found 179 genes upregulated by BTC-LPS-PGN and downregulated by DHA and 285 genes downregulated by BTC-LPS-PGN and upregulated by DHA (Fig 5A). Genes upregulated by BTC-LPS-PGN were enriched for several pathways related to monocyte/macrophage related immune functions, the cell cycle, and collagen binding. Among the most enriched categories for the downregulated genes were mitochondrion, NAD metabolic process, and endoplasmic reticulum membrane (Fig 5B–G).

To assess the relative contribution of BTC and of TLR2/4 agonists to the observed combined functional effect, we calculated a summary metric for each pathway (see Materials and Methods) and compared their values between each treatment group (Fig 5D).

To check if the expected effects of BTC are present in macrophages, we first verified expression of EGFR/cell cycle and epithelial-mesenchymal transition (EMT) pathways and observed their increase among the categories upregulated by BTC alone and in combination with TLR agonists (Fig 5E). Analysis of the downregulated pathways showed that in combination with TLR ligands BTC inhibited expression of genes involved in mitochondrial functions (TCA cycle, oxidative phosphorylation, Appendix Fig S2G) and NAD metabolic functions (critical pathways operating in mitochondria; Samuel & Shulman, 2018; Simões et al, 2018; Xie et al, 2020; Fig 5F).

As for the upregulated pathways, we observed diverse immune functions such as ‘innate immune response’, ‘regulation of interferon-gamma signaling’ with a similar or slightly enhanced expression of genes when BTC was added with TLR2/4 agonists (Fig 5G, Appendix Fig S2H). However, some genes showed clear interactions between BTC and TLR-agonists. For example, IL1B was upregulated by TLR2/4 agonists and slightly increased by BTC, but CSF1, a classical factor of macrophage growth and proliferation (Hume & MacDonald, 2012), was upregulated only when both BTC and TLR2/4- agonists were present (Fig 5H).

Among immune-related pathways, the strongest BTC effect either alone or in combination with TLR2/4 agonists was on the integrin-mediated signaling pathway, which partially overlapped with collagen binding (Fig 5I). Interestingly, while there was a trend for BTC increasing transcript levels of ITGα6 and ITGα9 and TLR2/4 agonists of ITGα1, only the combination of BTC with TLR-agonists significantly induced expression of all three integrins (Fig 5J, Appendix Fig S2I). Notably, the integrin pathway was not regulated by TLR2/4 agonists alone suggesting a possible unique role of the crosstalk between EGFR and TLR pathways in controlling fibrosis-related molecular function in macrophages.

Potential mechanism of BTC regulation by ω3 fatty acids

As a last step of this study, we sought to identify potential upstream regulators (transcription factors, TF) that could be targets of ω3 PUFA in control of betacellulin and its downstream processes. For this, we first established BTC-downstream profile by searching in the network genes that were regulated by BTC in vitro (Figs 4C and 5B–D) and by ω3 PUFA in the NASH mouse model (Fig 2A) in the opposite directions. We found 150 genes (62 up and 88 down regulated by BTC) and defined them as ω3-controlled BTC-regulated genes. Next, we investigated which genes (and which transcription factors among them) might mediate effects of ω3 PUFA on BTC-regulated genes. First, we ranked genes using the network degree and BiBC between ω3 PUFA and BTC-regulated genes (Fig EV5A). Among the 41 top-ranked genes, we found three transcription factors with two of them (Foxo3 and Kdm5d) expressed in cholangiocytes (the only expressing BTC in the liver). Further investigation of binding motifs identified that only Foxo3 had a motif in the Btc promoter (Fig EV5B), suggesting that regulation of this TF by ω3 PUFA can be implicated in the regulation of BTC. Interestingly, Foxo3 is a well-known suppressor of oncogenesis (Liu et al, 2018; Tsuji et al, 2021), and we found that this TF is downregulated by WD and upregulated by ω3 PUFA in mice (Fig EV5C).

Discussion

For many metabolic diseases such as diabetes (Jermendy et al, 2018; Dahlén et al, 2021), atherosclerosis (Bhatt et al, 2019), and obesity (Jastreboff et al, 2022) there are highly efficient drugs that treat and/or prevent development of these diseases. The only frequent representative of metabolic diseases that still lack pharmacological agents that would pass clinical trials is NAFLD/NASH. As this disease is often called “fatty liver”, most of attempted treatments target reduction of liver fat (Rinella, 2015). However, these treatments cannot resolve liver fibrosis, which is a more resistant aspect of the pathogenesis of this disease (Wattacheril et al, 2018). Fibrosis is the main cause of liver failure in patients with NASH and also precedes and leads to development of liver cancer (Anstee et al, 2019; Pfister et al, 2021). Therefore, revealing the mechanisms of action of DHA that reduce fibrosis in a preclinical mouse model may aid in the rational use of ω3 PUFA or help in developing new drugs that would act upon the same cellular/molecular targets (Jump et al, 2015).

In this study, we identified betacellulin (BTC), one of the less studied ligands of EGFR, as a master regulator whose reduction by DHA potentially leads to prevention/treatment of fibrosis. Indeed, we revealed that BTC induces the TGFβ-2, a critical contributor to liver fibrosis via collagen production by hepatic stellate cells (Dropmann et al, 2016). Moreover, in combination with TLR 2/4-agonists (also reduced by DHA), BTC induces integrin pathway in macrophages, the cell type in the liver most affected by DHA treatment and well-known to be involved in pathogenesis of fibrosis in different organs (Wynn & Vannella, 2016). Thus, BTC represents a candidate master regulator inducing two most important factors (collagens and integrins) contributing to liver fibrosis and consequently promoting liver cancer.

In addition to its effect on fibrosis, reduction of BTC seems to be also mediating another important effect of DHA, which is improvement of mitochondrial function-related pathways.

Indeed, mitochondrial damage is widely reported in NAFLD/NASH and its improvement by DHA is clearly seen even at the initial stage of our analysis (Fig 1E) indicating a strong impact of DHA on this pathway. Accordingly, BTC in combination with microbiota derived stimulation (represented by TLR agonists) has a negative effect on the expression of genes involved in mitochondrial functions.

Another robust effect of DHA observed at the initial step of our analysis was the inhibition of inflammation. Not surprisingly, a deeper investigation into this phenomenon led us to potential effects of DHA on microbiota-related molecules and on hepatic macrophages involved in sensing microbes. In fact, network analysis combined with scRNAseq data pointed to macrophages as main cellular targets of DHA in the liver (Fig 2B). Macrophages are also the primary cells in the liver that express both TLR 2/4 (Appendix Fig S2C) whose microbiota-derived agonists are decreased by DHA (Appendix Fig S2A; Lytle et al, 2015). This was an important observation considering that deficiency of either TLR2 or TLR4 attenuates severity of NADFL/NASH in mouse models (Spruss et al, 2009; Miura et al, 2013; Wu et al, 2020).

In contrast to TLRs, EGFR and other receptors that BTC binds are widely expressed across most cells in liver including macrophages. Moreover, EGFR deficiency in macrophages, but not in hepatocytes was shown to attenuate NASH in mice (Lanaya et al, 2014).

Thus, it is plausible that the impact of DHA on macrophage related to NAFLD/NASH pathogenesis can be explained by the fact that DHA simultaneously limits cell access to BTC and TLR2/4 agonists. Moreover, despite the limitations of our cell culture system and the difficulty of transition from mice to humans, we observed that BTC combined with TLR 2/4 agonists induce the integrin signaling pathway which was inhibited by DHA in the liver. Inspecting individual genes revealed that all detected integrins (ITGα1,6,9—Fig 5J, Appendix Fig S2I) required all three compounds for induction except for ITGβ1 that increased only with only BTC stimulation in THP-1 cells. Interestingly, ITGβ1 was also increased by BTC in hepatic stellate cells (Appendix Fig S1H). Of note, proteins coded to ITGαs and ITGβs form a complex that binds collagens (Bourgot et al, 2020). Accordingly, blocking integrin signals has been shown to attenuate fibrosis (Agarwal, 2014; Rahman et al, 2022). Hence, our results taken together with already existing literature about other pathologies (Bourgot et al, 2020) suggest that DHA inhibits the macrophage contribution to fibrosis by simultaneously inhibiting microbiota-derived signals and BTC. Furthermore, while BTC is a new player in this arena, the role of microbiota-derived signals in activating a profibrotic program in macrophages has been reported for different diseases in a few organs (He et al, 2021; Costa et al, 2022).

Leveraging complex multi-omic and single cell data (Xiong et al, 2019) and using a systems approach, our study constructed a statistical network model of cell–cell interactions affected by DHA (Fig EV4A). More importantly, network cause–effect related information flow (degree and BiBC) combined with a ligand–receptor database (Abugessaisa et al, 2021; Fig EV4A) allowed us to infer that a candidate master regulator molecule (BTC) produced by one cell (cholangiocytes) acts upon several other cells. DHA may affect several cell types involved in liver fibrosis. Our study, however, reveals that BTC inhibition by DHA simultaneously disrupts the integrin pathway in macrophages and TGFβ-2-driven collagen production by hepatic stellate cells, two processes that synergize in the development of liver fibrosis. Thus, we propose that removal of BTC and Tlr2/4 agonists prevents binding of integrins to collagen that is required for the scar development.

The main outcomes of our study, however, might be missing additional beneficial effects of DHA, especially those that are not related to BTC reduction and fibrosis. This is because we focused our investigations on cellular/molecular events relevant to NASH and its progression to liver cancer in humans. Specifically, while our network analysis modeled NASH in mice, we used it as a first step in identification of the most critical causal pathways upregulated in hepatic cancer meta-analysis and reversed by DHA. In this analysis we found ERBB as a top pathway with BTC being the top-ranked molecule in this pathway altered by DHA. Accordingly, our in vitro investigations of human hepatic stellate cells demonstrated that BTC promotes cell growth, which was an expected finding considering that BTC belongs to a family of growth factors (Olayioye et al, 2000; Wieduwilt & Moasser, 2008; Chen et al, 2016). Thus, in addition to promoting fibrosis by upregulation of TGFB2 in stellate cells, it may also increase numbers of these collagen-producing cells in the liver. Furthermore, our gene expression analyses of effects promoted by BTC in vitro and downregulated by DHA in the liver also support this growth function for BTC (Fig 4A), pointing to an activation of ERBB pathways, cell growth and even the epithelial mesenchymal transition. Dysregulation of these pathways are hallmarks of cancer (Dahlhoff et al, 2014; Lanaya et al, 2014; Chava et al, 2022). The growth-promoting actions of BTC are known to be mediated by epidermal growth factor receptors (ERBBs), namely ERBB1 (EGFR; Chava et al, 2022), ERBB2, ERBB3, and ERBB4. In liver however, the mechanism for BTC dependent cell proliferation has not been elucidated. The BTC-EGFR-ERBB4 pathway in pancreatic ductal adenocarcinoma has been well established by several groups and a BTC knock out can partially rescue the cancer progression (Hedegger et al, 2020). In the in vivo NASH model, we have a significant reversal of ERBB4 expression by DHA (Appendix Fig S2B). Furthermore, DHA also inhibits GRB2 and GSK3β, genes which are downstream in the ERBB signaling pathway found in our liver cancer meta-analysis (Fig 3A and C) and known to promote multiple types of cancers (He et al, 2021). These results suggest that DHA, besides its main effect as BTC reduction, might also inhibit residual activation of ERBB pathway via BTC-unrelated mechanisms. In addition, 49 EMT related genes, including the TGFβ-2 mediated (Dropmann et al, 2016) SMAD-EMT-cancer pathway were reversed by DHA treatments (Fig 3A). Therefore, our results suggest that the DHA attenuation of BTC-TGFβ-2-dependent molecular events might not be limited to reversal of fibrosis in NASH (Lytle et al, 2015) but also has a promise in preventing NASH's progression into liver cancer.

While investigating regulatory events upstream of BTC, we found FOXO3, a transcription factor with well-established function of tumor suppression (Liu et al, 2018; Tsuji et al, 2021), as a probable mediator of inhibitory effects of omega 3 fatty acids on BTC. These results provide additional details for understanding the molecular mechanism omega 3 action in the liver and contribute to rationale to use them for liver cancer prevention.

Clinical trials using EPA and/or DHA in NAFLD/NASH therapy have shown mixed, but promising results with significant improvement with disease severity, including hepatosteatosis and liver injury (Hodson et al, 2017; Zöhrer et al, 2017; Musa-Veloso et al, 2018; Okada et al, 2018; Tobin et al, 2018; Šmíd et al, 2022). Unfortunately, these studies did not address whether patients that did not respond to therapy represent a cellular/molecular subtype of NAFLD that is not treatable only by ω3 PUFA or perhaps have such severe disease that they are unlikely to respond to any therapeutic agent.

Despite the novelty and importance of our findings several questions remain to be investigated. For example, to evaluate the relative contribution of down regulation of BTC by DHA for treatment of NASH and prevention of cancer, experiments using mouse models are warranted.

In conclusion, with the discovery of BTC as a candidate to be one of the key mediators of ω3 PUFA therapeutic effects, our study opens a new avenue for investigation of NAFLD/NASH. In addition to finding new mechanisms of action of DHA, this study is the first to demonstrate that BTC can induce TGFβ-2 and synergize with microbial signals in the induction of integrins. Thus, while few earlier studies (Moon et al, 2006) showed increase of BTC in the liver tumors, our robust meta-analysis coupled with evidence for causal contributions shed a new light to this molecule in the pathogenesis of this cancer. Moreover, BTC's role in human NAFLD/NASH is entirely uncharted territory. Therefore, future studies should investigate if BTC-triggered gene expression signatures can serve as biomarkers guiding personalized ω3 PUFA therapy, as targets of new NAFLD/NASH drugs, and finally as a predictors of hepatic cancer risk in humans.

Materials and Methods Animals and diets

Study design for DHA mediated NASH prevention and remission in Male Ldlr^−/− mice.

This study was carried out in strict accordance with the recommendations in the Guide for Care and Use of Laboratory animals of the National Institutes of Health. All procedures for the use and care of animals for laboratory research were approved by the Institutional Animal Care and Use Committee at Oregon state University (Permit number: A3229-01). Anthropometric, plasma and liver samples used in this study were obtained from our previously published NASH prevention and treatment studies (Depner et al, 2013a; Lytle et al, 2015). Briefly, male mice (B6:129S7-Ldlr^tmHer/j, stock# 002207) were purchased from Jackson Labs and were group housed (4 mice/cage; n = 8 mice per group) and maintained on a 12-h light/dark cycle. Mice were acclimatized to the Oregon State University Linus Pauling science center vivarium for 2 weeks before proceeding with the experiments.

NASH prevention study

This study was designed to determine if EPA and DHA differed in their capacity to prevent western diet-induced NASH (Depner et al, 2013a). Mice consumed the one of the following 5 diets, ad libitum for 16 weeks; each group consisted of 8 male mice. Purina chow 5001 consisting of 13.5% energy as fat and 58.0% energy as carbohydrates was used as the Reference diet (RD). The western diet (WD) was obtained from Research Diets (12709B) and used to induce NASH. The WD consisted of 17% energy as protein, 43% energy as carbohydrate, and 41% energy as fat; cholesterol was at 0.2% wt:wt and does not contain either EPA or DHA (Depner et al, 2013a). The WD was supplemented with olive oil (WDO), eicosapentaenoic acid (WD + EPA), docosahexaenoic acid (WD + DHA), or both EPA and EHA (WD + EPA + DHA). EPA and DHA were added to the diets to yield 2% of total calories; for the EPA + DHA, each was added to yield 1% of total calories, i.e. 2% total calories as C_20-22 ω3 PUFA. Olive oil was added to the WD to have a uniform level of fat energy in all the WDs. Preliminary studies established that the addition of Olive oil had no effect on diet-induced fatty liver disease in Ldlr^−/− mice. EPA was purchased from Futurebiotics as Newharvest EPA, a DuPont product, while DHA was obtained as a generous gift from DSM, formally Martek Bioscience). The amount of EPA or DHA added to the diets is equivalent to the amount prescribed for treating patients for hypertriglyceridemia (Harris et al, 1997; Davidson et al, 2007). At the end of the 16-week feeding trial, mice were euthanized with CO₂, blood (RBC and plasma) and liver were collected and stored at −80°C for later extractions of RNA, lipids, proteins.

NASH treatment study

This study was designed to assess the capacity of DHA to reverse the effects of WD-induced NASH (Lytle et al, 2017). As such, male Ldlr^−/− mice were fed the WD for 22 weeks. These mice were separated into two groups: one group was maintained on a WD + olive oil, with the other group was maintained on the WD + DHA. The diet composition was as described above for the prevention study. Both groups were euthanized after 8 weeks of these diets. A control group was maintained on the RD for 30 weeks. At the end of the study, mice were euthanized, and samples collected and processed as above.

Liver histology

Approximately 100 mg of fresh mouse liver from each animal was fixed in buffered formalin, paraffin embedded, sliced and stained with hematoxylin–eosin (H & E), trichrome or Picro Sirius red (PCR; Nationwide Histology, Veradale, WA). Each slide contained 2–4 liver slices. Histological analysis and scoring for microsteatosis and macrosteatosis, inflammation (leukocytes) and fibrosis were provided by two investigators using the modified Kleiner scoring system established for mouse models of NAFLD as described previously (Garcia-Jaramillo et al, 2019). Histological samples were blinded as to diet, timepoint and diet. Digital images were taken with a Nikon Eclipse 6 microscope and digital camera (Mpixel) and NIS-BR Elements imaging software (v21.1; www.nikonmetrology.com). Digital images taken at 400× magnification were used for steatosis scoring. An effort was made to place the central vein of a lobule at one of the corners of the image so that each image covered at least one-quarter of that lobule, including all three zones of the hepatic lobule. Steatosis was objectively analyzed as the average percent surface area occupied by vacuoles using the image analysis software, ImageJ (https://imagej.nih.gov/ij/). Two images were taken of H&E-stained sections at 400× and the percent of affected surface area was calculated for each. The two values were then averaged. Steatosis was subjectively analyzed as the percent of affected surface area observed on H&E-stained slides at 100× (10× objective) and 400× magnifications. Vacuolization was characterized as macro vesicular, in which vacuoles displace hepatocyte nuclei, or micro vesicular. Macro vesicular and micro vesicular steatosis were scored separately. Severity was scored using the scale: 0 (0%), 1 (> 5% but < 33%), 2 (> 33% but < 66%), 3 (> 66%). When possible, the distribution of vacuoles was described as pericentral, midzonal, or periportal. Inflammation was defined as intralobular inflammatory foci of at least 5 leukocytes associated with disruption of hepatic plates or increased hepatocellular eosinophilia. Inflammation scores were calculated as the total number of clusters averaged over 5 fields in H&E-stained tissues examined at 100× (total number of clusters in 3.1 mm²). The following scale was used: normal 0 (< 0.5 foci), slight 1 (< 0.5–1 foci), moderate 2 (1–2 foci), severe 3 (> 2 foci). Fibrosis was also objectively quantified as percent surface area occupied by Sirius Red-stained collagen by image analysis using ImageJ. Two images were taken at 100× from the liver section of each mouse and the calculated percentage areas were averaged. Fibrosis was subjectively analyzed to determine severity and distribution patterns, perisinusoidal, periportal, pericentral, or bridging. The following scale was used: absent (0), mild (1), moderate (2), or severe (3).

Cell culture

LX2 cells were obtained from SL Friedman (Mount Sinai Medical School; Xu et al, 2005). LX2 cells are activated human hepatic stellate cells; they were maintained in DMEM with 10% FCS containing penicillin and streptomycin. Cells were plated onto 100 mm plastic petri dishes ∼ 100,000 cells/plate and treated with fatty acids (at 50 μM) in endotoxin-free BSA (at 20 μM) during the growth phase. Fatty acids [oleic acid (18:1 ω9) and DHA (22:6, ω3); Nu-Chek Prep] were used at 50 μM for 2–48-h cycle treatments. After this pre-treatment, cells were trypsinized, washed in PBS, counted and plated at 3,000 cells/well in a 96-well plate. Cells were fed DMEM with 1% FCS without or with betacellulin (human recombinant betacellulin [BTC], R&D Systems) for 72 h; the concentration ranged from 1.25 to 25 ng/ml. At the completion of BTC treatment, media was removed, washed with PBS and DNA/well was quantified using the CyQuant cell proliferation assay (Thermofisher); DNA fluoresces was quantified (excitation at 485 nm & emission at 530 nm). This experiment was repeated 3 times.

Collagen production

Pico Sirius red quantitation of collagen production: Cells were pretreated with fatty acids as described above, trypsinized, counted and plated onto 12-well cell culture plates at 80,000 cells/well in DMEM +1% FCS and treated without and with 20 μg/ml BTC for 72 h. At the end of treatment, media was removed, cells were washed with PBS and stained with pico Sirius Red (Abcam) for 2 h at room temperature. After staining cells were washed with an acetic acid solution, photographed and the stain was removed using 50 mM NaOH. The level of staining per well was quantified at 540 nm. The level of protein/well was quantified after solubilizing the protein in 50 mM NaOH and using the Pierce BSA kit. This experiment was repeated 4 times.

Cell viability

Alamar Blue (Creative Labs) assessment of cell vitality (NADH conversion to NAD⁺): Using the protocol described above, we assessed the vitality of cells after treatment with fatty acids and BTC. Alamar blue (50 μl/1.0 ml media) was added to the cells and fluorescence was measured at 590 nm 4 h after Alamar Blue addition.

RNASeq

RNA was extracted from LX2 cells using Trizol (Invitrogen) as previously described (Lytle et al, 2015) from cells treated with fatty acids (oleic acid and DHA at 50 μM for 72 h) as described above. Cells were seeded onto 6-well plates at 40,000 cells/well; and cells were treated with fatty acids for 96 h as described above. Afterward, media was removed, and cells were treated without and with BTC at 20 ng/ml for 72 h. Cells were harvested for RNA extraction (Trizol). cDNA was prepared for RNASeq analysis as described.

THP-1 cells BTC and TLR co-stimulation

THP-1 monocytes were cultured in RPMI 1640 medium adjusted to contain 4.5 g/l glucose and supplemented with 1% Penicillin/Streptomycin, 10% FBS, 1 mM sodium pyruvate, 10 mM HEPES. For experiments, monocytes were first seeded in 24 well plates at 400,000 cells/well in 1 ml of complete medium with 50 ng/ml PMA (phorbol 12-myristate 13-acetate) to induce polarization, for 24 h. Then, attached cells were washed with sterile PBS to remove residual serum and PMA containing media. Cells were stimulated with a TLR2 agonist at 40 ng/ml (PGN-BS, Invivogen, San Diego, CA), TLR4 agonist at 4 ng/ml (LPS-B5, Invivogen, San Diego, CA), BTC at 40 ng/ml (human recombinant Betacellulin protein, R&D Systems, MN), or combinations of all for 6 h. Treatments were prepared in serum-free semi-complete RPMI 1640 media (see above for other components) supplemented with 20 μM BSA and 50 μM BHT (Butylated hydroxytoluene). Cells were lysed with 300 μl RLT buffer (Qiagen) and cell lysates were stored at −80°C.

Primers qRT–PCR data analysis

THP-1 cells' response to TLR and BTC stimulation was assessed by qRT–PCR (Appendix Table S1). Briefly, raw Cycle Threshold (CT) values from the StepOnePlus Real Time PCR instrument for genes of interest were normalized to C_T values of a housekeeping gene, TMEM59, by delta C_T method and relativized by $2^{-\mathrm{\Delta }{C}_{\mathrm{T}}}$. Data were then median normalized and Log₂ transformed before being plotted in GraphPad Prism 9.

RNA extraction and RNA sequencing library preparation

RNA was extracted from cell lysates using the RNeasy Mini Kit and subjected to a DNase treatment according to manufacturer's protocols (Qiagen) then stored at −80°C until further use. mRNA libraries were prepared for sequencing with the Lexogen QuantSeq 3′ mRNA-Seq Library Prep Kit (FWD) HT for Illumina Sequencing platforms (Kit#k15.384) and sequenced on the Illumina NextSeq at Oregon State University.

cDNA and qRT–PCR

cDNA was prepared from 0.5 to 1 μg of RNA via reverse transcription using the qScript XLT cDNA SuperMix kit (Quantabio). qRT–PCR was performed for gene expression using the AzuraView GreenFast qPCR Blue Mix HR (Azura Genomics). 96-well plates were prepared with 10 ng of cDNA in triplicate reactions for each gene and sample and run on an Applied Biosystems StepOnePlus Real Time PCR instrument.

Sequencing of RNA (RNAseq)

RNA libraries were prepared with the QuantSeq 3′mRNA-Seq Library Prep Kit (Lexogen) for the Apollo 324 NGS Library Prep System and sequenced using Illumina NextSeq. The sequences were processed to remove the adapter, polyA and low-quality bases by BBTools (https://jgi.doe.gov/data-and-tools/bbtools/) using bbduk parameters of k = 13, ktrim = r, forcetrimleft = 12, useshortkmers = t, mink = 5, qtrim = r, trimq = 15, minlength = 20. Then the reads were aligned to mouse genome and transcriptome (ENSEMBL NCBIM37) using Tophat (v2.1.1) with default parameters. The number of reads per million for mouse liver genes were counted using HTSeq (v 0.6.0) and quantile normalized. Cell lines sequencing/analysis was performed similarly, but with the bbduk parameters of k = 13, ktrim = r, forcetrimleft = 11, useshortkmers = t, mink = 5, qtrim = r, trimq = 10, minlength = 20. The reads were aligned to the human genome and transcriptome (Gencode v40) using STAR v2.5.3a. BRB-ArrayTools was used to identify differentially expressed genes between treatments.

Metabolomes and lipidomes

Data were prepared, and analysis was carried out from the NASH Prevention study and Treatment study. Hepatic lipids and non-lipid metabolites were extracted and subject to UPLC/MS/MS analysis as previously described (Garcia-Jaramillo et al, 2019; Rodrigues et al, 2021) with minor modifications. Processed normalized data are available in Dataset EV8.

Treatment categorization of omics data

The genes and other parameters whose expressions or values are significantly changed by WD (FDR < 10%) were considered for treatment effects. Out of these parameters and genes, those that have treatment effect reversal with a P-value < 0.05 were categorized as DHA if uniquely DHA (opposite to WD + O/ND & WD + DHA/WD + O P-value < 0.05), uniquely EPA (opposite to WD + O/ND & WD + EPA/WD + O P-value < 0.05), similar in both EPA&DHA (opposite to WD + O/ND & both WD + EPA/WD + O and WD + DHA/WD + O, P-value < 0.05) and with no treatment effect (NA; with both WD + EPA/WD + O and WD + DHA/WD + O, P-value > 0.05). To be consistent with many previous studies, the treatment effects were also tested in the combination of preventive treatment though not used further in the analysis (EPA + DHA; WD + EPA + DHA/WD + O, P-value < 0.05).

Reconstructing the NASH liver multi-omic network

The network reconstruction was carried out as described in the previous papers from our group with minor dataset specific modifications (Dong et al, 2015, Li et al, 2022). The genes and other parameters whose expressions or values are significantly changed by WD (FDR < 10%) were chosen for constructing the NASH network. First, from liver tissue between all pairs of genes (GE) and metabolic parameters (phenotypes, P) spearman rank correlations were calculated by pooling the samples per diet (WD + O, WD + EPA and WD + DHA). Meta-analysis was performed to retain edges with same sign of correlation coefficient in all three diets. Edges were further filtered by the following criteria: individual P-value of correlation within each diet from pooled experiments < 20%, combined Fisher's P-value over diets from pooled experiments < 5% and FDR cutoff of 10% for edges within tissues and for phenotypes and between lipidomic, metabolomic and plasma biochemicals and edges needed to satisfy principles of causality (i.e., satisfied fold change relationship between the two partners in the WD + O vs. ND comparison). Next, correlations were calculated per diet for the experiment pairwise between parameters (Gene Expression + Lipidomic and Metabolomic data [LM]) and (P + LM). Finally, edges obtained from pooling were retained if they had the same sign of correlation coefficient as in 3 groups (3 diets, WD + O, WD + EPA and WD + DHA). False positive edges were removed (Yambartsev et al, 2016) while pooling the different diets in the creation of the network. The proportion of genes, metabolites and lipids that made it to the final network (following statistical cutoffs) was determined after applying selection for significance in differentially expressed parameters in liver prior to applying correlation cutoffs.

Single cell RNA (scRNA) sequencing data analysis Datasets

Single cell dataset (GSE129516) for mouse NASH model was obtained from single cell RNA-sequence on non-parenchymal cells of healthy vs. NASH mouse liver. These are then reanalyzed and used in our multi-omic network analysis to infer liver cell types. Human liver single cell RNA sequencing data (GSE136103) from normal and Cirrhosis patients were reanalyzed (Ramachandran et al, 2019; Xiong et al, 2019).

scRNA sequencing analysis

The raw gene expression matrix (UMI counts per gene per cell in the liver tissue) was filtered, normalized, and clustered using a standard Seurat version 3.1.0 in R [https://www.R-project.org/] (Stuart et al, 2019). Cell and gene filtering were performed as mentioned in previous publications. During quality filtering, cells with a very small library size (< 5,000) and a very high (> 12%) mitochondrial genome transcripts were removed. The genes detected (UMI count > 0) in less than three cells were also removed from further processes. Then log normalization and further clustering is performed using standard Seurat package procedures. Principal component analysis was used to reduce dimensions that represent cell clusters. The number of components from this analysis used for the elbow of a scree plot, which further aid in selecting the significant clusters. The different cell type clusters in a sample were visualized using t-distributed Stochastic Neighbor Embedding of the principal components as implemented in Seurat. The liver tissue specific cell-type identities for each cluster were determined manually using a compiled panel of available immune and other cell specific marker expression as per the previously published papers (Ramachandran et al, 2019; Xiong et al, 2019).

Single cell RNA sequence for assignment of gene in the NASH network to a specific cell type

The normalized UMI > 1.0, with a fold change significantly and uniquely expressed genes in a cell specific cluster were then assigned to that specific set of genes in the network to indicate they belong to that specific cell type. It is the primary rule to assign a gene to a cell type. Next, higher expression in the cell cluster (and optional fold change (log₂FC > 0.25, P value < 0.05) for a gene gets assigned with that specific cell type of the tissue. Here, additionally, ranking by average expression for each gene in the clusters helps to determine its cluster specificity by the higher expression in that cell type than another in the whole tissue. This is implemented for evaluating the highest cell cluster average expression of a gene, among all other cell clusters in network.

Detecting subnetworks and functional enrichment

Infomap (mapequation.org) was used to identify subnetworks using the default commands. Functional enrichment of clusters was then performed using metascape (http://metascape.org; Zhou et al, 2019).

Additionally single cell data overlay on NASH network as mentioned above allowed the cell type specific gene sub clusters.

Identifying of key nodes between subnetworks using BiBC analysis

Analysis of networks was performed using the python package NetworkX v2.2. Bipartite betweenness centrality (BiBC) was calculated between all cell cluster pairs (66 in total) of the 12 clusters previously identified based on single cell data overlay on NASH network. The nodes were then ranked by their resulting BiBC and scaled to range of 0–100. BiBC was also calculated and scaled similarly pairwise between all genes and anthropometric data, between all genes and Lipidomic/Metabolomic data.

Creation and analysis of random networks

Random networks were created similar to as was described in Kahalehili et al (2020). Briefly, 5,000 Erdos-Renyi random networks were created, utilizing the same number of nodes and edges that were present in the real, reconstructed network. BiBC was calculated both (i) between DHA controlled genes and anthropometric nodes and (ii) between DHA controlled genes and DHA controlled metabolomic/lipidomic data. Plotly (https://chart-studio.plotly.com/create/#/) was used to create the 2D contour histograms for visualization of the random network results.

Intrahepatic ligand–receptor interaction network

The knowledgebase of ligand–receptor interaction information available for mouse genes (Abugessaisa et al, 2021), was overlayed on to reconstructed NASH Multiomics network genes. Additionally, this allowed the interrogation of ligand–receptor network with respect to each cell in the network that a gene can be represented as ligand or receptor.

Human liver cancer meta-analysis

Human cancer data from hepatic cellular cancer (HCC) and cholangiocytes cancer (CC) were selected from GEO data sets (GSE14520, GSE26566, GSE102079, GSE56140, GSE98617, GSE76427, GSE84005). Meta-analysis using RankProd method described in the publicly available tool OMiCC (Shah et al, 2016) was carried out for these 7 sets of data with 32 healthy, 260 non-tumor samples and 544 tumor samples (Appendix Fig S3A). At first sample sets of both tumor types were compared against their respective paired non-tumor or healthy samples. Additionally, a standard (Fisher) approach of meta-analysis was also carried out and an overlapping set of genes were selected.

Then a signature set of genes were identified by matching (among orthologous) genes between mouse and human using Biomart (Cunningham et al, 2022) for the concordant fold change direction as WD + O/RD in the NASH network model and genes that are significant with FDR < 15% and Fisher P value < 5%. The human liver cancer meta-analysis signature genes overlapped with the NASH mouse model network genes with treatment effects to identify subset of genes.

Gene ontology analyses

The gene ontology and functional enrichment were carried out using Metascape and innateDB (Breuer et al, 2013; Zhou et al, 2019) with mouse or human reference databases. The reference database is determined depending on the species of the specific data in question.

Pathway summary metric

Using published approach (Levine et al, 2006) with minor adaptions described in our previous papers (Shulzhenko et al, 2011; Koscso et al, 2020), each top gene ontology pathway enriched in the THP-1 experiment (with BTC and TLR ligands and significantly revered by DHA in NASH preventive model) was recognized and the genes belonging to the pathways were identified. Using median normalization, a value for each gene in each treatment replicates and then additionally using their median, summary metric was calculated for the individual pathway. A paired t-test between normal and treatment condition was performed to evaluate significance.

Upstream regulation and transcription factor analysis of BTC

The overall analysis strategy consisted of: (i) establishing functional effects of BTC represented by BTC-dependent gene expression profile. (ii) Using network analysis among genes regulated by omega-3 we infer candidate regulators of BTC-dependent profile expressed in cholangiocytes. (iii) searching for transcription factors among group of genes found is step 2; (iv) identifying which TF has a binding motif in BTC gene.

Specifically, identification of top BiBC involved in the regulation of Betacellulin in the cholangiocytes were carried out by interrogating the network nodes with the Omega 3 reversal effect in NASH network derived from in vivo data and assigned to cholangiocytes while opposite in BTC treatment in vitro in both LX2 and THP1 experiments to the lipids with EPA/DHA acyl chain in them (Appendix Fig S6A and B). Then from this list of BiBC nodes, the genes with higher than 50 edges/degrees and top 1% < BiBC were further analyzed for transcription factor (TF) activity. Using the TF binding prediction analysis tools (SCENIC, TFBSPred; Aibar et al, 2017; Zogopoulos et al, 2021), the TF motif and binding site in the promoter of Btc were predicted and verified.

Statistical analysis

In mouse studies, data are expressed as geometric means of replicates. Data are shown as the mean ± standard deviation in animal studies and in vitro experiments. Group comparisons were performed using an unpaired t test and ordinary one-way analysis of variance (ANOVA), followed by Tukey's post hoc or Dunnett's multiple comparisons tests, where P < 0.05 indicates statistical significance. In cell line and RNA sequence analysis, comparisons between groups were performed using Student's t test or the Mann–Whitney U test and Kruskal-Wallis test when appropriate. Categorical variables are shown as counts and percentages. Differences between categorical variables were assessed with the chi-squared test or Fisher's exact test. Spearman's rank correlation rho coefficients were calculated for network edges between all parameters using R statistical packages. Details of statistical analyses are described additionally in the corresponding figure legends. GraphPad Prism 9 was used for all analyses.

Data availability

RNA-Seq data: Gene Expression Omnibus GSE215223 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE215223).

RNA-Seq data: Gene Expression Omnibus GSE215224 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE215224).

RNA-Seq data: Gene Expression Omnibus GSE215225 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE215225).

RNA-Seq data: Gene Expression Omnibus GSE215227 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE215227).

Network Model: NASH_DHA_Betacellulin_Multi-omics_Network (https://www.ndexbio.org/#/network/60fc7ad9-3d93-11ee-aa50-005056ae23aa?accesskey=764ecdb2a2d2a6178586a7afd0755c9c8fab9e9a531b4336cf26b46432cac4a1).

Acknowledgements

We would like to thank Dr. Scott Friedman from the Division of Liver Diseases at the Icahn School of Medicine at Mount Sinai for providing the LX2 cell line and Christiane V. Löhr, MedVet, Dr. med. vet., PhD for the evaluation of liver histology. And the personnel of CQLS at the Oregon State University for IT support. Funding: US Department of Agriculture, National Institute of food and agriculture grant 2009-65200-05846 (DBJ), National Institutes of Health grants, R01 DK094600 (DBJ), R01 DK112360 (DBJ) and R01 DK103761 (NS).

Author contributions

Jyothi Padiadpu: Data curation; software; formal analysis; supervision; validation; investigation; visualization; methodology; writing – original draft; writing – review and editing. Manuel Garcia-Jaramillo: Data curation; formal analysis; investigation; methodology. Nolan K Newman: Data curation; software; formal analysis; validation; investigation; visualization; methodology; writing – original draft; writing – review and editing. Jacob W Pederson: Validation; investigation; visualization; methodology. Richard Rodrigues: Software; formal analysis; investigation; methodology. Zhipeng Li: Formal analysis; investigation; methodology. Sehajvir Singh: Formal analysis; investigation. Philip Monnier: Formal analysis; investigation; methodology. Giorgio Trinchieri: Resources; investigation; methodology. Kevin Brown: Resources; software; methodology. Amiran K Dzutsev: Resources; investigation; methodology. Natalia Shulzhenko: Resources; supervision; funding acquisition; investigation; visualization; methodology; project administration; writing – review and editing. Donald B Jump: Conceptualization; resources; formal analysis; supervision; funding acquisition; investigation; methodology; writing – original draft; project administration; writing – review and editing. Andrey Morgun: Conceptualization; resources; software; formal analysis; supervision; visualization; methodology; writing – original draft; project administration; writing – review and editing.

Disclosure and competing interests statement

The authors declare that they have no conflict of interest.