Candidate genes investigation for severe nonalcoholic fatty liver disease based on bioinformatics analysis.

BACKGROUND: Nonalcoholic fatty liver disease (NAFLD) is the most common chronic liver condition worldwide. However, its etiology and fundamental pathophysiology for the disease process are poorly understood. In this study, we thus used bioinformatics to identify candidate genes potentially causative of severe NAFLD.
METHODS: Gene expression profile data GSE49541 were downloaded from the Gene Expression Omnibus database. Tissues samples from 32 severe and 40 mild NAFLD patients were evaluated to identify differentially expressed genes (DEGs) between the 2 groups, followed by analyses of Gene Ontology (GO) functions and Kyoto Encyclopedia of Genes and Genomes pathways. Then, a weighted protein-protein interaction (PPI) network was constructed, and subnetworks and candidate genes were screened. Moreover, the GSE48452 data (14 normal liver tissue samples and 18 nonalcoholic steatohepatitis samples) were used to verify the results obtained from the above analyses.
RESULTS: A total of 100 upregulated genes and 24 downregulated ones were identified in severe NAFLD. Functional enrichment and pathway analyses showed that these DEGs were mainly associated with cell adhesion, inflammatory response, and chemokine activity. The top 5 subnetworks were selected based on the PPI network. A total of 5 hub genes, including ubiquilin 4 (UBQLN4), amyloid-beta precursor protein (APP), sex hormone-binding globulin (SHBG), cadherin-associated protein beta 1 (CTNNB1) and collagen type I alpha 1 (COL1A1), were considered to be candidate genes for NAFLD. In addition, the verification data confirmed the status of COL1A1, SHBG, and APP as candidate genes.
CONCLUSION: UBQLN4, APP, CTNNB1, SHBG, and COL1A1 might be involved in the development of NAFLD, and are proposed as the potential markers for predicting the development of this condition.

Introduction

