WGCNA-Based Identification of Hub Genes and Key Pathways Involved in Nonalcoholic Fatty Liver Disease

Abstract

Background. The morbidity of nonalcoholic fatty liver disease (NAFLD) has been rising, but the pathogenesis of NAFLD is still elusive. This study is aimed at determining NAFLD-related hub genes based on weighted gene coexpression network analysis (WGCNA). Methods. GSE126848 dataset based construction of coexpression networks was performed based on WGCNA. Database for Annotation, Visualization, and Integrated Discovery (DAVID) was utilized for Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis. Hub genes were identified and validated in independent datasets and mouse model. Results. We found that the steelblue module was most significantly correlated with NAFLD. Total 15 hub genes (NDUFA9, UQCRQ, NDUFB8, COPS5, RPS17, UBL5, PSMA3, PSMA1, SF3B5, MRPL27, RPL26, PDCD5, PFDN6, SNRPD2, PSMB3) were derived from both the coexpression and PPI networks and considered “true” hub genes. Functional enrichment analysis showed that the hub genes were related to NAFLD pathway and oxidative phosphorylation. Independent dataset-based analysis and the establishment of NAFLD mouse model confirmed the involvement of two hub genes NDUFA9 and UQCRQ in the pathogenesis of NAFLD. Conclusions. Oxidative phosphorylation and NAFLD pathway may be crucially involved in the pathogenesis of NAFLD, and two hub genes NDUFA9 and UQCRQ might be diagnostic biomarkers and therapeutic targets for NAFLD.

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

NAFLD is one of the most prevalent chronic liver diseases in the world, affecting approximately 30% of the adult population globally [1]. Based on disease severity, NAFLD is classified as nonalcoholic fatty liver (NAFL) or nonalcoholic steatohepatitis (NASH). Patients with NASH exhibit histological lobular inflammation and hepatocyte ballooning, while NASH may progress to fibrosis faster than NAFL [2]. Patients with NAFLD often develop cirrhosis and numerous liver-related complications. In the USA, NASH is the third leading cause of end-stage liver disease and hepatocellular malignancy [3].

NAFLD has drawn remarkable attention due to its widespread prevalence and socioeconomic burden. Currently, no drug or surgery has been approved for the treatment of NAFLD, and weight loss is the only proven option for the management of NAFLD [4]. Recently, bioinformatics tools based on gene expression profiling via high-throughput microarray technology have been applied to elucidate the pathogenesis of a variety of diseases, including NAFLD [5–8]. However, most of the studies only focused on the screening of differentially expressed genes (DEGs) rather than the functional connections among these genes by gene expression pattern analysis [6].

NAFLD-related DEGs can be identified using traditional bioinformatics methods, including the Limma software package. However, these procedures could lead to neglect of the genes that have little difference in fold change but contribute to the pathogenesis of NAFLD. Given that differential gene expression levels may not reflect the complexity of NAFLD, we employed a new analytical method WGCNA (Weighted Gene Co-expression Network Analysis) to determine essential genes and signaling pathways involved in the pathogenesis of NAFLD.

WGCNA uses data from thousands of the most varied genes or all the genes to identify gene sets of interest compared to genes that are only concerned with differential expression while analyzing significant associations with phenotypes. WGCNA has the following two advantages: one is to lose fewer genes; the other is to collect a large number of genes into several gene sets and associate them with phenotypes without multiple hypothesis tests [9, 10]. Therefore, we used WGCNA to analyze the gene expression synthesis database of gene expression data set GSE126848 to identify hub genes and molecular pathways involved in NAFLD.

2. Materials and Methods

2.1. Data Collection

Data of GSE126848 based on the GPL18573 platform (Illumina NextSeq 500 (Homo sapiens)) were obtained from Gene Expression Omnibus (GEO) database of NCBI (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE126848). GSE126848 contained 57 liver biopsy samples from 14 normal healthy and 12 obese individuals as well as 15 NAFL and 16 NASH subjects. The dataset used a second-generation sequencing method to provide count data after the map. The workflow of this study involving coexpression network construction, hub genes identification, functional analysis, and validation is shown in Figure 1.

[figure omitted; refer to PDF]

Next, we detected the expression of three candidate hub genes in the livers of NAFLD and control groups. We found significantly decreased expression of NDUFA9 and UQCRQ in the NAFLD model group compared to the control group, but no significant difference in the expression of NDUFB8 between NAFLD and control groups (Figure 6).

[figure omitted; refer to PDF]

4. Discussion

