This article is mentioned in:

Introduction

Hepatocellular carcinoma (HCC) is a fatal disease (1). Several types of potentially curative treatments have become available due to the advances in technology and surgical techniques (2). Despite these advances in the clinical treatment of HCC, the mortality rate remains high (3). Thus, novel potential therapies for HCC are urgently required.

As potential therapy targets in malignant cells microRNA (miRNA/miR) perform important roles in regulating gene expression (4). Dysregulation of miRNA is widely involved in human cancer, including HCC (5). A previous study indicated that the mutation of the miR-122 binding site in the interleukin-1α 3′ untranslated region is associated with the incidence rate of HCC (6). As a potential tumor suppressor miRNA, miR-1285-3p has been shown to regulate the expression of JUN (Jun Proto-Oncogene, AP-1 Transcription Factor Subunit) in the disease progression of HCC (7). The regulation of feedback between miRNAs and target genes performs important roles in the progression of HCC (8). Downregulated miR-148b has been shown to be a biomarker for the early detection of HCC (9). Previous miRNA-based HCC gene therapy studies have demonstrated significant inhibition of miRNA in HCC cells, indicating a promising alternative to current therapeutics (10–12). However, the association between miRNA and their target genes, as well as the influence of these associations during the process of HCC is rarely studied.

miRNA expression profiling has proven useful in diagnosing and understanding the development and progression of several diseases including HCC (13–15). The use of miRNA expression profiling as a prognostic biomarker for the detection of HCC is practicable (16). Since chronic hepatitis B virus (HBV) infection is a significant cause of HCC, the present study performed a bioinformatics analysis based on the miRNA expression profile of the HBV infected HCC tissue samples and chronic Hepatitis B (CHB) tissue samples with no fibrosis from 12 patients with HCC. The present study may help to explore the potential disease progression of HCC resulting from CHB and provide information regarding HCC-related miRNA/marker genes for the gene therapy of HCC.

Materials and methods

Samples and microarray data