Nonalcoholic fatty liver disease (NAFLD) is the most common chronic liver condition worldwide.[ Its incidence is as high as 30% in developed countries and nearly 10% in developing countries.[ NAFLD is a clinicopathologic syndrome that encompasses several clinical entities, ranging from simple steatosis to steatohepatitis, fibrosis, and end-stage liver disease.[ Because the etiology and fundamental pathophysiology underlying NAFLD are poorly understood,[ it is challenging to diagnose and treat NAFLD patients before symptomatic cirrhosis or arises.

Gene expression profiling is considered as a powerful tool for exploring diagnostic and predictive biomarkers, especially in the targeting therapy for diseases such as cancer.[ In a previous study, the differences in expression of the transcriptome and proteome in the liver between NAFLD and normal were determined using a gene expression profile, and several key pathways involved xenobiotic and lipid metabolism, inflammatory response, and cell-cycle control were identified.[ A study on the differential gene expression in nonalcoholic steatohepatitis (NASH) also showed that the uptake transporter genes were coordinately targeted for downregulation at the global level during the pathological development of NASH.[ In addition, using microarray analysis, Moylan et al[ indicated that the expression of certain metabolism-related genes was induced in severe NAFLD. Although many genes have been reportedly related to the process of NAFLD development, candidate genes potentially causative of severe NAFLD have not yet been screened, and it has remained unclear whether such genes induce severe NAFLD by interacting with each other.

Against the above background, the present study involved an exploration of the candidate genes of NAFLD, the establishment of a weighted regulatory network and the mining of hub genes in severe NAFLD samples compared with mild NAFLD samples were explored by using bioinformatic methods. We aimed to provide molecular mechanisms underlying NAFLD, and investigate new therapeutic targets of NAFLD.

Materials and methods

Samples

Gene expression profile data GSE49541[ were downloaded from Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geo/) with platform GPL570 (HG-U133_Plus_2) Affymetrix Human Genome U133 Plus 2.0 Array. Data on 72 tissue samples were downloaded, including those from 32 severe NAFLD patients (fibrosis stages 3–4) and 40 mild NAFLD patients (fibrosis stage 0–1). The groups were matched for sex, age (±5 years) and body mass index ( kg/m2) (±3 points). The samples were collected from NAFLD cases in the Duke University Health System NAFLD Biorepository, with approval from the Institutional Review Board of Duke University. Biorepository liver samples are remnants from clinically indicated liver biopsies.

Data preprocessing and differential expression analysis

The normalization of gene expression profile data was performed using the Robust Multichip Averaging (RMA) method[ of the affy package[ in R (v.3.0.0) (http://bioconductor.org/biocLite.R), and the Linear Models for Microarray Data (limma, http://www.bioconductor.org/packages/release/bioc/html/limma.html) package[ was applied to identify the differentially expressed genes (DEGs) by comparing the gene expression levels in samples between mild and severe NAFLD cases. Resampling-based empirical Bayes multiple testing procedures[ were also conducted to correct the P value. Subsequently, an adjusted P value <.05 and |logFC (log fold change)| >0.58 were selected as the thresholds for DEG screening.

Gene Ontology annotation and pathway analysis

The Database for Annotation, Visualization, and Integrated Discovery (DAVID)[ (http://david.abcc.ncifcrf.gov/) is a tool for the functional classification of genes that provides a comprehensive set of functional annotation tools enabling investigators to understand the biological meaning behind large lists of genes. Here, DAVID was used for GO annotation analysis (including the 3 categories of biological process, cellular component, and molecular function 3 aspects) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis. A P value <.05 was considered to indicate a significant difference.

Weighted regulatory network construction

The protein–protein interactions (PPIs) that were related to genes in GSE49541 were selected according to the Human Protein Reference Database (HPRD, (http://www.hprd.org/).[ The average value of rank correlation coefficient and the difference for pairs of regulatory relationship in PPIs were calculated according to the Eq. (1). The average absolute value of the rank correlation coefficient in control samples was considered as the weight of PPI,[ and then the weighted PPI network was constructed.

Here, E is the edge between gene V and gene V, k is the kth sample, V and V are ranked by their expression in the samples, respectively; X is the rank of V of kth sample, X is the rank of V of kth sample, x and x are the average rank s of V and V in the samples, respectively.

Here, r and r represent the Spearman coefficients of r in 2 samples, respectively.

Finally, based on the permutation test, the random of each PPI was calculated. Subsequently, the sample labels were permuted for 10,000 times and a random was generated. The PPIs with a value >90% random value were filtered out.[

Subnetwork investigation and protein–protein interaction score calculation

In the PPI network, the nodes with a degree >15 were defined as candidate genes potentially causative of disease, and the subnetworks consisted of candidate genes and the genes with which they interact. The score of PPIs in the subnetworks were also calculated. Briefly, all PPIs in the weighted network were ranked from large to small according to their weight coefficient and this was defined as the background set (E), whereas the subnetwork was defined as the objective set (S). Then, the enrichment score (ES) of the subnetwork was calculated by walking down background set using the Gene Set Enrichment Analysis method[ (http://www.broadinstitute.org/gsea/index.jsp). The formula is listed as below:

where E is the jth PPI in the ranked regulatory pairs; r is the weight of the jth PPI pair in background set; P is a parameter and set as 1; N is the number of PPI in E; NH is the number of PPI in the subnet S. The ES was equal to the maximum deviation between Phit and Pmiss.

The PPI without contribution to ES was removed from subnetwork.[ To estimate the significance of ES of the subnetwork, pairs of background regulatory relationships were rearranged randomly for 1000 times, the random ES of subnetwork was calculated, while the ES was transformed into Z value based on the equation[:

where ES (bar) is the mean of the random ES set and S′ is the standard deviation of the random ES set.

Calculating the enrichment score of differentially expressed genes in subnetwork

On the basis of the subnetwork obtained above, the ES of DEGs were calculated. Briefly, the genes in gene microarray and subnetwork were considered as background set and objective set, respectively; then, these genes were arranged from large to small; the ES of subnetwork walked by background set was calculated using the following equation[:

Here, g is the jth gene in the ranked genes, r is the magnitude of differential expression of the jth gene, P is a parameter and set as 1, M is the number of genes in L, and MH is the number of genes in Strimmed.

To estimate the significance of ES of the subnetwork, the background genes were rearranged randomly for 1000 times. The random ES of the subnetwork was calculated, followed by the ES being transformed into a Z value based on the Eq. (5).

Candidate genes screening

The 2 Z values obtained as described above were normalized and summed, followed by the acquisition of the combined these 2 parts. Here, using this combined score, the top 5 subnetworks were considered as candidate subnetworks of the disease, whereas the hub genes were defined as candidate genes of the disease.[

Data verification

Gene expression profile data GSE48452 were downloaded from the GEO database (http://www.ncbi.nlm.nih.gov/geo/) with the platform (HuGene-1_1-st) Affymetrix Human Gene 1.1 ST Array [transcript (gene) version], including data on 14 normal liver tissue samples and 18 NASH samples. The liver samples had been obtained from NAFLD patients, who had provided written informed consent. The study protocol was also approved by the institutional review board (“Ethikkommissionder Medizinischen Fakult der University Kiel,” D425/07, A111/99). Empirical Bayes analysis using t test procedures in the limma package (Version 3.10.3, http://www.bioconductor.org/packages/2.9/bioc/html/limma.html) was applied to identify the DEGs by comparing gene expression levels between normal and NASH samples. Subsequently, an adjusted P value <.05 and |logFC| >0.58 were selected as the significance thresholds for DEGs. Subsequently, the Pearson correlation analysis was performed to confirm the differences between the key candidate genes in severe NAFLD and the target DEGs with which they interacted.

Results

Identification of differentially expressed genes

The gene expression profile data were normalized by using the RMA method (Fig. 1). Upon applying the thresholds of an adjusted P value <.05 and |logFC| >0.58, a total of 124 DEGs were identified, including 100 upregulated genes and 24 downregulated ones (Supplementary Table 1). The heatmap of the DEGs is shown in Figure 2. The top 10 upregulated genes and downregulated genes are listed in Table 1.

Figure 1

Box plot of gene expression profile data before and after normalization. The top box represents the distribution of data before normalization, and the bottom box represents the distribution of data after normalization. Horizontal axis represents 72 nonalcoholic fatty liver disease samples.

Figure 2

The expression of differentially expressed genes in nonalcoholic fatty liver disease samples. The bottom horizontal axis represents nonalcoholic fatty liver disease samples, the first 40 samples were from mild nonalcoholic fatty liver disease samples and the other 32 samples were from severe nonalcoholic fatty liver disease. The right vertical axis represents genes. Red indicates upregulated genes and green indicates downregulated genes.

Table 1

The top 10 upregulated and downregulated differentially expressed genes.

Gene Ontology annotation analysis and pathway analysis

GO annotation analysis of the DEGs was performed by using DAVID, the main results of which are listed in Table 2 (Supplementary Table 2). The biological processes that were particularly commonly associated with NAFLD were cell adhesion, chemotaxis, collagen fibril organization, and inflammatory response; the cellular components were extracellular matrix and collagen; and the molecular function was chemokine activity. Furthermore, the identified DEGs were significantly clustered in cell adhesion and chemokine signaling pathways.

Table 2

Gene Ontology annotation analysis of differentially expressed genes in 3 aspects.

Weighted protein–protein interaction network analysis

Based on the HPRD, the PPIs related to genes in the gene microarray were selected, and a PPI network was constructed. The entire network consisted of 8964 nodes and 34,915 edges (Fig. 3).

Figure 3

The top 5 subnetworks with the highest Z scores. Red spots represent differentially expressed genes, and blue or yellow spots represent nondifferentially expressed genes.

The enrichment score of subnetwork and candidate genes

By calculating the weighted PPIs of the subnetwork and the score of DEGs, the significance of ES was estimated with Z value. Then, the top 5 subnetworks with the highest Z values were selected (Fig. 4). The corresponding hub genes were considered as candidate genes, which included ubiquilin 4 (UBQLN4), also known as ataxin-1-interacting protein (A1UP), amyloid-beta precursor protein (APP), sex hormone–binding globulin (SHBG), cadherin-associated protein beta 1 (CTNNB1) and collagen type I alpha 1 (COL1A1) (Table 3).

Figure 4

Pearson correlation analysis between APP and the target-DEGs which they interact. A, the correlation between APP and COL1A2 (R = 0.519, P = .002); B, the correlation between APP and COL4A1 (R = 0.638, P < .001); C, the correlation between APP and CRYAB (R = 0.443, P = .011). R represents Pearson correlation coefficient, R > 0 indicates a positive correlation; R < 0 indicates a negative correlation; if |R| is further away from 0, the correlation is stronger. The P value indicates the significant difference. APP = amyloid-beta precursor protein, DEG = differentially expressed gene.

Table 3

The candidate genes and their interacted genes.

The results of verification data

To confirm the above results, the DEGs between normal and NASH samples were identified through data verification. A total of 181 DEGs were identified, including 119 upregulated DEGs and 62 downregulated DEGs (Supplementary Table 3). As shown in Table 4, COL1A1 and SHBG exhibited significant differences in their expression level. AKR1B10 and CYP2C19, which were included in the original top10 DEGs, were also verified. UBQLN4, APP, and SHBG interacted with the target DEGs in the PPI networks. Pearson correlation analysis revealed that there were significant positive correlations between APP and its target-DEGs (COL1A2, CRYAB, COL4A1) (P < .05, Fig. 4A–C).

Table 4

The results of verification data including the top 10 differentially expressed genes and the 5 candidate genes.

Discussion

NAFLD, is an emerging public health problem, that may be a highly chronic liver condition.[ Although the histologic method has been approved, a useful therapy method for treating NAFLD is yet to be established.[ Gene expression profiling is valuable for researching diagnostic and predictive biomarkers of disease, including NAFLD,[ which may facilitate the development of new therapeutic drugs. In the present study, a total of 124 DEGs associated with severe NAFLD were identified by comparing the gene expression profile with that in mild NAFLD. These DEGs were mainly associated with cell adhesion, inflammatory response, and chemokine signaling pathways. Selection of top 5 subnetworks based on the PPI network indicated that UBQLN4, APP, CTNNB1, SHBG, and COL1A1 may be involved in the development of NAFLD, and these were considered to be markers with potential utility for predicating NAFLD.

In a previous study, inflammation was proved to be associated with the process of NAFLD development.[ Mikolasevic et al[ demonstrated that in patients maintained on hemodialysis, there is probably some interaction between NAFLD and inflammation, malnutrition, and atherosclerosis. In support of this, a significant correlation between the intima-media thickness of the carotid artery and hepatic inflammation score was identified in NAFLD rats.[ Moreover, Browning and Horton[ indicated that the histological hallmarks of NASH, such as inflammation, cell death, and fibrosis promoted the progression of NAFLD. In fact, findings showed that inflammation resulting in a stress response of hepatocytes, might lead to lipid accumulation, and therefore could precede steatosis in NASH.[ In another study based on gene expression in human cases of NAFLD, Greco et al[ asserted that cell adhesion was significantly associated with liver fat content. In severe liver injury, neural cell adhesion molecules weaken the cell–cell and cell–matrix interactions, thereby allowing ductular reactions/hepatic progenitor cells to migrate for normal development and regeneration.[ In fact, each cell adhesion molecule may play an important role during development in hepatic histogenesis, including hepatoblast/hepatocyte-stellate cell interactions.[ In the present study, the selected DEGs were particularly associated with processes such as cell adhesion and inflammatory response. Thus, we speculated that these biological processes such as cell adhesion, fibrosis, and inflammatory response might play important roles in the development of NAFLD.

COL1A1 has been demonstrated to be upregulated in liver fibrosis by the activation of stellate cells and the progression of liver fibrosis.[ Zhao et al[ showed that COL1A1 gene polymorphism was associated with liver fibrogenesis, since the T allele at 1997 of COL1A1 was crucial to the increased transcriptional activity. COL1A2 is an independent predictor of survival in diseases.[ Many diseases that are inherited in an autosomal dominant fashion are caused by mutations in the COL1A1/COL1A2 genes.[ Because COL1A1 and COL1A2 were revealed to be DEGs associated with NAFLD in this study, we speculated that not only COL1A2, but also the dysregulation of COL1A1/COL1A2 participated in the progression of NAFLD. Furthermore, the present study based on GO and KEGG pathway analyses showed that COL1A1 might be related to cell adhesion and inflammation; thus, COL1A1/COL1A2 might play important roles in diseases by regulating genes associated with cell adhesion and inflammation. In addition, COL1A1 was confirmed to be upregulated in NASH samples compared with the level in normal controls. SHBG, a glycoprotein expressed predominantly in the hepatocytes, regulates the transport of sex steroid hormones in the bloodstream to their target tissues.[ As one of the circulation factors released from fatty liver, SHBG was reported to be directly involved in the pathogenesis of local and systemic inflammation, and peripheral as well as hepatic insulin resistance.[ In the present study, APP and SHBG were shown to be connected by COL1A2, which further indicates that the APP and SHBG genes have a close relationship in the process of NAFLD development.

Ubiquilins (UBQLN), a family of ubiquitin-binding proteins, are involved in several protein degradation pathways and have been implicated in various diseases.[ Previous studies indicated that the members of this family-mediated degradation of misfolded proteins and they were implicated in a number of pathological and physiological conditions.[ For example, Matsuda et al[ indicated that UBQLN4 was highly expressed in organs such as liver. Unfortunately, to date, few details of the relationship between UBQLN4 and NAFLD have been clarified. In this study, UBQLN4 was revealed to be a hub gene in the PPI network. Thus, we speculated that UBQLN4 might participate in the progression of via the degradation of misfolded proteins. CTNNB1 is located on the short arm of chromosome 3, as determined by in situ fluorescence analysis,[ and has been reported to be commonly involved in benign liver tumorigenesis.[ Mutation of the CTNNB1 gene mutation is also an important indicator of prognosis in primary sporadic aggressive fibromatosis.[ In addition, CTNNB1 was also found to often display point mutations resulting in loss of function in a range of cancers, with the notable exception in hepatocellular carcinoma,[ revealing that the expression level of CTNNB1 might have great significance for liver function. Kubota et al[ indicated that mutational analysis of CTNNB1 for predicting disease such as solid-pseudopapillary neoplasm is feasible. In the present study, CTNNB1 investigated given its status here as a hub protein. The findings suggest that CTNNB1 is useful as a target gene for further investigation of NAFLD. The analyses demonstrated the detailed action of these 2 candidate genes in NAFLD, and also showed that the interaction between them is rare. Therefore, further experiments are needed to explain the functions of CTNNB1 and UBQLN4 in NAFLD.

Despite all these results, there are some limitations in the present study. First, verification experiments were not performed because of the lack of sufficient liver tissue samples. Second, we did not obtain the gene expression data, including for normal liver tissue, severe NAFLD tissue, and mild NAFLD tissue, categorized according to fibrosis stage in the NCBI database. Therefore, more experiments and data are required to verify the above results.

In conclusion, cell adhesion, fibrosis, and inflammatory response might play crucial roles in severe NAFLD. Similarly, 5 candidate genes, namely UBQLN4, APP, CTNNB1, SHBG, and COL1A1 might be involved in the development of NAFLD, and might be considered as potential markers for predicting the development of NAFLD. The findings of this study might explain the development of NAFLD and reveal new target genes for treating NAFLD.