In recent years, a number of candidate causal genes for NAFLD, including PNPLA3, PPP1R3B, SAMM50, and TRIB1, have been identified based on the genome-wide association studies [24–26]. However, the exact mechanisms underlying the development of NAFLD remain to be determined. Most studies focused on the comparison of the gene expression between NAFLD patients and normal individuals. For example, Hoang et al. [27] and Frades et al. [7] employed DEG-seq to analyze the differentially expressed genes but did not perform dataset or laboratory verification. In contrast, in this study, we verified the selected important hub genes not only in an independent dataset but also in a high-fat feed-induced NAFLD mouse model.

To determine key genes and pathways in the pathogenesis of NAFLD, we performed WGCNA and identified fifteen “true” hub genes that are shared in both the coexpression and PPI networks. Furthermore, KEGG pathway analysis of these “true” hub genes revealed that three genes NDUFB8, NDUFA9, and UQCRQ were enriched in NAFLD and oxidative phosphorylation highly related to NAFLD pathogenesis. Therefore, we speculated that these three genes may be the key genes involved in the pathogenesis of NAFLD. In another independent dataset GSE89632, although the difference in logFC values between NDUFA9 and UQCRQ was weak, the downward trend was obvious. qRT-PCR assay verified that NDUFA9 and UQCRQ were the hub genes. Furthermore, we demonstrated significantly decreased expression of NDUFA9 and UQCRQ in NAFLD mouse model compared to the control.

NDUFA9 (NADH: ubiquinone oxidoreductase subunit A9) gene encodes a subunit of the mitochondrial membrane respiratory chain NADH dehydrogenase (complex I). Instead of being implicated in the catalysis, NDUFA9 is required for proper assembly of the complex I [28]. Wang et al. [29] reported that metformin protected the mitochondrial structure of mouse colorectal epithelial cells by increasing NDUFA9 expression to inhibit colitis and colitis-related cancers. Another study [ 30] showed that pioglitazone could bind to the complex I subunits NDUFA9, NDUFB6, and NDUFA6, inducing the complex disassembly, decreasing its activity, and promoting the expression of nuclear DNA-encoded subunits of the complex I in mice and HepG2 cells. Houtkooper et al. [31] reported that NR supplementation in mammalian cells and mouse tissues led to an increase in NAD(+) levels as well as the activation of SIRT1 and SIRT3, facilitating oxidative metabolism and protecting against metabolic abnormalities induced by a high-fat diet. These data suggest that NDUFA9 may be involved in NAFLD by regulating the activity of complex I subunits.

KEGG analysis revealed that the three hub genes identified in this study are predominantly enriched in the oxidative phosphorylation pathway. Oxidative phosphorylation represents a process by which the energy of adenosine triphosphate (ATP) can be efficiently produced. Oxidative phosphorylation was decreased in individuals with NAFLD [32]. Moreover, HFD-fed mice showed reduced activities of oxidative phosphorylation enzymes, which can be attributed to decreased amount of fully assembled complexes [33]. Oxidative phosphorylation dysfunction has been shown to cause oxidative stress involving cytochrome P4502E1, xanthine oxidase, NADPH oxidase, and liver mitochondria [34–36].

UQCRQ (ubiquinol-cytochrome c reductase, complex III subunit VII) is a gene encoding a low-molecular mass ubiquinone-binding protein. UQCRQ is identified as a small core-associated protein and a subunit of ubiquinol-cytochrome c reductase complex III [37]. Although no report showed that UQCRQ is directly related to NAFLD, UQCRQ is associated with oxidative phosphorylation pathway [38]. Based on the relationship between oxidative phosphorylation and NAFLD, we speculate that UQCRQ is linked to NAFLD and may be a key gene in the pathogenesis of NAFLD.

In this study, we identified the key regulatory genes and pathways involved in NAFLD using the integrated method of bioinformatics, including WGCNA, functional genomics, and gene regulatory network. We found that key genes and pathways such as NDUFA9, UQCRQ, and ribosome, as well as proteasome and oxidative phosphorylation are essential to the pathogenesis of NAFLD. The present study has the following limitations: First, our data obtained from the WGCNA remain to be validated in an independent cohort due to the single platform of dataset used in the study. Second, the absence of clinical traits such as serum biomarker profiling, liver biopsy histology for scoring hepatic steatosis, and fibrosis in the raw data might affect the assessment of NAFLD phenotypes. Third, the use of mouse rather than human samples may restrict the capacity to detect the differences in gene expression in NAFLD patients.