Gene expression profile data GSE67882 was downloaded from the Gene Expression Omnibus (GEO) database (17) (http://www.ncbi.nlm.nih.gov/geo/) based on the platform GPL10850; Agilent-021827 Human miRNA Microarray (V3; Agilent Technologies, Inc., Santa Clara, CA, USA). Since the study was performed based on the raw microarray data retrieved from the GEO database, all differential expressed miRNAs meeting the threshold of adjusted P<0.05 and |log2FC (fold change)|>0.52 were analyzed. A total of 14 upregulated miRNAs and 16 downregulated miRNAs were obtained. The present study was performed by using data from 4 HBV infected HCC tissue samples (HCC group) and 8 CHB tissue samples with no fibrosis (control group) of GSE67882 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE67882). The microarray data was downloaded on July 24, 2015, and the bioinformatic analysis was completed on August 12, 2015.

Data preprocessing and differential expression analysis

The normalization (robust estimation of variance-stabilizing and calibrating transformations) of gene expression profile data was performed using Vsn (18) (http://bioconductor.org/help/search/index.html?q=vsn/) of Bioconductor in R (v.3.0.0). The Linear Models for Microarray Data (limma, http://www.bioconductor.org/packages/release/bioc/html/limma.html) package (19) in Bioconductor was applied to identify differentially expressed miRNAs by comparing the expression levels between samples in the HCC group and the control group. The corresponding P-value of miRNAs subsequent to an unpaired t-test was defined as the adjusted P-value. Subsequently, the adjusted P<0.05 and |log2FC|>0.52 were selected as the threshold values for screening.

Analysis for the target genes regulated by differentially expressed miRNAs

Since the inhibition of gene expression following transcription is realized by the combination of miRNA and mRNA, differentially expressed miRNAs in HCC samples were further analyzed to explore the significance of miRNA in liver cancer. In the present study, two experimental certificate databases including miRecords (20) and miRWalk (21) were used to explore suitable miRNA-target gene relations. The relations that existed in at least one kind of database aforementioned were selected, then the target genes of differentially expressed miRNAs identified by miRecords or miRWalk were obtained.

Gene Ontology (GO) annotation and pathway analysis

The Database for Annotation, Visualization and Integrated Discovery (DAVID; David Bioinformatics Resource) (22) is a gene functional classification tool that provides a comprehensive set of functional annotation tools for investigators to understand the biological meaning behind large lists of genes. GO functional enrichment analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways (http://www.genome.jp/kegg/pathway.html) enrichment analysis (23) of DEGs were performed using DAVID. GO has 3 functional categories: Molecular function (MF); biological process (BP); and cellular component (CC). P<0.05 was considered to indicate a statistically significant difference.

Transcription factor-miRNA association investigation

Each miRNA may have several target genes, but outstanding target genes are more likely to assemble transcription factors (24). Thus, the analysis between transcription factors and miRNA is beneficial to understand this pathomechanism. iRegulon can be used to detect transcription factors, motifs and their optimal sets of direct targets from a set of genes (25). In the present study, the iRegulon in Cytoscape software was used to explore the potential associations between miRNA and transcription factors. The minimum identity between orthologous genes was 0.05, while the maximum false discovery rate on motif similarity was 0.001. The normalized enrichment score (NES) of >4 was considered as the threshold value for the selection of potential associations.

HCC-associated miRNA investigation

The Comparative Toxicogenomics Database (CTD) is a robust, publicly available database that is used to advance understanding regarding how environmental exposures affect human health (26). In the present study, the investigation of HCC-associated miRNA was performed using the CTD database. The biomarkers of a disease are generally identified by gene expression traits, abnormalities in gene structure or mutation which could perform roles in the etiology of a disease. Thus miRNAs which regulate these genes may be important for HCC.

Results

Investigation into differentially expressed miRNAs

The original data were analyzed and filtered. A total of 14 upregulated miRNAs and 16 downregulated miRNAs were obtained with thresholds of P<0.05 and |log2FC|>0.52. The heat map of differentially expressed miRNAs is displayed in Fig. 1.

Figure 1. Heat map for the differentially expressed miRNAs. Green represents the low expression level of miRNA, red represents the high expression level of miRNA and blank represents that the expression levels of miRNAs do not change. HCC, hepatocellular carcinoma; hsa, human; miRNA/miR, microRNA.

Analysis of differentially expressed miRNAs and their associated target genes

Based on two experimental certificate databases (miRecords and miRWalk), a total of 380 suitable miRNA-target gene relations were obtained. A total of 22 miRNAs including hsa-miR-15a, hsa-miR-125b, hsa-miR-122, hsa-miR-146b-5p and hsa-miR-142-3p were revealed to have experimental certificate target genes.

GO and KEGG pathway analysis were used to gain further insights into the pathways significantly involved by miRNAs. GO functional enrichment analysis demonstrated that differentially expressed miRNAs, including hsa-miR-15a and hsa-miR-125b, were mainly involved in the positive regulation of histone H3-K9 methylation (GO, 0051574; BP), fibrillar center (GO, 0001650; CC) and mismatched DNA binding (GO, 0030983; MF; Table I). KEGG pathway analysis in the present study demonstrated that differentially expressed miRNAs were mainly enriched in the pathways in cancer map05200 (http://www.genome.jp/dbget-bin/www_bget?pathway:map05200; Fig. 2). Based on the experimental certificate miRNA-target gene databases aforementioned, the miRNA-target gene regulatory network was constructed using Cytoscape software. The results demonstrate that there were 13 significantly upregulated miRNAs and 9 significantly downregulated miRNAs (Fig. 3).

Figure 2. Results of KEGG pathway enrichment analysis. X-axis represents the number of miRNAs, Y-axis represents the KEGG pathways significantly enriched by differentially expressed miRNAs, P<0.05. KEGG, Kyoto Encyclopedia of Genes and Genomes; miRNA/miR, microRNA; MAPK, mitogen activated protein kinase; ECM, extra cellular matrix.

Figure 3. Regulatory network for miRNA and target genes. Green squares represent the upregulated miRNAs, red squares represent the downregulated miRNAs and blue circles represent the target genes. hsa, human; miRNA/miR, microRNA.

Table I. Results of GO functional enrichment analysis of differentially expressed miRNAs.

Table I.

Results of GO functional enrichment analysis of differentially expressed miRNAs.

GO ID	Term	Count	P-value
BP
0051574	Positive regulation of histone H3-K9 methylation	13	6.04×10−7
0051570	Regulation of histone H3-K9 methylation	13	6.48×10−6
0031061	Negative regulation of histone methylation	12	3.76×10−5
0031062	Positive regulation of histone methylation	13	5.34×10−5
0010001	Glial cell differentiation	21	7.58×10−5
CC
0001650	Fibrillar center	6	1.32×10−4
0035098	ESC/E(Z) complex	12	8.27×10−4
0005869	Dynactin complex	5	1.09×10−3
0044452	Nucleolar part	7	4.18×10−3
0000791	Euchromatin	8	8.40×10−3
MF
0030983	Mismatched DNA binding	7	1.84×10−4
0008263	Pyrimidine-specific mismatch base pair	6	2.74×10−4
	DNA N-glycosylase activity
0043398	HLH domain binding	9	3.92×10−4
0035198	miRNA binding	13	4.56×10−4
0000700	Mismatch base pair DNA N-glycosylase activity	6	6.17×10−4

[i] The top 5 in BP, CC and MF are listed. P<0.05 was considered to indicate a statistically significant difference. GO, Gene Ontology; BP, biological process; CC, cellular component; MF, molecular function.

Transcription factor-miRNA relations analysis

The iRegulon in Cytoscape software was used to explore the potential association between miRNA and transcription factors. With NES >4, the 3 miRNAs with the highest number of target genes were hsa-miR-15a (Fig. 4), hsa-miR-125b (Fig. 5) and hsa-miR-122 (Fig. 6). There were 61 target genes that were upregulated by hsa-miR-15a, including transcription factors such as zinc finger protein 143, poly (ADP-ribose) polymerase 1, V-maf musculoaponeurotic fibrosarcoma oncogene homolog K, DNA polymerase epsilon subunit 3, heat shock factor protein 2 and interferon regulatory factor 8. A total of 54 target genes were downregulated by hsa-miR-125b, including transcription factors such as general transcription factor II–I repeat domain-containing protein 1, serum response factor, polymerase (RNA) II (DNA directed) polypeptide A, B-cell lymphoma 6 protein and transcription factor E2F1. A total of 49 target genes were downregulated by hsa-miR-122, including transcription factors such as forkhead box protein A2, CCAAT/enhancer-binding protein α and estrogen-associated receptor γ.

Figure 4. Sub-modules of regulatory networks for hsa-miR-15a and associated transcription factors. Yellow circles represent the transcription factors, purple circles represent the co-regulated target genes for transcription factors and blue circles represent the target genes of miR-15a. hsa, human; miRNA/miR, microRNA.

Figure 5. Sub-modules of regulatory network for hsa-miR-125b and associated transcription factors. Yellow circles represent the transcription factors, purple circles represent the co-regulated target genes for transcription factors and blue circles represent the target genes of miR-125b. hsa, human. miRNA/miR, microRNA.

Figure 6. Sub-modules of regulatory network for hsa-miR-122 and associated transcription factors. Yellow circles represent the transcription factors, purple circles represent the co-regulated target genes for transcription factors and blue circles represent the target genes of miR-122. hsa, human; miRNA/miR, microRNA.

HCC-associated miRNA analysis

A total of 11 differentially expressed miRNA that regulated the marker genes of HCC were identified (Table II). Among them, hsa-miR-146-5p (n=13), hsa-miR-15a (n=5) and hsa-miR-125b (n=4) were the miRNAs that regulated the highest number of marker genes of HCC. A further investigation was made on other HCC-associated genes (not the marker genes but recorded in databases) in the CTD database. The results indicated that a total of 49 target genes of miR-122 were HCC-associated genes.

Table II. Differentially expressed miRNAs that regulated the marker genes of hepatic cell carcinoma (HCC).

Table II.

Differentially expressed miRNAs that regulated the marker genes of hepatic cell carcinoma (HCC).

miRNA	Marker gene
hsa-miR-99a	MTOR, IGF1R
hsa-miR-424	KDR
hsa-miR-34a	CCND1, MYC, MET
hsa-miR-23b	MET
hsa-miR-21	PPARA, E2F1, BTG2
hsa-miR-199a-5p	CCND1, MET, FASN
hsa-miR-15a	UGDH, BCL2, CCND1, HSDL2, JUN
hsa-miR-148a	CDKN1B
hsa-miR-146b-5p	JUN, IL6, FOS, MET, TNF, SOCS1, CCND1, MYC, BRAF, MMP9, IGF1R, BCL2, FASN
hsa-miR-142-3p	APCS, KRAS, JUN, CCND1
hsa-miR-125b	ADAMTS1, TP53, CDKN2A, CYP1A1

[i] hsa, human; miR/miRNA, microRNA.

Discussion

HCC is one of the most common types of human cancer worldwide (1). However, the molecular mechanisms behind the pathogenesis and progression of HCC remain unclear. The present study performed a bioinformatics analysis on the miRNA expression profile of HCC. By screening the differentially expressed miRNAs, a total of 14 upregulated miRNAs and 16 downregulated miRNAs between HCC samples and control samples were obtained. These differentially expressed miRNAs were mainly involved in BP like positive regulation of histone H3-K9 methylation, and were observed in KEGG pathways in cancer map05200. A total of 3 outstanding regulatory networks of miRNA (hsa-miR-15a, hsa-miR-125b and hsa-miR-122) were further investigated. A total of 11 differentially expressed miRNAs that regulated the candidate marker genes of HCC were explored.

In the present study, GO enrichment analysis demonstrated that the target genes that were regulated by the differentially expressed miRNA, such as hsa-miR-15a and hsa-miR-125b, were mainly enriched in functions such as histone H3-K9 methylation. DNA methylation and histone H3-K9 modifications perform important roles in the process of tumor formation by the dysregulation of target genes (27). A previous study demonstrated that H3-K9 methylation contributed to chromosomal targeting of polycomb group proteins (28). DNA methylation profiling for serum DNA of patients with HCC revealed a potential association between HCC and DNA methylation (29). Peng et al (30) indicated that heterochromatin (essential for chromosome organization and inheritance) stability requires regulators of histone H3-K9 methylation, which indicates that regulation of histone H3-K9 methylation is vital for cancer-associated gene expression. Thus, the present study speculated that histone H3-K9 modifications may relate to the progression of HCC. KEGG pathway analysis in the present study revealed that target genes of differentially expressed miRNAs were mainly enriched in multiple cellular processes related to cancer and MAPK signaling. Daroqui et al (31) proved that MAPK signaling promotes cell invasiveness and in vivo tumor progression via actin remodeling in response to transforming growth factor-β. Tumor necrosis factor has been revealed to activate the cyclic adenosine monophosphate response element-binding protein through the p38 MAPK/mitogen and stress activated protein kinase 1 signaling pathway (32). MAPK signaling and pathways in cancer had a close association with HCC processes, and so the differentially expressed miRNA that relate with the MAPK signaling pathway and cancer-related pathways (as determined by KEGG analysis) may take part in the development of HCC.

It is considered that the regulation of transcription factors and miRNAs is essential for the appropriate execution of biological events and developmental processes (33). By incorporating prior information from predicted regulatory interactions with parallel miRNA/mRNA expression datasets, Peng et al (34) indicated that the transcription factor-miRNA co-regulatory based on feed-forward loops served an important role in the development of prostate cancer. In the present study, the regulatory networks of 3 outstanding miRNAs (hsa-miR-15a, hsa-miR-125b and hsa-miR-122) were obtained. This result demonstrates that these transcription factors had close associations with the 3 miRNAs in the co-regulation of target genes, which further indicated the complex regulatory roles of them during the development of HCC. Thus, it is speculated that hsa-miR-15a, hsa-miR-125b and hsa-miR-122 may be the miRNAs to be investigated for their role in HCC gene therapy.

Human miR-146b-5p is an miRNA which is located on chromosome 10q24.3 (35). miR-146b-5p is frequently downregulated in solid tumors, including prostate cancer, pancreatic cancer and glioblastoma (36). In the present study, HCC-associated miRNA analysis revealed that there were 11 differentially expressed miRNAs that regulate the marker genes of HCC including miR-146b-5p, which had the highest number of marker genes. Furthermore, the HCC-associated gene (not recorded marker genes) investigation based on the CTD database revealed that miR-122 had the highest number of HCC-associated genes. A previous study indicated that miR-122 negatively correlates with liver fibrosis as detected by histology and FibroScan software (37). In animal models, miR-122 has been proven to be a marker of liver cell damage, inflammation and regeneration in acetaminophen-induced liver injury (38). In the present study, miR-122 was revealed to be closely associated with HCC, thus the target gene of miR-122 may be the potential novel marker gene for HCC. Therefore, it was speculated that hsa-miR-15a, hsa-miR-125b, hsa-miR-122 and hsa-miR-146b-5p may perform important roles in the progression of HCC, and may be new biomarkers for gene therapy of HCC. However, there were limitations to the present this study, such as lack of experimental verification and joint analysis of expression profile data. An investigation based on experimental verification is required in future studies.

In conclusion, miRNAs such as hsa-miR-15a, hsa-miR-125b, hsa-miR-122 and hsa-miR-146b-5p may be new biomarkers for gene therapy of HCC. Histone H3-K9 methylation and miRNAs identified through KEGG analysis to be involved in cancer formation may be closely associated with the development of HCC.

Acknowledgements

The present project was supported by the social development fund of Nantong (grant no. HS2014063) and The 12th Talent Summit of Top Six Industries in Jiangsu Province (grant no. WSW-080).

References