Cheng Xiao
Siliang Chen
Chunru Yang1
Jieying Liu
Miao Yu- 1Department of Endocrinology, Key Laboratory of Endocrinology, National Health Commission, Peking Union Medical College Hospital, Peking Union Medical College, Chinese Academy of Medical Sciences, Beijing, China
- 2Department of Dermatology and Venereology, West China Hospital, Sichuan University, Chengdu, China
- 3Laboratory of Dermatology, Clinical Institute of Inflammation and Immunology, Frontiers Science Centre for Disease-related Molecular Network, West China Hospital, Sichuan University, Chengdu, China
- 4Department of Medical Research Center, Peking Union Medical College Hospital, Peking Union Medical College, Chinese Academy of Medical Sciences, Beijing, China
Polyunsaturated fatty acids (PUFAs) play important roles in the aetiology and pathogenesis of metabolic dysfunction-associated fatty liver disease (MAFLD). However, the underlying molecular mechanisms are not understood. We analysed a public GEO dataset, GSE89632, to identify differentially expressed genes (DEGs) in MAFLD. Weighted gene coexpression network analysis (WGCNA) was used to reveal the core gene regulation network and to explore the PUFA-related hub genes in MAFLD. We experimentally verified these genes by quantitative reverse transcription PCR in high-fat diet (HFD)-fed mice. A total of 286 common DEGs (89 upregulated; 197 downregulated), mostly related to inflammatory and immune responses, were identified. Six modules were constructed using WGCNA, and 2 modules showed significant correlations with PUFAs. After combining these 2 modules with DEGs, the top 10 hub genes were identified. We further established a MAFLD mouse model with liver steatosis, as proved by HE and Oil Red O staining. Of the hub genes, ADAM metallopeptidase with thrombospondin type 1 motif 1 (adamts1) (p = 0.005) and transforming growth factor β3 (tgfβ3) (p < 0.001) showed significantly lower mRNA expression in MAFLD in vivo. adamts1 and tgfβ3 bridged PUFAs and MAFLD, which might be potential causative genes and therapeutic targets of MAFLD.
Introduction
Metabolic dysfunction-associated fatty liver disease (MAFLD), which is also referred to as nonalcoholic fatty liver disease (NAFLD) (Eslam et al., 2020), includes a series of chronic liver diseases, ranging from steatosis (SS) to nonalcoholic steatohepatitis (NASH), cirrhosis, and hepatocellular carcinoma. Previous studies reported that MAFLD was closely associated with an increased risk of developing cardiovascular diseases, type 2 diabetes, and other diseases (Anstee et al., 2013; Park et al., 2013; Adams et al., 2017). Currently, one-quarter of the world’s population suffers from MAFLD (Younossi et al., 2019), which is a major health problem worldwide. However, the molecular aetiology and appropriate pharmacotherapeutic approaches for MALFD have not yet been elucidated (Lim et al., 2021). Therefore, there is an urgent need to elucidate the molecular mechanism underlying the pathogenesis of MALFD.
A “two-hit hypothesis” was proposed to explain the pathogenesis of MAFLD (Day and James, 1998; Dowman et al., 2010). According to this theory, lipid accumulation or SS in the liver in the form of triglyceride (TG) and free fatty acids are considered the first hit, rendering the liver more vulnerable to inflammatory insult, which acts as the second hit. Cytokines and chemokines (Peverill et al., 2014) secreted by hepatocytes or activated neutrophils further enhance injury to hepatic tissue, which leads to MAFLD (Machado and Diehl, 2016; Koyama and Brenner, 2017; Schuster et al., 2018). Multiple studies have suggested that polyunsaturated fatty acids (PUFAs) play important roles in the aetiology of MAFLD. It has been reported that dietary PUFA intake is lower in individuals with MAFLD (Puri et al., 2007; Da Silva et al., 2014; Arendt et al., 2015), and PUFA supplementation may ameliorate the clinical symptoms of NASH(Li et al., 2015). One of the explanations for the relationship between PUFAs and MAFLD is that PUFAs contribute to the improvement in lipid metabolism, increase in insulin sensitivity, and amelioration of inflammation (Di Minno et al., 2012; Imamura et al., 2016; Schulze et al., 2020), which further lead to a lower risk of MAFLD. However, the underlying mechanism needs to be further elucidated.
To identify crucial PUFAs-related genes in MAFLD, we used a public GEO dataset to analyse differentially expressed genes (DEGs) between MAFLD patients and healthy controls (HCs) and to determine their biological functions. Now many bioinformatics methods have been developed to identify potential therapeutic targets (Wang et al., 2022) (LIANG et al., 2022). Weighted gene coexpression network analysis (WGCNA) is a tool used to identify modules as candidate regulators and drivers of disease states (Langfelder and Horvath, 2008), and is widely used in liver diseases including porto-sinusoidal vascular disease, cholangiocarcinoma and hepatocellular carcinoma (Hernández-Gea et al., 2021; Long et al., 2021; XU et al., 2022). Besides, a recent study applied WGCNA to identify two hub genes NDUFA9 and UQCRQ which may be involved in the pathogenesis of MAFLD (Zeng et al., 2021). In this study, we applied WGCNA to construct gene modules and to identify hub genes in the PUFAs-related modules. The bioinformatics results were validated by comparing the expression of these hub genes between mice with MAFLD and control mice using quantitative reverse transcription PCR (RT–qPCR).
Materials and methods
Data preprocessing
Figure 1 shows the flowchart of the study design. The mRNA microarray dataset GSE89632 from the GEO database was used in our study. We downloaded the series matrix file and GPL14951 platform data table. We sorted the clinical information from the series matrix file, annotated the probe IDs in this matrix and generated an expression matrix with gene symbols using the microarray platform data table. Samples from 24 HCs, 20 patients with simple SS, and 19 patients with NASH were used for further analysis.
FIGURE 1. Flowchart of study.
DEG screening
We used the “limma” R package to screen DEGs between HC and patients with SS/NASH. We set the thresholds to |log2 (fold-change)| >1 and adjusted p < 0.05 (Lin et al., 2019; Chan et al., 2020; Dobson et al., 2020; Xiu et al., 2020). Then, we combined the DEGs from HC-SS and HC-NASH. We identified the common DEGs of these 2 comparisons. To determine the biological function of these common DEGs, we conducted Gene Ontology (GO) and Kyoto Encyclopedia of Genes Genomes (KEGG) pathway analyses on these DEGs using the “clusterProfiler” R package. We set the cut-off value for the functional annotation analysis to adjusted p < 0.05.
Weighted correlation network construction
We used the “WGCNA” R package for the network analysis. To construct a network with biological significance, we filtered genes with the top 25% expression variance to construct the gene coexpression network. We first built a correlation matrix based on Pearson’s correlation and transformed this correlation matrix into an adjacency matrix with the formula amn = |cmn|β, where amn represents the element in the adjacency matrix, cmn represents the element in the correlation matrix and β represents the soft threshold. To detect gene modules from the coexpression network, we constructed a topological overlap measure (TOM) matrix. We performed average linkage hierarchical clustering and used a dynamic tree cut to detect gene modules using the TOM matrix. We set the minimal module size as 50 and merged similar modules with a threshold of 0.25.
DEG coexpression network construction and hub gene mining
After gene module detection, we correlated the gene modules with clinical traits using the eigengene of each module. We identified modules with the strongest correlation with the levels of at least one PUFA (i.e., arachidonic acid, eicosapentaenoic acid, and docosahexaenoic acid), and we considered these modules to be PUFAs-related gene modules. We also performed GO analysis on the genes in the PUFAs-related modules to better understand their biological function. To construct DEG coexpression networks for hub gene mining, we mapped common DEGs onto these PUFAs-related modules and constructed DEG coexpression networks. Genes with the top 10 connectivities were considered hub genes in each DEG coexpression network.
Experimental validation of hub genes
C57BL/6 female mice (HFK Bioscience Co., Ltd., Beijing, China, SCXK-2016-0006) were housed in a specific pathogen-free environment (12-h light/dark cycle) with ad libitum access to food and water. After 1 week of adaptation, 6-week-old mice with similar body weights were randomly assigned to 2 groups: the control diet group (CD, n = 6), which was fed a standard rodent diet (AIN-93G) (15.8% of the calories as fat); and the high-fat diet group (HFD, n = 6), which was fed an HFD (D12492) (60% of the calories as fat). High fat diet formulations were listed in Supplementary Table S5. The oral glucose tolerance test (OGTT, oral administration of 1.0 g/kg of glucose, after 12 h fasting) was performed at the end of the study. All the animal procedures were conducted in compliance with the National Institutes of Health Guide for the Care and Use of Laboratory Animals.
Pathological staining
Mice were euthanized after 13 weeks to evaluate the pathological changes in their livers. HE staining was conducted according to the protocol of the HE Staining Kit (G1120, Solarbio). Oil Red O staining was performed using Oil Red O kits (G1261, Solarbio). Masson’s staining and Sirius red staining were performed using a Sirius red staining kit (S8060, Solarbio) and Masson’s staining kit (G1340, Solarbio), respectively, according to the manufacturer’s protocol.
RT–qPCR
Total RNA was extracted from liver tissues using a commercially available kit (Vazyme, RC112). The RNA concentration and quality were determined according to the 260/280 nm ratio, which was measured using a NanoDrop spectrophotometer (ND-100, Thermo Scientific). Total cDNA was synthesized from isolated RNA samples by HiScript III RT SuperMix (Vazyme, R323-01). The mRNA expression of target genes was quantified by RT-qPCR using ChamQ SYBR qPCR Master Mix (Vazyme, Q331-02) on a Light Cycler 480 (Roche, Basel, Switzerland). The GAPDH housekeeping gene (Sangon Biotech, B661304-0001) was used as the internal control. Relative target gene expression was calculated by the 2−△△CT method. The sequences of the primers used in RT-qPCR are listed in Supplementary Table S4.
Western blot
Proteins from mouse liver samples were extracted using RIPA buffer (APPLYGEN, China) containing protease inhibitors (Solarbio, China) and PMSF (Solarbio, China) Samples containing equal amounts of protein were denatured and subjected to electrophoresis in 10% SDS-PAGE gels followed by transfer to PVDF membrane and probed with specific antibodies, including ADAMTS1 (1:500), TGFβ3 (1:500), and GAPDH (1:1000) (Proteintech, Inc.). Blots bands were visualized using the horseradish peroxidase conjugated secondary antibodies and chemiluminescent substrate.
Glucose and lipid metabolism parameters measurements
Liver tissues were homogenized in triglyceride assay buffer, and serum was tested directly. TG levels were measured using the Triglyceride Quantification Assay Kit (ab65336, Abcam) according to the manufacturer’s instructions. Fasting blood glucose levels were determined using an automated analyzer (ADVIA1800 Siemens, Germany). Fasting insulin level was measured using the Mouse Insulin ELISA (ALPCO, America) according to the manufacturer’s instructions.
Statistical analysis
Bioinformatical analyses, such as DEGs screening, WGCNA, and functional annotation, were conducted in R v3.6.2. DEG coexpression network construction and hub gene mining were performed in Cytoscape v3.7.0. The student’s t-test was used to compare data from two groups. A p value < 0.05 was considered statistically significant.
Results
Filtering of DEGs
A total of 409 genes were significantly differentially expressed between the HC and SS groups (Figure 2A), whereas 410 genes were differentially expressed between the HC and the NASH groups (Figure 2C). The heatmap and the volcano plot show the relative expression patterns of these DEGs (Figures 2B,D). A total of 286 DEGs were common in both comparisons, of which 89 were upregulated and 197 were downregulated (Supplementary Figure S1; Supplementary Table S1). The functional annotation analysis of these common DEGs is shown in Figure 3 and Supplementary Table S2. The key GO terms in biological processes were mainly related to macrophage activation (GO: 0042116), positive regulation of inflammatory response (GO: 0050729), and cytokine production involved in immune response (GO: 0002367) (Figures 3A,B). After KEGG analysis, DEGs were mapped to KEGG pathways (Figures 3C,D). The most affected pathways were the tumour necrosis factor (TNF) signalling pathway (hsa04668), interleukin (IL)-17 signalling pathway (hsa04657), and nuclear factor κ-B (NF-κB) signalling pathway (hsa04064). Generally, all of these common DEGs might be primarily involved in inflammatory and immune responses.
FIGURE 2. DEGs screening for HC-SS and HC-NASH. In heatmap, red color indicates up-regulation while blue color indicates down-regulation. In volcano plot, red dots represent up-regulated genes while green dots represent down-regulated genes. (A,B) Heatmap and volcano plot for HC-SS. (C,D) Heatmap and volcano plot for HC-NASH.
FIGURE 3. Functional annotation for DEGs. (A,B) Bar plot and dot plot of GO analysis for DEGs. (C,D) Bar plot and dot plot for KEGG pathway analysis for DEGs.
Construction of the gene coexpression network
Based on clustering analysis, 2 outliers (GSM2385767 and GSM2385782) were identified and excluded before analysis (Supplementary Figure S2). A cut-off of R2 = 0.85 was used to select the soft-threshold β, and β = 16 was selected for network construction (Figures 4A,B). On the log-log plot the distribution approximately follows a straight line (R2 = 0.86), which is referred to as approximately scale-free topology. The constructed network met the requirements of scale-free topology. (Figures 4C,D). After module detection, we identified 6 gene modules (Figure 4E), and the correlation of these modules with clinical traits is presented in a heatmap (Figure 4F). Among the 6 gene modules, the brown module had the highest positive correlation with PUFAs levels, while the yellow module had the strongest negative correlation with PUFAs levels (Supplementary Figure S3). Functional annotation analysis showed that the brown module was mainly associated with inflammatory and immune responses (Figures 5A,B; Supplementary Table S3), and the yellow module was mainly associated with lipid metabolism (Figures 5C,D; Supplementary Table S3).
FIGURE 4. Weighted correlation network construction. (A,B) The soft-threshold β was chosen to be 16. (C,D) Log-log plot of whole-network connectivity distribution. The x-axis shows the logarithm of the whole network connectivity, and the y-axis is the logarithm of the corresponding frequency distribution. The network met the requirements of scale-free topology. (E) A total of 6 module was detected. (F) The correlation heatmap for gene modules and traits.
FIGURE 5. Functional annotation for PUFAs-related modules. (A,B) Bar plot and dot plot of GO analysis for brown module. (C,D) Bar plot and dot plot for GO analysis for yellow module.
DEG coexpression network and mining of hub genes
To combine DEGs with PUFAs-related modules, we mapped common DEGs onto the brown and yellow gene modules and acquired 2 DEG coexpression networks (Figures 6A,B). In coexpression network, a hub is a node with several links with other nodes that greatly exceed the average. The top 10 connectivity genes were the genes with the greatest number of links and were considered hub genes (Figures 6C,D; Table 1). To verify the hub genes, we established an HFD-fed mouse model to investigate the mRNA expression of these genes.
FIGURE 6. DEG co-expression network and hub genes. (A,B) DEG co-expression networks for brown and yellow modules. (C,D) hub genes for brown and yellow modules.
TABLE 1. Crucial genes for the brown and the yellow module.
HFD-fed mice developed liver SS and metabolic dysfunction
The body weight and hepatic weight significantly increased in the HFD group compared with the CD group (Figures 7A,B). The liver tissues of mice were stained with HE to observe lipid accumulation in hepatocytes as vacuoles (Figure 7C). Oil Red O staining confirmed the presence of massive lipid deposition in the livers of HFD-fed mice but not in those of control mice (Figure 7C). Masson’s staining and Sirius red staining, which were used to detect fibrosis, revealed no significant differences between the groups (Figure 7C). The HFD group also showed significantly higher hepatic TG contents (p = 0.031, Figure 7D), although there was no significance in the serum TG levels (Figure 7E). In addition, OGTT and area under curve of OGTT further showed impaired glucose tolerance in HFD-fed mice (Figure 7F). Both levels of serum glucose and insulin were significantly higher in HFD-fed mice (Figures 7G,H), which meant the HFD-fed mice developed insulin resistance. Overall, all of the above results confirmed the development of MAFLD in HFD-fed mice without obvious fibrosis.
FIGURE 7. Evaluation of metabolic phenotype in the HFD-fed mice. (A,B) The body weight and hepatic weight between HFD and CD groups, respectively. n = 6. (C)Representative HE staining, Oil-red O staining, Masson staining, and Sirius red staining images of liver from each group. Black scale bar, 200 um. (D,E)The hepatic TG and serum TG between HFD and CD group. n = 5. (F) The oral glucose intolerance test (oral administration of 1.0 g/kg of glucose, after 12 h fasting). (G,H) Fasting serum glucose and insulin were measured. n = 5. Data are mean ± SD. *p < 0.05 versus CD group.
Validation of hub genes in HFD-fed mice
To verify whether the hub genes identified in the two DEG coexpression networks were variably expressed, mRNA was extracted from liver tissues, and gene expression was quantified by RT-qPCR. It has been shown that ADAM metallopeptidase with thrombospondin type 1 motif 1 (adamts1), which is a gene in the brown module, was expressed at significantly lower levels in HFD-fed mice, and these results matched expectations well (p = 0.005, Figure 8A). However, the expression of the socs3, epha2, and zc3h12a genes showed no significant group differences (all p > 0.05, Figure 8B). In the yellow module, significantly lower tgfβ3 expression was observed in the HFD-fed mice (p < 0.001, Figure 8C), although rmnd1 and gper expression did not differ significantly between the groups (all p > 0.05) (Figure 8D). The changes of TGFβ3 and ADAMTS1 at the protein levels were further confirmed by WB assays. HFD-fed mice showed lower expression of TGFβ3 and ADAMTS1 (Figure 8E).
FIGURE 8. (A,B) Expression of the hub genes in the brown modules. n = 5. (C,D) Expression of the hub genes in the yellow modules. Results represent the mean ± SD of five to six independent experiments. n = 5. (E) One representative western blot of three independent western blot experiments is shown. Data are mean ± SD.*p < 0.05 versus CD group.
Discussion
In the present study, 286 common DEGs, which were mainly involved in inflammation and immune responses, were identified between HCs and patients with MAFLD. Based on the WGCNA, we constructed 6 modules, of which one inflammation-related and one lipid metabolism-related module were most strongly associated with PUFAs. The genes with the top 10 connectivities within the 2 modules were mined as hub genes. For the first time, we confirmed that two hub genes, namely, adamts1 and tgfβ3, were significantly differentially expressed between healthy mice and HFD-fed mice with MAFLD. These findings suggest that adamts1 and tgfβ3 play important roles in the process of MAFLD and may serve as biomarkers for the disease.
In our study, functional annotation analysis showed that pathways related to the inflammatory response, especially the macrophage activation pathway, TNF-α signalling pathway, NF-κB signalling pathway, and IL-17 signalling pathway, might play important roles in MAFLD. Similarly, multiple studies have reported that inflammation is crucial for the process of MAFLD. According to previous studies, the secretion of cytokines and chemokines, including TNF-α, IL-6, and IL-1β, by hepatic macrophages was increased in fatty livers (De Taeye et al., 2007; Stienstra et al., 2010; Chawla et al., 2011). Additionally, NF-κB and TNF-α appear to form positive feedback loops and exacerbate the inflammation and injury of hepatocytes (Rolo et al., 2012). In addition, IL-17, which is produced by pro-inflammatory T helper 17 (Th17) cells, can promote liver inflammation and fibrosis by facilitating the production of IL-6, IL-1, and TNF-α by inflammatory cells and activating hepatic stellate cells to produce collagen type I (Meng et al., 2012; Harley et al., 2014). Conversely, depletion of hepatic macrophages can effectively mitigate the progression of diet-induced SS and hepatic insulin resistance in rodent experiments (Huang et al., 2010; Lanthier et al., 2010; Lanthier et al., 2011). These findings suggested that the inflammatory response was activated in individuals with MAFLD, while macrophages might play a mediator role.
Although the aetiology of MAFLD is far from elucidated, growing evidence suggests that decreased intake of PUFAs is one of the causes of the disease (Puri et al., 2007; Da Silva et al., 2014; Arendt et al., 2015; Li et al., 2015). Based on the WGCNA, our results suggested that PUFAs were associated with MAFLD by regulating the coexpression network of inflammation and lipid metabolism (Figure 3). Consistent with our findings, many studies reported that PUFAs could regulate the TNF-α signalling pathway and NF-κB pathway, which were found to be closely associated with liver injury in MAFLD in our study (Figure 3) and other studies (Zhao et al., 2004; Schmöcker et al., 2007; Tapia et al., 2014). PUFAs, primarily n-3 PUFAs, can inhibit activation of the NF-κB pathway and suppress the inflammatory response (Calder, 2015). Moreover, additional supplementation with n-3 PUFAs can significantly reduce the C-reactive protein, IL-6, and TNF-a levels and further block the progression of western-diet induced MAFLD (Li et al., 2014; Lytle et al., 2017). Regarding lipid metabolism, dietary PUFAs were found to stimulate fatty acid oxidation by binding directly to peroxisome proliferator-activated receptor-α(Pawar et al., 2002). A randomized controlled trial also assessed the association of dietary PUFAs with lipid metabolism and peripheral insulin sensitivity in vivo (Hodson et al., 2017). After a 15-18 months intervention, individuals who achieved higher enrichment of erythrocyte docosahexaenoic acid, an PUFA, presented significantly improved hepatic insulin sensitivity, decreased fasting and postprandial plasma triglyceride concentrations, and reduced fasting hepatic de novo lipogenesis (Hodson et al., 2017).
Due to the therapeutic role of PUFAs in treating MAFLD, further studies on the core gene regulatory network underlying the PUFAs-MAFLD association would help reveal therapeutic targets for the disease. For the first time, we mined the hub genes of the PUFAs-related gene modules (Table 1) and validated the results in mouse experiments. We found that the expression levels of adamts1 in the brown inflammation-related module and tgfβ3 in the yellow lipid metabolism-related module were significantly decreased in mice with MAFLD. These results suggest that adamts1 and tgfβ3 may play important roles in the association between PUFAs and MAFLD.
ADAMTS1, the first identified Adamts family member, is characterized by its ability to cleave proteoglycans, aggrecan, collagen, elastin, and other extracellular matrix proteins via its metalloprotease-dependent catalytic (Rodriguez-Manzaneque et al., 2009; Esselens et al., 2010) and thrombospondin-dependent regions (Luque et al., 2003); thus, it plays roles in degrading extracellular matrix (ECM) components or inhibiting angiogenesis. A recent study reported that ablation of adamts1 in adipose tissue led to enlarged adipose tissue mass, reduced insulin sensitivity, and dysregulated lipid metabolism (Chen et al., 2016). In addition, lower mRNA levels of adamts1 in both subcutaneous and visceral white adipose tissues were associated with higher body mass index in humans (Chen et al., 2016). The reason underlying the above phenomenon might be due to adamts1 enzymatically impairing adipogenesis via ECM remodelling (Chen et al., 2016), which indicated the role of ADAMTS1 in maintaining lipids homeostasis. It was also reported that ADAMTS1 served as an extracellular protease that could be activated by prostaglandin F2α, a metabolite of the arachidonic acid (a type of PUFAs) metabolism (Keightley et al., 2010). In addition, ADAMTS1 also identified as an inflammatory associated protein is required for a balanced immune response (Kuno et al., 1997; Rodríguez-Baena et al., 2018). Thus, we speculate that ADAMTS1 decreasing in the liver might indicate abnormal immune response and reflect adipogenesis inhibition via ECM remodelling and PUFA regulation. However, further studies are warranted to prove the hypothesis.
TGFβ3, a member of the TGFβ superfamily, plays a multifunctional role in the immune response and ECM formation and regulates cell fate. The mRNA level of tgfβ3 was shown to be significantly lower in the HFD-fed mice in our study. Consistent with our study, TGFβ3 was reported to maintain lipid homeostasis. For example, overexpression of TGFβ3 was reported to reduce the histopathological damage observed in liver fibrosis (Zhang et al., 2010). Another study suggested that TGFβ3 participates in adipose tissue hypertrophy by regulating adipocyte precursor cell proliferation (Petrus et al., 2018). Moreover, mice with tgfβ3 haploinsufficiency that are fed an HFD show increased glucose intolerance and weight gain in vivo compared to their wild-type littermates. TGFβ3 was also proven to play an anti-inflammatory role and might participate in obesity and insulin resistance. TGFβ3 secretion was also found to be decreased in a genetic model of mice with CD4+ T cell-specific KLF10 knockout in vitro and in vivo, and these mice have a predisposition to obesity, insulin resistance, and fatty liver (Wara et al., 2020). It is quite interesting that TGFβ3 could mediate the secretion of prostaglandin E2, which is also an important metabolite in the arachidonic acid metabolic pathway, to affect macrophage polarization. Interestingly, TGFβ3 not only participates in the inflammatory response (Okamura et al., 2015) but also improves glucose tolerance and phenotypic changes in adipocyte morphology (Hall et al., 2013), suggesting that tgfβ3 may offer potential therapeutic benefits for MAFLD.
Our study has many strengths. First, we used bioinformatics tools, including DEGs analyses and WGCNA, combined with wet-lab experiments to propose and verify our hypothesis. Second, for the first time, we provide evidence that two hub genes, namely, adamts1 and tgfβ3, might play important roles in the association of PUFAs with MAFLD. Nevertheless, this study has limitations. First, although we mined hub genes that might play a role in the PUFAs-MAFLD association, the MAFLD mouse model we used in the validation experiment was established by feeding on HFD. A model of PUFAs diet supplementation could help us further validate our hypothesis. Second, the specific relationship between DEGs and the development of MAFLD should be clarified. Further insight into the molecular function of adamts1 and tgfβ3 in PUFA metabolism and MAFLD pathogenesis is needed.
Conclusion
In conclusion, PUFAs were associated with a brown inflammation-related gene module and a yellow lipid metabolism-related module of MAFLD. The expression levels of adamts1 in the brown module and tgβ3 in the yellow module were downregulated in MAFLD mice. These findings provide evidence regarding the role of adamts1 and tgfβ3 in the PUFAs-MAFLD association. Overall, the results of this comprehensive analysis have implications for personalized medicine and are of great clinical significance.
Data availability statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: http://www.ncbi.nlm.nih.gov/geo/; GSE89632.
Ethics statement
The animal study was reviewed and approved by the Animal Care and Ethics Committee at Peking Union Medical College Hospital.
Author contributions
CX and SC developed the concept of the research and wrote the original draft of the manuscript; CX, CY, and SC collected data and performed the analysis, JL and MY supervised and secured funds for the research. All authors read and approved the final manuscript.
Funding
This study was funded by the grants from National Key Research and Development Program of China (2020YFC2004505, 2018YFC2001105), National Natural Science Foundation of China (No. 82170855, No. 81900723), and Fundamental Research Funds for the Central Universities (3332021094).
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2022.951224/full#supplementary-material
SUPPLEMENTARY FIGURE S1 | Combination of DEGs between HC-SS and HC-NASH.
SUPPLEMENTARY FIGURE S2 | Sample clustering before network construction. Two outliers were shown by sample clustering. (B) Outliers were excluded in further analysis.
SUPPLEMENTARY FIGURE S3 | Module membership (MM) versus gene significance (GS) plot. Genes with higher module membership showed stronger gene significance for PUFAs.
References
Adams, L. A., Anstee, Q. M., Tilg, H., and Targher, G. (2017). Non-alcoholic fatty liver disease and its relationship with cardiovascular disease and other extrahepatic diseases. Gut 66 (6), 1138–1153. doi:10.1136/gutjnl-2017-313884
Anstee, Q. M., Targher, G., and Day, C. P. (2013). Progression of NAFLD to diabetes mellitus, cardiovascular disease or cirrhosis. Nat. Rev. Gastroenterol. Hepatol. 10 (6), 330–344. doi:10.1038/nrgastro.2013.41
Arendt, B. M., Comelli, E. M., Ma, D. W., Lou, W., Teterina, A., Kim, T., et al. (2015). Altered hepatic gene expression in nonalcoholic fatty liver disease is associated with lower hepatic n-3 and n-6 polyunsaturated fatty acids. Hepatology 61 (5), 1565–1578. doi:10.1002/hep.27695
Calder, P. C. (2015). Marine omega-3 fatty acids and inflammatory processes: Effects, mechanisms and clinical relevance. Biochim. Biophys. Acta 1851 (4), 469–484. doi:10.1016/j.bbalip.2014.08.010
Chan, T. W., Fu, T., Bahn, J. H., Jun, H. I., Lee, J. H., Quinones-Valdez, G., et al. (2020). RNA editing in cancer impacts mRNA abundance in immune response pathways. Genome Biol. 21 (1), 268. doi:10.1186/s13059-020-02171-4
Chawla, A., Nguyen, K. D., and Goh, Y. P. (2011). Macrophage-mediated inflammation in metabolic disease. Nat. Rev. Immunol. 11 (11), 738–749. doi:10.1038/nri3071
Chen, S. Z., Ning, L. F., Xu, X., Jiang, W. Y., Xing, C., Jia, W. P., et al. (2016). The miR-181d-regulated metalloproteinase Adamts1 enzymatically impairs adipogenesis via ECM remodeling. Cell Death Differ. 23 (11), 1778–1791. doi:10.1038/cdd.2016.66
Da Silva, H. E., Arendt, B. M., Noureldin, S. A., Therapondos, G., Guindi, M., and Allard, J. P. (2014). A cross-sectional study assessing dietary intake and physical activity in Canadian patients with nonalcoholic fatty liver disease vs healthy controls. J. Acad. Nutr. Diet. 114 (8), 1181–1194. doi:10.1016/j.jand.2014.01.009
Day, C. P., and James, O. F. (1998). Steatohepatitis: A tale of two "hits. Gastroenterology 114 (4), 842–845. doi:10.1016/s0016-5085(98)70599-2
De Taeye, B. M., Novitskaya, T., McGuinness, O. P., Gleaves, L., Medda, M., Covington, J. W., et al. (2007). Macrophage TNF-alpha contributes to insulin resistance and hepatic steatosis in diet-induced obesity. Am. J. Physiol. Endocrinol. Metab. 293 (3), E713–E725. doi:10.1152/ajpendo.00194.2007
Di Minno, M. N., Russolillo, A., Lupoli, R., Ambrosino, P., Di Minno, A., and Tarantino, G. (2012). Omega-3 fatty acids for the treatment of non-alcoholic fatty liver disease. World J. Gastroenterol. 18 (41), 5839–5847. doi:10.3748/wjg.v18.i41.5839
Dobson, S. M., García-Prat, L., Vanner, R. J., Wintersinger, J., Waanders, E., Gu, Z., et al. (2020). Relapse-fated latent diagnosis subclones in acute B lineage leukemia are drug tolerant and possess distinct metabolic programs. Cancer Discov. 10 (4), 568–587. doi:10.1158/2159-8290.Cd-19-1059
Dowman, J. K., Tomlinson, J. W., and Newsome, P. N. (2010). Pathogenesis of non-alcoholic fatty liver disease. Qjm 103 (2), 71–83. doi:10.1093/qjmed/hcp158
Eslam, M., Sanyal, A. J., George, J., and International Consensus, P. (2020). Mafld: A consensus-driven proposed nomenclature for metabolic associated fatty liver disease. Gastroenterology 158 (7), 1999–2014. doi:10.1053/j.gastro.2019.11.312
Esselens, C., Malapeira, J., Colome, N., Casal, C., Rodriguez-Manzaneque, J. C., Canals, F., et al. (2010). The cleavage of semaphorin 3C induced by ADAMTS1 promotes cell migration. J. Biol. Chem. 285 (4), 2463–2473. doi:10.1074/jbc.M109.055129
Hall, B. E., Wankhade, U. D., Konkel, J. E., Cherukuri, K., Nagineni, C. N., Flanders, K. C., et al. (2013). Transforming growth factor-β3 (TGF-β3) knock-in ameliorates inflammation due to TGF-β1 deficiency while promoting glucose tolerance. J. Biol. Chem. 288 (44), 32074–32092. doi:10.1074/jbc.M113.480764
Harley, I. T., Stankiewicz, T. E., Giles, D. A., Softic, S., Flick, L. M., Cappelletti, M., et al. (2014). IL-17 signaling accelerates the progression of nonalcoholic fatty liver disease in mice. Hepatology 59 (5), 1830–1839. doi:10.1002/hep.26746
Hernández-Gea, V., Campreciós, G., Betancourt, F., Pérez-Campuzano, V., Seijo, S., Díaz, A., et al. (2021). Co-expression gene network analysis reveals novel regulatory pathways involved in porto-sinusoidal vascular disease. J. Hepatol. 75 (4), 924–934. doi:10.1016/j.jhep.2021.05.014
Hodson, L., Bhatia, L., Scorletti, E., Smith, D. E., Jackson, N. C., Shojaee-Moradie, F., et al. (2017). Docosahexaenoic acid enrichment in NAFLD is associated with improvements in hepatic metabolism and hepatic insulin sensitivity: A pilot study. Eur. J. Clin. Nutr. 71 (8), 1251–1979. doi:10.1038/ejcn.2017.145
Huang, W., Metlakunta, A., Dedousis, N., Zhang, P., Sipula, I., Dube, J. J., et al. (2010). Depletion of liver Kupffer cells prevents the development of diet-induced hepatic steatosis and insulin resistance. Diabetes 59 (2), 347–357. doi:10.2337/db09-0016
Imamura, F., Micha, R., Wu, J. H., de Oliveira Otto, M. C., Otite, F. O., Abioye, A. I., et al. (2016). Effects of saturated fat, polyunsaturated fat, monounsaturated fat, and carbohydrate on glucose-insulin homeostasis: A systematic review and meta-analysis of randomised controlled feeding trials. PLoS Med. 13 (7), e1002087. doi:10.1371/journal.pmed.1002087
Keightley, M. C., Sales, K. J., and Jabbour, H. N. (2010). PGF2α-F-prostanoid receptor signalling via ADAMTS1 modulates epithelial cell invasion and endothelial cell function in endometrial cancer. BMC Cancer 10, 488. doi:10.1186/1471-2407-10-488
Koyama, Y., and Brenner, D. A. (2017). Liver inflammation and fibrosis. J. Clin. Invest. 127 (1), 55–64. doi:10.1172/jci88881
Kuno, K., Kanada, N., Nakashima, E., Fujiki, F., Ichimura, F., and Matsushima, K. (1997). Molecular cloning of a gene encoding a new type of metalloproteinase-disintegrin family protein with thrombospondin motifs as an inflammation associated gene. J. Biol. Chem. 272 (1), 556–562. doi:10.1074/jbc.272.1.556
Langfelder, P., and Horvath, S. (2008). Wgcna: an R package for weighted correlation network analysis. BMC Bioinforma. 9, 559. doi:10.1186/1471-2105-9-559
Lanthier, N., Molendi-Coste, O., Cani, P. D., van Rooijen, N., Horsmans, Y., and Leclercq, I. A. (2011). Kupffer cell depletion prevents but has no therapeutic effect on metabolic and inflammatory changes induced by a high-fat diet. Faseb J. 25 (12), 4301–4311. doi:10.1096/fj.11-189472
Lanthier, N., Molendi-Coste, O., Horsmans, Y., van Rooijen, N., Cani, P. D., and Leclercq, I. A. (2010). Kupffer cell activation is a causal factor for hepatic insulin resistance. Am. J. Physiol. Gastrointest. Liver Physiol. 298 (1), G107–G116. doi:10.1152/ajpgi.00391.2009
Li, K., Huang, T., Zheng, J., Wu, K., and Li, D. (2014). Effect of marine-derived n-3 polyunsaturated fatty acids on C-reactive protein, interleukin 6 and tumor necrosis factor α: A meta-analysis. PLoS One 9 (2), e88103. doi:10.1371/journal.pone.0088103
Li, Y. H., Yang, L. H., Sha, K. H., Liu, T. G., Zhang, L. G., and Liu, X. X. (2015). Efficacy of poly-unsaturated fatty acid therapy on patients with nonalcoholic steatohepatitis. World J. Gastroenterol. 21 (22), 7008–7013. doi:10.3748/wjg.v21.i22.7008
Liang, W., Luo, Q., Zhang, Z., Yang, K., Yang, A., Chi, Q., et al. (2022). An integrated bioinformatics analysis and experimental study identified key biomarkers CD300A or CXCL1, pathways and immune infiltration in diabetic nephropathy mice. BIOCELL 46 (8), 1989–2002. doi:10.32604/biocell.2022.019300
Lim, S., Kim, J. W., and Targher, G. (2021). Links between metabolic syndrome and metabolic dysfunction-associated fatty liver disease. Trends Endocrinol. Metab. 32 (7), 500–514. doi:10.1016/j.tem.2021.04.008
Lin, C. Y., Lee, C. H., Chuang, Y. H., Lee, J. Y., Chiu, Y. Y., Wu Lee, Y. H., et al. (2019). Membrane protein-regulated networks across human cancers. Nat. Commun. 10 (1), 3131. doi:10.1038/s41467-019-10920-8
Long, J., Huang, S., Bai, Y., Mao, J., Wang, A., Lin, Y., et al. (2021). Transcriptional landscape of cholangiocarcinoma revealed by weighted gene coexpression network analysis. Brief. Bioinform. 22 (4), bbaa224. doi:10.1093/bib/bbaa224
Luque, A., Carpizo, D. R., and Iruela-Arispe, M. L. (2003). ADAMTS1/METH1 inhibits endothelial cell proliferation by direct binding and sequestration of VEGF165. J. Biol. Chem. 278 (26), 23656–23665. doi:10.1074/jbc.M212964200
Lytle, K. A., Wong, C. P., and Jump, D. B. (2017). Docosahexaenoic acid blocks progression of Western diet-induced nonalcoholic steatohepatitis in obese Ldlr-/- mice. PLoS One 12 (4), e0173376. doi:10.1371/journal.pone.0173376
Machado, M. V., and Diehl, A. M. (2016). Pathogenesis of nonalcoholic steatohepatitis. Gastroenterology 150 (8), 1769–1777. doi:10.1053/j.gastro.2016.02.066
Meng, F., Wang, K., Aoyama, T., Grivennikov, S. I., Paik, Y., Scholten, D., et al. (2012). Interleukin-17 signaling in inflammatory, Kupffer cells, and hepatic stellate cells exacerbates liver fibrosis in mice. Gastroenterology 143 (3), 765765–765776. doi:10.1053/j.gastro.2012.05.049
Okamura, T., Morita, K., Iwasaki, Y., Inoue, M., Komai, T., Fujio, K., et al. (2015). Role of TGF-β3 in the regulation of immune responses. Clin. Exp. Rheumatol. 33 (4), S63–S69.
Park, S. K., Seo, M. H., Shin, H. C., and Ryoo, J. H. (2013). Clinical availability of nonalcoholic fatty liver disease as an early predictor of type 2 diabetes mellitus in Korean men: 5-year prospective cohort study. Hepatology 57 (4), 1378–1383. doi:10.1002/hep.26183
Pawar, A., Xu, J., Jerks, E., Mangelsdorf, D. J., and Jump, D. B. (2002). Fatty acid regulation of liver X receptors (LXR) and peroxisome proliferator-activated receptor alpha (PPARalpha ) in HEK293 cells. J. Biol. Chem. 277 (42), 39243–39250. doi:10.1074/jbc.M206170200
Petrus, P., Mejhert, N., Corrales, P., Lecoutre, S., Li, Q., Maldonado, E., et al. (2018). Transforming growth factor-β3 regulates adipocyte number in subcutaneous white adipose tissue. Cell Rep. 25 (3), 551–560. e555. doi:10.1016/j.celrep.2018.09.069
Peverill, W., Powell, L. W., and Skoien, R. (2014). Evolving concepts in the pathogenesis of NASH: Beyond steatosis and inflammation. Int. J. Mol. Sci. 15 (5), 8591–8638. doi:10.3390/ijms15058591
Puri, P., Baillie, R. A., Wiest, M. M., Mirshahi, F., Choudhury, J., Cheung, O., et al. (2007). A lipidomic analysis of nonalcoholic fatty liver disease. Hepatology 46 (4), 1081–1090. doi:10.1002/hep.21763
Rodríguez-Baena, F. J., Redondo-García, S., Peris-Torres, C., Martino-Echarri, E., Fernández-Rodríguez, R., Plaza-Calonge, M. D. C., et al. (2018). ADAMTS1 protease is required for a balanced immune cell repertoire and tumour inflammatory response. Sci. Rep. 8 (1), 13103. doi:10.1038/s41598-018-31288-7
Rodriguez-Manzaneque, J. C., Carpizo, D., Plaza-Calonge Mdel, C., Torres-Collado, A. X., Thai, S. N., Simons, M., et al. (2009). Cleavage of syndecan-4 by ADAMTS1 provokes defects in adhesion. Int. J. Biochem. Cell Biol. 41 (4), 800–810. doi:10.1016/j.biocel.2008.08.014
Rolo, A. P., Teodoro, J. S., and Palmeira, C. M. (2012). Role of oxidative stress in the pathogenesis of nonalcoholic steatohepatitis. Free Radic. Biol. Med. 52 (1), 59–69. doi:10.1016/j.freeradbiomed.2011.10.003
Schmöcker, C., Weylandt, K. H., Kahlke, L., Wang, J., Lobeck, H., Tiegs, G., et al. (2007). Omega-3 fatty acids alleviate chemically induced acute hepatitis by suppression of cytokines. Hepatology 45 (4), 864–869. doi:10.1002/hep.21626
Schulze, M. B., Minihane, A. M., Saleh, R. N. M., and Risérus, U. (2020). Intake and metabolism of omega-3 and omega-6 polyunsaturated fatty acids: Nutritional implications for cardiometabolic diseases. Lancet. Diabetes Endocrinol. 8 (11), 915–930. doi:10.1016/s2213-8587(20)30148-0
Schuster, S., Cabrera, D., Arrese, M., and Feldstein, A. E. (2018). Triggering and resolution of inflammation in NASH. Nat. Rev. Gastroenterol. Hepatol. 15 (6), 349–364. doi:10.1038/s41575-018-0009-6
Stienstra, R., Saudale, F., Duval, C., Keshtkar, S., Groener, J. E., van Rooijen, N., et al. (2010). Kupffer cells promote hepatic steatosis via interleukin-1beta-dependent suppression of peroxisome proliferator-activated receptor alpha activity. Hepatology 51 (2), 511–522. doi:10.1002/hep.23337
Tapia, G., Valenzuela, R., Espinosa, A., Romanque, P., Dossi, C., Gonzalez-Mañán, D., et al. (2014). N-3 long-chain PUFA supplementation prevents high fat diet induced mouse liver steatosis and inflammation in relation to PPAR-α upregulation and NF-κB DNA binding abrogation. Mol. Nutr. Food Res. 58 (6), 1333–1341. doi:10.1002/mnfr.201300458
Wang, C., Han, J., Liu, M., Huang, Y., Zhou, T., Jiang, N., et al. (2022). RNA-sequencing of human aortic valves identifies that miR-629-3p and TAGLN miRNA-mRNA pair involving in calcified aortic valve disease. J. Physiol. Biochem. 1, 1. doi:10.1007/s13105-022-00905-5
Wara, A. K., Wang, S., Wu, C., Fang, F., Haemmig, S., Weber, B. N., et al. (2020). KLF10 deficiency in CD4(+) T cells triggers obesity, insulin resistance, and fatty liver. Cell Rep. 33 (13), 108550. doi:10.1016/j.celrep.2020.108550
Xiu, Y., Dong, Q., Fu, L., Bossler, A., Tang, X., Boyce, B., et al. (2020). Coactivation of NF-κB and Notch signaling is sufficient to induce B-cell transformation and enables B-myeloid conversion. Blood 135 (2), 108–120. doi:10.1182/blood.2019001438
Xu, R., Wu, Q., Gong, Y., Wu, Y., Chi, Q., and Sun, D. (2022). A novel prognostic target-gene signature and nomogram based on an integrated bioinformatics analysis in hepatocellular carcinoma. BIOCELL 46 (5), 1261–1288. doi:10.32604/biocell.2022.018427
Younossi, Z., Tacke, F., Arrese, M., Chander Sharma, B., Mostafa, I., Bugianesi, E., et al. (2019). Global perspectives on nonalcoholic fatty liver disease and nonalcoholic steatohepatitis. Hepatology 69 (6), 2672–2682. doi:10.1002/hep.30251
Zeng, F., Shi, M., Xiao, H., and Chi, X. (2021). WGCNA-based identification of hub genes and key pathways involved in nonalcoholic fatty liver disease. Biomed. Res. Int. 2021, 5633211. doi:10.1155/2021/5633211
Zhang, Y., Liu, P., Gao, X., Qian, W., and Xu, K. (2010). rAAV2-TGF-β(3) decreases collagen synthesis and deposition in the liver of experimental hepatic fibrosis rat. Dig. Dis. Sci. 55 (10), 2821–2830. doi:10.1007/s10620-009-1119-3
Keywords: WGCNA, metabolic dysfunction-associated fatty liver disease, polyunsaturated fatty acids, differentially expressed genes, ADAMTS1 gene, TGFβ3 gene
Citation: Xiao C, Chen S, Yang C, Liu J and Yu M (2022) Identification of polyunsaturated fatty acids related key modules and genes in metabolic dysfunction-associated fatty liver disease using WGCNA analysis. Front. Genet. 13:951224. doi: 10.3389/fgene.2022.951224
Received: 23 May 2022; Accepted: 26 October 2022;
Published: 08 November 2022.
Edited by:
Babak Behnam, National Sanitation Foundation International, United StatesReviewed by:
Xuelian Xiong, Fudan University, ChinaQingjia Chi, Wuhan University of Technology, China
Copyright © 2022 Xiao, Chen, Yang, Liu and Yu. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Jieying Liu, liujieying@vip.163.com; Miao Yu, yumiaoxh@163.com
†These authors have contributed equally to this work and share first authorship