In conclusion, based on WGCNA analyses, we identified 15 NAFLD-related candidate hub genes in the steelblue module. Moreover, we found that the proteasome, oxidative phosphorylation, NAFLD, and ribosome may be involved in the pathogenesis of NAFLD. Importantly, we established mouse model of NAFLD and verified two hub genes NDUFA9 and UQCRQ, which may act as biomarkers and therapeutic targets for NAFLD.

Authors’ Contributions

FZ, MS, and HX performed the experiments. XC designed the study. All authors read and approved the manuscript.

References

[1] K. Garber, "The new liver epidemic," Nature biotechnology, vol. 37 no. 3, pp. 209-214, 2019.

[2] A. M. Diehl, C. Day, "Cause, pathogenesis, and treatment of nonalcoholic steatohepatitis," The New England journal of medicine, vol. 377 no. 21, pp. 2063-2072, DOI: 10.1056/NEJMra1503519, 2017.

[3] Z. M. Younossi, M. Otgonsuren, L. Henry, C. Venkatesan, A. Mishra, M. Erario, S. Hunt, "Association of nonalcoholic fatty liver disease (NAFLD) with hepatocellular carcinoma (HCC) in the United States from 2004 to 2009," Hepatology, vol. 62 no. 6, pp. 1723-1730, DOI: 10.1002/hep.28123, 2015.

[4] A. Mazzotti, M. T. Caletti, L. Brodosi, S. Di Domizio, M. L. Forchielli, S. Petta, E. Bugianesi, G. Bianchi, G. Marchesini, "An internet-based approach for lifestyle changes in patients with NAFLD: Two- year effects on weight loss and surrogate markers," Journal of hepatology, vol. 69 no. 5, pp. 1155-1163, DOI: 10.1016/j.jhep.2018.07.013, 2018.

[5] X. Wang, T. Yamamoto, M. Kadowaki, Y. Yang, "Identification of key pathways and gene expression in the activation of mast cells via calcium flux using bioinformatics analysis," Biocell, vol. 45 no. 2, pp. 395-415, 2021.

[6] B. M. Arendt, E. M. Comelli, D. W. Ma, W. Lou, A. Teterina, T. Kim, S. K. Fung, D. K. Wong, I. McGilvray, S. E. Fischer, J. P. Allard, "Altered hepatic gene expression in nonalcoholic fatty liver disease is associated with lower hepatic n-3 and n-6 polyunsaturated fatty acids," Hepatology, vol. 61 no. 5, pp. 1565-1578, DOI: 10.1002/hep.27695, 2015.

[7] H. Shen, S. Ren, W. Wang, C. Zhang, H. Hao, Q. Shen, Y. Duan, Z. Wang, W. Ge, "Profiles of immune status and related pathways in sepsis: evidence based on GEO and bioinformatics," Biocell, vol. 44 no. 4, pp. 583-589, 2020.

[8] P. Langfelder, S. Horvath, "WGCNA: an R package for weighted correlation network analysis," BMC Bioinformatics, vol. 9, 2008.

[9] X. Xiao, A. Moreno-Moral, M. Rotival, L. Bottolo, E. Petretto, "Multi-tissue analysis of co-expression networks by higher-order generalized singular value decomposition identifies functionally coherent transcriptional modules," PLoS genetics, vol. 10 no. 1,DOI: 10.1371/journal.pgen.1004006, 2014.

[10] W. Gao, Q. Ma, C. Tang, Y. Zhan, Y. Duan, H. Ni, Y. Xu, "Microenvironment and related genes predict outcomes of patients with cervical cancer: evidence from TCGA and bioinformatic analysis," Biocell, vol. 44 no. 4, pp. 597-605, 2020.

[11] C. Chen, L. Cheng, K. Grennan, F. Pibiri, C. Zhang, J. A. Badner, E. S. Gershon, C. Liu, "Two gene co-expression modules differentiate psychotics and controls," Molecular Psychiatry, vol. 18 no. 12, pp. 1308-1314, DOI: 10.1038/mp.2012.146, 2013.

[12] M. K. Terkelsen, S. M. Bendixen, D. Hansen, E. A. H. Scott, A. F. Moeller, R. Nielsen, S. Mandrup, A. Schlosser, T. L. Andersen, G. L. Sorensen, A. Krag, "Transcriptional dynamics of hepatic sinusoid-associated cells after liver injury," Hepatology, vol. 72 no. 6, pp. 2119-2133, DOI: 10.1002/hep.31215, 2020.

[13] E. Radulescu, A. E. Jaffe, R. E. Straub, Q. Chen, J. H. Shin, T. M. Hyde, J. E. Kleinman, D. R. Weinberger, "Identification and prioritization of gene sets associated with schizophrenia risk by co-expression network analysis in human brain," Molecular Psychiatry, vol. 25 no. 4, pp. 791-804, 2020.

