Screening the key genes of hepatocellular adenoma via microarray analysis of DNA expression and methylation profiles.

The aim of the present study was to identify the biomarkers involved in the development of hepatocellular adenoma (HCA) through integrated analysis of gene expression and methylation microarray. The microarray dataset GSE7473, containing HNF1α-mutated HCA and their corresponding non-tumor livers, 5 HNF1α-mutated HCA and 4 non-related non-tumor livers, was downloaded from the Gene Expression Omnibus (GEO) database. The DNA methylation profile GSE43091, consisting of 50 HCA and 4 normal liver tissues, was also downloaded from the GEO database. Differentially expressed genes (DEGs) were identified by the limma package of R. A t-test was conducted on the differentially methylated sites. Functional enrichment analysis of DEGs was performed through the Database for Annotation, Visualization and Integrated Analysis. The genes corresponding to the differentially methylated sites were obtained by the annotation files of methylation chip platform. A total of 182 DEGs and 3,902 differentially methylated sites were identified in HCA. In addition, 238 enriched GO terms, including organic acid metabolic process and carboxylic acid metabolic process, and 14 KEGG pathways, including chemical carcinogenesis, were identified. Furthermore, 12 DEGs were identified to contain differentially methylated sites, among which, 8 overlapped genes, including pregnancy zone protein and solute carrier family 22 member 1 (SLC22A1), exhibited inverse associations between gene expression levels and DNA methylation levels. The DNA methylation levels may be potential targets of HCA. The present study revealed that the 8 overlapped genes, including annexin A2, chitinase 3-like 1, fibroblast growth factor receptor 4, mal, T-cell differentiation protein like, palladin, cytoskeletal associated protein, plasmalemma vesicle associated protein and SLC22A1, may be potential therapeutic targets of HCA.

Introduction

Hepatocellular adenoma (HCA) is a benign liver tumor that occurs mainly in young females subsequent to the long-term use of oral contraceptives (1,2). It is a type of rare tumor with low morbidity (3). However, HCA may lead to hemorrhage, and even malignant transformation to hepatocellular carcinoma (HCC), which is the third leading cause of cancer-associated mortality worldwide (4–8). A rising incidence has been reported due to improved application of diagnostic imaging techniques (9). HCA is rare in children, men and post-menopausal women. The use of androgenic steroids for Fanconi's anemia or acquired aplastic anemia is a risk factor for the development of HCA (10). Other drugs are involved in its development, such as clomiphene, barbiturates and recombinant human growth hormone (11–13). Obesity and alcohol abuse have also been reported as risk factors for developing HCA (14).

Previously, ultrasound and magnetic resonance imaging (MRI) have been suggested as effective tools for the diagnosis of HCA and explore its biological mechanisms (15). According to the different genetic mutations, HCA is divided into 4 major molecular subgroups: Hepatocellular nuclear factor-1α (HNF1α)-mutated type HCA; β-catenin-mutated type HCA; inflammatory type HCA; and unclassified type HCA (8,14,16). The transcription factor 1 (TCF1) gene, which encodes HNF1α, has been identified in the liver (16). Overexpression of certain genes, including erb-b2 receptor tyrosine kinase 2, mechanistic target of rapamycin, platelet-derived growth factor α polypeptide, platelet-derived growth factor β polypeptide and cyclin D1, has been identified in HCA, and the products of these genes are associated with cell proliferation, cell cycle activation and angiogenesis (17).

In the human genome, GC-rich DNA sequences, also known as CpG islands, are frequently enriched in the first exon and the promoter (18). DNA methylation at CpG islands located upstream of a gene promoter is associated with differential expression of the gene. DNA methylation regulates gene silencing by directly inhibiting the binding of methylation-dependent transcriptional activators or indirectly altering the affinity of proteins, including methylated DNA binding domain protein, involved in chromatin remodeling (19–23). At present, DNA methylation is widely identified in human cancer, including HCC. It was reported that long interspersed nuclear element-1 (LINE-1) has lower DNA methylation levels in hepatitis virus and aflatoxin-associated HCC compared with normal liver tissue (24,25). In addition, hypomethylation of LINE-1 was associated with advanced disease and poorer survival in HCC (26). The DNA methylation level of spermidine/spermine N1-acetyltransferase family member 2 (SAT2) also has a significant role in liver carcinogenesis. It has been suggested that decreased SAT2 methylation of white blood cell DNA was significantly associated with increased HCC risk later in life (27).