[14] Y. Hu, L. Cheng, W. Zhong, M. Chen, Q. Zhang, "Bioinformatics analysis of gene expression profiles for risk prediction in patients with septic shock," Medical science monitor: international medical journal of experimental and clinical research, vol. 25, pp. 9563-9571, DOI: 10.12659/MSM.918491, 2019.

[15] L. Zhang, Y. Gong, S. Wang, F. Gao, "Anti-colorectal cancer mechanisms of formononetin identified by network pharmacological approach," Medical science monitor: international medical journal of experimental and clinical research, vol. 25, pp. 7709-7714, 2019.

[16] G. Yang, Y. Zhang, J. Yang, "A five-microRNA signature as prognostic biomarker in colorectal cancer by bioinformatics analysis," Frontiers in oncology, vol. 9,DOI: 10.3389/fonc.2019.01207, 2019.

[17] S. K. Saha, T. I. Jeon, S. B. Jang, S. J. Kim, K. M. Lim, Y. J. Choi, H. G. Kim, A. Kim, S. G. Cho, "Bioinformatics approach for identifying novel biomarkers and their signaling pathways involved in interstitial cystitis/bladder pain syndrome with Hunner lesion," Journal of Clinical Medicine, vol. 9 no. 6,DOI: 10.3390/jcm9061935, 2020.

[18] X. Zhang, H. Feng, Z. Li, D. Li, S. Liu, H. Huang, M. Li, "Application of weighted gene co-expression network analysis to identify key modules and hub genes in oral squamous cell carcinoma tumorigenesis," OncoTargets and therapy, vol. Volume 11, pp. 6001-6021, DOI: 10.2147/OTT.S171791, 2018.

[19] W. Huang da, B. T. Sherman, R. A. Lempicki, "Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources," Nature protocols, vol. 4 no. 1, pp. 44-57, 2009.

[20] D. W. Huang, B. T. Sherman, Q. Tan, J. R. Collins, W. G. Alvord, J. Roayaei, R. Stephens, M. W. Baseler, H. C. Lane, R. A. Lempicki, "The DAVID gene functional classification tool: a novel biological module-centric algorithm to functionally analyze large gene lists," Genome biology, vol. 8 no. 9,DOI: 10.1186/gb-2007-8-9-r183, 2007.

[21] J. Qin, T. Yang, N. Zeng, C. Wan, L. Gao, X. Li, L. Chen, Y. Shen, F. Wen, "Differential coexpression networks in bronchiolitis and emphysema phenotypes reveal heterogeneous mechanisms of chronic obstructive pulmonary disease," Journal of cellular and molecular medicine, vol. 23 no. 10, pp. 6989-6999, 2019.

[22] D. A. Giles, M. E. Moreno-Fernandez, T. E. Stankiewicz, S. Graspeuntner, M. Cappelletti, D. Wu, R. Mukherjee, C. C. Chan, M. J. Lawson, J. Klarquist, A. Sünderhauf, "Thermoneutral housing exacerbates nonalcoholic fatty liver disease in mice and allows for sex-independent disease modeling," Nature medicine, vol. 23 no. 7, pp. 829-838, 2017.

[23] M. I. Hernandez-Alvarez, D. Sebastian, S. Vives, S. Ivanova, P. Bartoccioni, P. Kakimoto, N. Plana, S. R. Veiga, V. Hernández, N. Vasconcelos, G. Peddinti, "Deficient endoplasmic reticulum-mitochondrial phosphatidylserine transfer causes liver disease," Cell, vol. 177 no. 4, pp. 881-95 e17, 2019.

[24] J. Kozlitina, E. Smagris, S. Stender, B. G. Nordestgaard, H. H. Zhou, A. Tybjaerg-Hansen, T. F. Vogt, H. H. Hobbs, J. C. Cohen, "Exome-wide association study identifies a TM6SF2 variant that confers susceptibility to nonalcoholic fatty liver disease," Nature genetics, vol. 46 no. 4, pp. 352-356, 2014.

[25] T. Kitamoto, A. Kitamoto, M. Yoneda, H. Hyogo, H. Ochi, T. Nakamura, H. Teranishi, S. Mizusawa, T. Ueno, K. Chayama, A. Nakajima, "Genome-wide scan revealed that polymorphisms in the PNPLA3, SAMM50, and PARVB genes are associated with development and progression of nonalcoholic fatty liver disease in Japan," Human Genetics, vol. 132 no. 7, pp. 783-792, 2013.