It is challenging to diagnose HCC and HCA at an early stage. Numerous therapies are limited when HCC enters the advanced stage, as the advanced stage is accompanied by severe liver dysfunction (28). Therefore, it is necessary to identify early biomarkers of HCA. In the present study, data of the DNA methylation profile and gene expression profile was extracted from the Gene Expression Omnibus (GEO) database. Certain therapeutic targets and the related pathways that may be associated with the development of HCA were identified by microarray analysis. This may contribute to promoting available biomarkers for the early diagnosis, therapy and prognosis of HCA.

Materials and methods

Microarray data

The gene expression profile and DNA methylation profile were both downloaded from the GEO (http://www.ncbi.nlm.nih.gov/geo/) database. The gene expression profile (GSE7473) contained 41 samples, including 8 HNF1α-mutated HCA and the corresponding non-tumor liver samples (each sample was assessed four times using 11K_VJF-ARRAY; GPL3282), and 5 HNF1α-mutated HCA and 4 non-related non-tumor liver samples (the 9 samples were assessed using GPL96 Affymetrix Human Genome U133A Array). In the present study, the 9 samples that were assessed via GPL96 were used as the objects, and the 5 HNF1α-mutated HCA and 4 non-related non-tumor liver samples were classified as the case and control groups, respectively. The DNA methylation profile (GSE43091), provided by Pilati et al (29), contained 50 HCA and 4 normal liver tissues. These 54 samples were achieved by GPL13534 Illumina HumanMethylation450 BeadChip (HumanMethylation450_15017482).

Data preprocessing

The raw microarray data were converted into expression data using the affy package of R. The values of multiple probes that correspond to the same gene were summarized. For original DNA methylation data, the β value of every methylated site was calculated and normalized using the IMA package of R.

Identification of differentially methylated sites and differentially expressed genes (DEGs)

DEGs were identified using the limma package of R with P<0.05 and |log2 (fold-change)|>1. A paired Student's t-test was conducted on the methylation levels between HCA samples and normal samples, and the differentially methylated sites with adjusted P<0.05 and |Δβ|>0.2 were selected.

Functional enrichment analysis of DEGs

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis of DEGs was performed using the Database for Annotation, Visualization and Integrated Discovery (DAVID). DAVID was used to perform functional annotation for a list of genes, gene functional classification or gene ID conversion. All GO terms and KEGG pathways with P<0.05 that contained at least five genes were selected for subsequent analysis.

Comprehensive analysis of gene expression profile and DNA methylation profile

The genes in which differentially methylated sites were located were identified using the annotation files of the methylation chip platform.

Results

Differentially methylated sites and differentially expressed genes

In total, 182 DEGs (53 upregulated and 129 downregulated) were identified in HCA. The volcano plot (Fig. 1) showed the distribution of DEGs. From the heatmap (Fig. 2), the case samples were found to be distinguished from the control samples. Additionally, a total of 3,902 differentially methylated sites were obtained, including 3,715 downregulated methylated sites and 187 upregulated methylated sites. These methylated sites were mostly located in the intergenic and gene-coding regions of genes (Fig. 3).

Figure 1.

A volcano plot of differentially expressed genes. The red plus signs represent upregulated differentially expressed genes, the blue triangles represent downregulated differentially expressed genes, and the black circles represent non-differentially expressed genes. FC, fold-change

Figure 2.

Two-way clustering analysis for metastasis in 5 hepatocellular nuclear factor-1α-mutated hepatocellular adenoma and 4 non-related non-tumor livers. Expression level was normalized per gene, and the relative value to the median among nine samples is shown by color. Red and black-red indicate high expression, and green represents relatively low expression.

Figure 3.

Location distribution of differentially methylated sites in genes. Opensea indicates the intergenic region; body indicates the gene coding region; 5′UTR indicates the 5′UTR region; 3′UTR indicates the 3′UTR region; 1st exon indicates the first expressed region; TSS200 indicates the 200 bp upstream of the transcription start site; and TSS1500 indicates the 1,500 bp upstream of the transcription start site. UTR, untranslated region.

Enriched GO terms and KEGG pathways

In the present study, a total of 238 enriched GO terms and 14 KEGG pathways were identified according to the criteria P<0.05. The top 20 enriched GO terms are listed in Table I. The majority of the enriched GO terms were involved in the organic acid metabolic process. The enriched KEGG pathways of the DEGs are shown in Table II. Certain KEGG pathways, for example the chemical carcinogenesis pathway, mineral absorption pathway and Bile secretion, were directly associated with HCA, and they may affect the development of HCA.

Table I.

Top 20 enriched GO terms for differentially expressed genes.

Category	GO ID	GO name	Gene number	P-value
BP	GO:0006082	Organic acid metabolic process	54	1.89×10−23
BP	GO:0019752	Carboxylic acid metabolic process	51	2.08×10−23
BP	GO:0043436	Oxoacid metabolic process	51	2.88×10−21
BP	GO:0032787	Monocarboxylic acid metabolic process	36	3.97×10−20
BP	GO:0042493	Response to drug	26	4.53×10−14
BP	GO:0006805	Xenobiotic metabolic process	17	2.54×10−13
BP	GO:0071466	Cellular response to xenobiotic stimulus	17	2.82×10−13
BP	GO:0010038	Response to metal ion	21	3.60×10−13
BP	GO:0009410	Response to xenobiotic stimulus	17	5.16×10−13
BP	GO:0055114	Oxidation-reduction process	38	1.05×10−12
BP	GO:0010035	Response to inorganic substance	24	1.90×10−12
MF	GO:0016491	Oxidoreductase activity	31	2.60×10−12
BP	GO:0006629	Lipid metabolic process	40	1.54×10−11
MF	GO:0004497	Monooxygenase activity	13	2.02×10−11
BP	GO:0071294	Cellular response to zinc ion	7	4.12×10−11
CC	GO:0005615	Extracellular space	39	4.44×10−11
BP	GO:0044282	Small molecule catabolic process	19	4.76×10−11
BP	GO:0008202	Steroid metabolic process	19	1.20×10−10
BP	GO:0071248	Cellular response to metal ion	12	2.21×10−10
MF	GO:0020037	Heme binding	13	4.63×10−10

GO, Gene Ontology; CC, cellular component; MF, molecular function; BP, biological process.

Table II.

Enriched KEGG pathways for differentially expressed genes.

KEGG pathway name	Gene number	P-value
Chemical carcinogenesis	11	7.14×10−09
Mineral absorption	9	1.84×10−08
Bile secretion	10	3.21×10−08
Tryptophan metabolism	8	4.18×10−08
Linoleic acid metabolism	7	7.54×10−08
Retinol metabolism	9	1.44×10−07
Drug metabolism-cytochrome P450	9	2.47×10−07
Arginine and proline metabolism	8	7.43×10−07
Steroid hormone biosynthesis	8	7.43×10−07
Metabolism of xenobiotics by
cytochrome P450	8	5.61×10−06
Serotonergic synapse	8	1.31×10−04
Arachidonic acid metabolism	6	1.94×10−04
Carbon metabolism	6	1.44×10−03
Glycolysis/gluconeogenesis	5	1.78×10−03

KEGG, Kyoto Encyclopedia of Genes and Genomes.

Key genes in hepatocellular adenoma

In total, 12 DEGs that contained differentially methylated sites were identified in the case groups compared with the control groups (Fig. 4). Among the DEGs, 8 genes with inverse associations between gene expression level and DNA methylation level were identified, consisting of annexin A2 (ANXA2), chitinase 3-like 1 (CHI3L1), fibroblast growth factor receptor 4 (FGFR4), mal, T-cell differentiation protein like (MALL), palladin, cytoskeletal associated protein (PALLD), plasmalemma vesicle associated protein (PLVAP), pregnancy zone protein (PZP) and solute carrier family 22 member 1 (SLC22A1).

Figure 4.

In total, 12 genes were differentially expressed and methylated in the case groups compared with the control groups. LogFC_Gene, fold-change in gene differential expression; Beta_differ, fold-change in gene methylation level; ANXA2, annexin A2; CHI3L1, chitinase 3-like 1; FGFR4, fibroblast growth factor receptor 4; MALL, mal, T-cell differentiation protein like; PALLD, palladin, cytoskeletal associated protein; PLVAP, plasmalemma vesicle associated protein; PZP, pregnancy zone protein; SLC22A1, solute carrier family 22 member 1.

Discussion

The development of human cancers is associated with two factors: Gradual accumulation and mutual interactions of genetic and epigenetic alterations (30). As one of the major characteristics in human cancers, epigenetic alterations may suggest the molecular mechanisms underlying malignant transformation (30–33). DNA methylation is one of these epigenetic mechanisms, and is widely found in human cancers, including HCA. At present, studies have reported that a family of DNA methyltransferase enzymes (DNMTs) mediates DNA methylation. DNMT1 was found to maintain methylation, whereas DNMT3A and DNMT3B induced de novo methylation (34). Promoter hypermethylation was associated with gene expression and resulted in transcriptional inhibition and loss of gene function. Certain studies have revealed that dysregulation of the removal and establishment of DNA methylation was involved in hepatocarcinogenesis (35,36). In HCC, hypermethylation mainly affects the expression of certain tumor suppressor genes, particularly the genes involving in cell differentiation, cell proliferation, cell adhesion, cellular metabolism, and DNA repair. Hypermethylated genes, including adenomatosis polyposis coli, Ras association domain family member 1 and suppressor of cytokine signaling 1, have been identified in chronic hepatitis and cirrhosis (20,37–40). Also, genes such as glutathione S-transferase pi 1, cyclin dependent kinase inhibitor 2A, cytochrome c oxidase subunit II, HIC ZBTB transcriptional repressor 1, and runt related transcription factor 3, which were frequently methylated, have been identified in dysplastic liver nodules (20,41–43).

Previously, the emergence of new diagnostic methods such as ultrasound, and novel treatment methods such as liver-directed therapy, has improved the prognosis of HCA (44,45). However, a limit remains in terms of its early diagnosis and curative potential. The present study used microarray technology to identify key factors, such as the genes, biological process and signal pathways, involved in HCA.

In the present study, DNA expression and methylation profiles were obtained by bioinformatics to identify the differentially methylated sites and DEGs in HCA compared with normal liver tissues. A total of 182 DEGs (53 upregulated and 129 downregulated) and 3,902 differentially methylated sites (187 upregulated and 3,715 downregulated) were identified. In addition, 8 overlapped genes with inverse correlations between methylation levels and gene expression levels were identified, including PZP and SLC22A1.

Important genes may be potential targets for HCA diagnosis or treatment. PZP is a major pregnancy-associated plasma protein that is strongly associated with α 2-macroglobulin (46). PZP is considered an auxiliary index for the identification of gynecological tumors. Berne (47) revealed that estrogen induced PZP expression, and PZP was found in the serum of women who usually took oral contraceptives containing estrogen. As the long-term use of oral contraceptives may lead to HCA, this indicates that PZP may be a major gene in HCA. Polyspecific organic cation transporters are involved in the uptake of positively charged and neutral small molecules, certain drugs and environmental toxins into organs including the liver, kidney and intestine (48). In 1967, the human organic cation transporter 1, which is encoded by SLC22A1, was first reported to remove >70% of the serotonin in the portal blood via filtration and metabolism in the liver (49). Solute carrier family 22 (organic cation transporter), member 1 is the major active influx protein responsible for the transport of imatinib mesylate into cells (50). The present study hypothesized that SLC22A1 may play a critical role in removing endogenous substances, drugs and other toxins associated with HCA.

According to the functional enrichment of DEG analysis, certain GO terms associated with metabolic process and response to drugs were identified, and the KEGG pathways closely associated with chemical carcinogenesis and mineral absorption were obtained. Chemical carcinogenesis has become synonymous with genotoxic events, which lead to DNA damage and genetic mutations (51). In addition, epigenetic effects, such as aberrant DNA methylation, have been identified as one of the key contributors to carcinogenesis. Aberrant DNA methylation has been widely studied in carcinogenesis (19,52,53). Although certain animal models, which reflect the association between exposure to chemical carcinogens and epigenetic effects, have been established, the specific mechanism of this process has yet to be clarified (54).

Overall, using bioinformatics analysis, DEGs, differentially methylated sites, significant GO terms and KEGG pathways were obtained. It was found that DEGs were mainly involved in acid metabolic processing and chemical carcinogenesis in HCA. Based on comprehensive bioinformatics analysis, 8 important DEGs, consisting of SLC22A1, PZP, ANXA2, CHI3L1, FGFR4, MALL, PALLD and PLVAP, were identified. These DEGs, which are related to HCA, may potentially act as biomarkers for detection, prognosis, monitoring and predicting therapeutic responses in HCA. However, additional experiments are required to confirm their function in HCA.