[26] J. C. Chambers, W. Zhang, J. Sehmi, X. Li, M. N. Wass, P. Van der Harst, H. Holm, S. Sanna, M. Kavousi, S. E. Baumeister, L. J. Coin, "Genome-wide association study identifies loci influencing concentrations of liver enzymes in plasma," Nature genetics, vol. 43 no. 11, pp. 1131-1138, 2011.

[27] S. A. Hoang, A. Oseini, R. E. Feaver, B. K. Cole, A. Asgharpour, R. Vincent, M. Siddiqui, M. J. Lawson, N. C. Day, J. M. Taylor, B. R. Wamhoff, "Gene expression predicts histological severity and reveals distinct molecular profiles of nonalcoholic fatty liver disease," Scientific Reports, vol. 9 no. 1, 2019.

[28] F. Baertling, L. Sanchez-Caballero, M. A. M. van den Brand, C. W. Fung, S. H. Chan, V. C. Wong, D. M. E. Hellebrekers, I. F. M. de Coo, J. A. M. Smeitink, R. J. T. Rodenburg, L. G. J. Nijtmans, "NDUFA9 point mutations cause a variable mitochondrial complex I assembly defect," Clinical genetics, vol. 93 no. 1, pp. 111-118, DOI: 10.1111/cge.13089, 2018.

[29] S. Q. Wang, S. X. Cui, X. J. Qu, "Metformin inhibited colitis and colitis-associated cancer (CAC) through protecting mitochondrial structures of colorectal epithelial cells in mice," Cancer biology & therapy, vol. 20 no. 3, pp. 338-348, 2019.

[30] I. Garcia-Ruiz, P. Solis-Munoz, D. Fernandez-Moreira, T. Munoz-Yague, J. A. Solis-Herruzo, "Pioglitazone leads to an inactivation and disassembly of complex I of the mitochondrial respiratory chain," BMC Biology, vol. 11, 2013.

[31] C. Canto, R. H. Houtkooper, E. Pirinen, D. Y. Youn, M. H. Oosterveer, Y. Cen, P. J. Fernandez-Marcos, H. Yamamoto, P. A. Andreux, P. Cettour-Rose, K. Gademann, "The NAD + Precursor Nicotinamide Riboside Enhances Oxidative Metabolism and Protects against High-Fat Diet-Induced Obesity," Cell metabolism, vol. 15 no. 6, pp. 838-847, DOI: 10.1016/j.cmet.2012.04.022, 2012.

[32] M. Perez-Carreras, P. Del Hoyo, M. A. Martin, J. C. Rubio, A. Martin, G. Castellano, F. Colina, J. Arenas, J. A. Solis-Herruzo, "Defective hepatic mitochondrial respiratory chain in patients with nonalcoholic steatohepatitis," Hepatology, vol. 38 no. 4, pp. 999-1007, 2003.

[33] I. Garcia-Ruiz, P. Solis-Munoz, D. Fernandez-Moreira, M. Grau, F. Colina, T. Munoz-Yague, J. A. Solís-Herruzo, "High-fat diet decreases activity of the oxidative phosphorylation complexes and causes nonalcoholic steatohepatitis in mice," Disease models & mechanisms, vol. 7 no. 11, pp. 1287-1296, 2014.

[34] H. K. Seitz, "The role of cytochrome P4502E1 in the pathogenesis of alcoholic liver disease and carcinogenesis," Chemico-biological interactions, vol. 316,DOI: 10.1016/j.cbi.2019.108918, 2020.

[35] J. Nanduri, D. R. Vaddi, S. A. Khan, N. Wang, V. Makerenko, N. R. Prabhakar, "Xanthine oxidase mediates hypoxia-inducible factor-2alpha degradation by intermittent hypoxia," PLoS One, vol. 8 no. 10, 2013.

[36] A. J. Sanyal, "Insulin resistance and nonalcoholic steatohepatitis: fat or fiction?," The American journal of gastroenterology., vol. 96 no. 2, pp. 274-276, 2001.

[37] O. Barel, Z. Shorer, H. Flusser, R. Ofir, G. Narkis, G. Finer, H. Shalev, A. Nasasra, A. Saada, O. S. Birk, "Mitochondrial complex III deficiency associated with a homozygous mutation in UQCRQ," American journal of human genetics, vol. 82 no. 5, pp. 1211-1216, 2008.

[38] R. Tian, S. Xu, S. Chai, D. Yin, H. Zakon, G. Yang, "Stronger selective constraint on downstream genes in the oxidative phosphorylation pathway of cetaceans," Journal of evolutionary biology, vol. 31 no. 2, pp. 217-228, 2018.