Screening therapeutic targets of ribavirin in hepatocellular carcinoma.

The objective of the present study was to screen the key genes of ribavirin in hepatocellular carcinoma (HCC) and provide novel therapeutic targets for HCC treatment. The mRNA expression datasets of GSE23031 and GSE74656, as well as the microRNA (miRNA) expression dataset of GSE22058 were downloaded from the Gene Expressed Omnibus database. In the GSE23031 dataset, there were three HCC cell lines treated with PBS and three HCC cell lines treated with ribavirin. In the GSE74656 dataset, five HCC tissues and five carcinoma adjacent tissues were selected. In the GSE22058 dataset, 96 HCC tissues and 96 carcinoma adjacent tissues were selected. The differentially expressed genes (DEGs) and differentially expressed miRNAs were identified via the limma package of R. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis was performed with the Database for Annotation, Visualization and Integrated Discovery. The target mRNAs of DEMs were obtained with TargetScan. A total of 559 DEGs (designated DEG-Ribavirin) were identified in HCC cells treated with ribavirin compared with PBS and 632 DEGs (designated DEG-Tumor) were identified in HCC tissues compared with carcinoma adjacent tissues. A total of 220 differentially expressed miRNAs were identified in HCC tissues compared with carcinoma adjacent tissues. In addition, 121 GO terms and three KEGG pathways of DEG-Ribavirin were obtained, and 383 GO terms and 25 KEGG pathways of DEG-Tumor were obtained. A total of five key miRNA-mRNA regulated pairs were identified, namely miR-183→CCNB1, miR-96→DEPDC1, miR-96→NTN4, miR-183→NTN4 and miR-145→NTN4. The present study indicated that certain miRNAs (including miR-96, miR-145 and miR-183) and mRNAs (including NAT2, FBXO5, CCNB1, DEPDC1 and NTN4) may be associated with the effects of ribavirin on HCC. Furthermore, they may provide novel therapeutic targets for HCC treatment.

Introduction

Hepatocellular carcinoma (HCC) is the fifth most common malignancy and arises most frequently in patients with cirrhosis (1). It is the second most common cause of cancer-associated mortality globally with 1.6 million mortalities per year, and it is hypothesized that the high global incidence rate and late presentation of HCC may be responsible for this (2,3). Additionally, the general prognosis was poor with an overall survival rate between 3 and 5% in 2006 (4). Symptoms of HCC include yellow skin, bloating from fluid in the abdomen, easy bruising from blood clotting abnormalities, loss of appetite, unintentional weight loss, nausea, vomiting and tiredness (5,6). The primary risk factors for HCC were hepatitis C, hepatitis B, alcoholism, aflatoxin and cirrhosis of the liver (7–10). Liver transplantation, tyrosine kinase inhibitors and surgical resection are currently the primary treatment options (11–13). The treatment of HCC has not been fundamentally improved, which may be seen in the increasing morbidity and mortality each year (14). Ribavirin is an anti-viral drug used to treat hepatitis C, respiratory syncytial virus and other viral infections. If infection is persistent, ribavirin is often used in combination with peginterferon α-2b or peginterferon α-2a (15,16). It has been reported that hepatitis C infection was globally associated with 25% of HCC cases in 2006 (15). Therefore, ribavirin, by itself or in conjunction with peginterferon α-2b or pegylated interferon, has been used to treat HCC in patients with viral infections (17–20). Exploration of the genetic changes in HCC cells is necessary for the study of the pathogenesis and progression of HCC, as well as to develop effective treatments. In the present study, a microarray analysis of mRNA and microRNA (miRNA) was performed in the treatment of ribavirin on HCC, in order to identify possible biomarkers and provide novel potential therapeutic targets for HCC.

Materials and methods

Microarray data and data prx10-processing

The mRNA expression datasets of GSE23031 and GSE74656, as well as the miRNA expression dataset of GSE22058, (21–23) were downloaded from the Gene Expressed Omnibus database (http://www.ncbi.nlm.nih.gov/geo/). They were analyzed using the platforms GPL570 [HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array (Thermo Fisher Scientific, Inc., Waltham, MA, USA), GPL16043 GeneChip® PrimeView™ Human Gene Expression Array (with External spikx10-in RNAs; Thermo Fisher Scientific, Inc.) and GPL10457 Rosetta human miRNA qPCR array (Rosetta Inpharmatics; Merck Sharp & Dohme, Hoddesdon, UK), respectively. The mRNA data (GSE23031) contained three HCC cell lines treated with PBS and three HCC cell lines treated with ribavirin. In the GSE74656, five HCC tissues and five carcinoma adjacent tissues were selected for the study. In the GSE22058, 96 HCC tissues and 96 carcinoma adjacent tissues were selected to study. Robust Multi-Array Average (RMA) was an algorithm used to create an expression matrix from Affymetrix data (24). The raw data were converted into a recognizable format by R, and the RMA was used for correction and normalization.

Differential expression analysis

The differentially expressed genes (DEGs) were identified via the limma package V3.32.10 (http://www.bioconductor.org/packages/3.5/bioc/html/limma.html) (25). According to the criteria: P<0.05 and |log(fold change)|>1, the DEGs were identified in HCC cells treated with ribavirin compared with PBS and designated DEG-Ribavirin. With the same criteria, the DEGs were identified in HCC tissues compared with their matched adjacent tissues and designated DEG-Tumor. Additionally, the differentially expressed miRNAs (DEMs) were obtained in HCC tissues compared with carcinoma adjacent tissues with P<0.05 and |log(fold change)|>0.3.

Functional and pathway enrichment analysis

The Database for Annotation, Visualization and Integrated Discovery (https://david.ncifcrf.gov/) (26) is a widely-used web-based tool for functional and pathway enrichment analysis. In the present study, it was used to perform Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis of DEG-Ribavirin and DEG-Tumor data. The GO terms and KEGG pathways were selected with P<0.05.

Comparison of DEGs and screening of miRNA-mRNA regulated pairs

The overlapped DEGs of DEG-Ribavirin and DEG-Tumor were selected, and the overlapped DEGs with opposite expression between DEG-Ribavirin and DEG-Tumor were also selected. The TargetScan database was used to predict biological target mRNAs of miRNAs that matched the seed region of each miRNA (27). The target mRNAs of DEMs were then selected using TargetScan. The key miRNAs, which regulated the overlapped DEGs with opposite expression between DEG-Ribavirin and DEG-Tumor, were identified. Subsequently, the miRNA-mRNA regulated pairs were constructed.

Results

DEGs. A total of 559 DEGs (269 upregulated and 290 downregulated) and 623 DEGs (272 upregulated and 351 downregulated) were identified in DEG-Ribavirin and DEG-Tumor. The heat map of them and the top 30 most significant DEGs are presented in Figs. 1 and 2, Tables I and II, respectively. A total of 220 DEMs were obtained. The 30 most significant DEMs are presented in Table III.

Figure 1.

Heatmap of differentially expressed genes in hepatocellular carcinoma cells treated with ribavirin compared with PBS.

Figure 2.

Heatmap of differentially expressed genes in hepatocellular carcinoma tissues compared with their matched adjacent tissues.

Table I.

Top 30 most significant differentially expressed genes in hepatocellular carcinoma cells treated with ribavirin compared with PBS.

Gene	logFC	Average expression	P-value
NCF2	3.636237	10.08812	4.01×10−13
TXNIP	3.04897	9.867782	5.95×10−12
PRORSD1P	2.527183	6.729178	4.91×10−11
NPPB	3.0072	12.09607	8.26×10−11
CYR61	2.018572	10.7489	9.19×10−11
PLA2G4C	1.93915	5.61348	9.29×10−11
CDKN2B	2.282063	7.523593	9.52×10−11
LEAP2	2.52559	9.937438	1.39×10−10
DUSP4	1.86831	8.419143	3.18×10−10
FLVCR1-AS1	2.01682	8.90536	4.17×10−10
ASH1L-AS1	1.803117	9.666342	4.20×10−10
CXCL3	1.797497	8.092698	5.80×10−10
GLIPR2	1.688522	9.089844	5.91×10−10
CYP1A1	1.560867	11.576	7.61×10−10
EID2B	1.667057	8.101158	9.18×10−10
JUNB	1.574027	8.65205	1.10×10−09
BTG2	1.498662	9.091506	1.41×10−09
PPL	1.94965	8.769882	2.08×10−09
ZNF436-AS1	2.570947	8.267133	2.12×10−09
TIGD7	1.654807	6.61165	2.17×10−09
LOC284513	1.59588	6.487373	2.40×10−09
TUFT1	1.50993	10.31677	2.43×10−09
CTSE	1.439267	11.22817	2.48×10−09
ERP27	1.376543	10.6632	2.51×10−09
UCA1	1.611053	10.73581	2.60×10−09
HSD3B1	1.699597	5.372037	2.61×10−09
GATA6-AS1	1.971443	9.336978	2.99×10−09
LOC100134822	1.45604	7.840603	3.38×10−09
THUMPD3-AS1	1.365701	7.84944	3.41×10−09
GDA	1.482742	8.376296	3.62×10−09

logFC, log fold-change.

Table II.

Top 30 most significant DEGs in HCC tissues compared with carcinoma adjacent tissues.

Gene	logFC	Average expression	P-value
CTHRC1	2.856945	5.825291	6.69×10−07
PEA15	1.468145	8.69977	1.04×10−06
CENPE	2.043658	6.056165	1.29×10−06
C21orf56	1.237448	5.670084	1.54×10−06
DBN1	1.383559	6.53977	1.76×10−06
GLA	1.380323	7.172235	2.00×10−06
DDX39	1.17872	7.721654	2.63×10−06
MPV17	1.190142	9.079669	2.67×10−06
RFX5	1.417855	7.558042	2.89×10−06
TMEM144	1.379456	5.655336	3.59×10−06
ASNS	2.515006	6.603396	3.63×10−06
SLC38A6	1.678498	7.591919	3.93×10−06
GRAMD1A	1.013492	7.228709	4.65×10−06
COMMD8	1.109678	8.353807	7.57×10−06
YWHAZ	1.008572	8.92205	7.94×10−06
PLXNC1	1.451381	5.772042	1.18×10−05
PLVAP	1.001818	5.796379	1.20×10−05
SHCBP1	1.327673	4.917339	1.39×10−05
LAMC1	1.104793	6.905647	1.46×10−05
ANXA2P2	1.89594	11.61567	1.52×10−05
ACTR3	1.094921	9.606995	1.63×10−05
CCDC88A	1.002649	6.195715	1.64×10−05
E2F3	1.093946	6.541921	1.75×10−05
FAM118B	1.006248	6.9915	1.76×10−05
ZNF354A	1.069078	5.918941	1.81×10−05
RASGEF1A	2.021872	4.908126	1.84×10−05
NUP37	1.362635	6.839194	2.01×10−05
ESM1	1.291572	4.762991	2.03×10−05
KPNA2	2.606398	9.341649	2.13×10−05
DCUN1D5	1.38889	8.463629	2.34×10−05

logFC, log fold-change; HCC, hepatocellular carcinoma; DEGs, differentially expressed genes.

Table III.

Top 30 most significant differentially expressed miRNA in HCC tissues compared with carcinoma adjacent tissues.

miRNA	logFC	P-value
hsa-mir-188	0.36842	2.41×10−39
hsa-mir-106b	0.32071	1.52×10−36
hsa-mir-214	−0.54861	7.34×10−36
hsa-mir-93	0.32346	1.73×10−35
hsa-mir-10a	−0.51702	6.33×10−35
hsa-mir-199a-1	−0.78765	6.28×10−33
hsa-mir-199a-2	−0.73035	5.57×10−31
hsa-mir-301	0.5629	1.16×10−28
hsa-mir-424	−0.31246	1.72×10−28
hsa-mir-33	0.29308	7.92×10−26
hsa-mir-324-5p	0.37868	1.01×10−23
hsa-mir-25	0.20282	1.53×10−22
hsa-mir-125b	−0.37741	1.95×10−22
hsa-mir-339	0.21825	3.36×10−21
hsa-mir-145	−0.35693	4.27×10−21
hsa-mir-148b	0.20955	6.76×10−21
hsa-mir-151	0.25728	9.53×10−21
hsa-mir-221	0.35829	1.39×10−20
hsa-mir-18a	0.41731	3.09×10−20
hsa-mir-130b	0.38724	3.75×10−20
hsa-mir-195	−0.34663	6.72×10−20
hsa-mir-99a	−0.38524	8.55×10−20
hsa-mir-15b	0.27273	1.80×10−19
hsa-mir-183	0.42758	2.31×10−19
hsa-mir-222	0.30098	2.62×10−18
hsa-mir-125a	−0.26924	4.47×10−18
hsa-mir-378	−0.33767	7.60×10−18
hsa-mir-101	−0.24801	1.85×10−17
hsa-mir-331	0.20462	2.33×10−17
hsa-mir-200b	−0.6996	2.38×10−17

HCC, hepatocellular carcinoma; logFC, log fold-change; mir, microRNA.

GO terms and KEGG pathways

A total of 121 GO terms and 3 KEGG pathways (cell cycle pathway, p53 signaling pathway and glycine, serine and threonine metabolism pathway) of DEG-Ribavirin were obtained. A total of 383 GO terms and 25 KEGG pathways of DEG-Tumor were obtained. The top 20 enriched GO terms of DEG-Ribavirin and DEG-Tumor are presented in Tables IV and V, respectively. The enriched KEGG pathways of DEG-Ribavirin and DEG-Tumor are presented in Tables VI and VII, respectively.

Table IV.

Top 20 enriched GO terms of DEG-Ribavirin.

Category	GO ID	Go name	Gene number	P-value
BP	GO:0022403	Cell cycle phase	29	1.66×10−06
BP	GO:0007049	Cell cycle	43	1.66×10−06
BP	GO:0000279	M phase	24	8.49×10−06
BP	GO:0007067	Mitosis	19	1.02×10−05
BP	GO:0000280	Nuclear division	19	1.02×10−05
BP	GO:0000087	M phase of mitotic cell cycle	19	1.29×10−05
BP	GO:0048285	Organelle fission	19	1.73×10−05
BP	GO:0000278	Mitotic cell cycle	25	1.88×10−05
BP	GO:0022402	Cell cycle process	32	3.23×10−05
BP	GO:0008283	Cell proliferation	27	3.62×10−05
BP	GO:0031497	Chromatin assembly	11	5.88×10−05
BP	GO:0065004	Protein-DNA complex assembly	11	8.65×10−05
BP	GO:0006334	Nucleosome assembly	10	2.36×10−04
BP	GO:0051726	Regulation of cell cycle	21	2.43×10−04
BP	GO:0051301	Cell division	19	4.39×10−04
BP	GO:0034728	Nucleosome organization	10	5.08×10−04
BP	GO:0006323	DNA packaging	11	6.79×10−04
CC	GO:0000786	Nucleosome	8	9.25×10−04
BP	GO:0010033	Response to organic substance	33	0.001109
CC	GO:0005819	Spindle	12	0.001114

GO, Gene Ontology; DEGs, differentially expressed genes; CC, cellular component; BP, biological process.

Table V.

Top 20 enriched GO terms of DEG-tumor.

Category	GO ID	Go name	Gene number	P-value
BP	GO:0000279	M phase	49	1.22×10−17
BP	GO:0022403	Cell cycle phase	54	7.41×10−17
BP	GO:0007067	Mitosis	39	1.19×10−16
BP	GO:0000280	Nuclear division	39	1.19×10−16
BP	GO:0000087	M phase of mitotic cell cycle	39	2.21×10−16
BP	GO:0000278	Mitotic cell cycle	50	4.09×10−16
BP	GO:0048285	Organelle fission	39	5.79×10−16
BP	GO:0022402	Cell cycle process	61	5.37×10−15
MF	GO:0048037	Cofactor binding	39	3.37×10−14
BP	GO:0016054	Organic acid catabolic process	26	4.73×10−14
BP	GO:0046395	Carboxylic acid catabolic process	26	4.73×10−14
CC	GO:0005819	Spindle	30	7.20×10−14
BP	GO:0007049	Cell cycle	71	9.42×10−14
BP	GO:0055114	Oxidation reduction	62	3.81×10−13
BP	GO:0007059	Chromosome segregation	21	2.69×10−12
MF	GO:0009055	Electron carrier activity	34	3.23×10−12
CC	GO:0000793	Condensed chromosome	26	5.97×10−12
CC	GO:0000777	Condensed chromosome kinetochore	18	1.41×10−11
MF	GO:0050662	Coenzyme binding	29	6.32×10−11
CC	GO:0000779	Condensed chromosome, centromeric region	18	1.37×10−10

GO, Gene Ontology; DEGs, differentially expressed genes; CC, cellular component; BP, biological process; MF, molecular foundation.

Table VI.

Enriched KEGG pathways of DEG-Ribavirin.

Category	Pathway name	Gene number	P-value
KEGG_PATHWAY	hsa04110: Cell cycle	15	1.19×10−05
KEGG_PATHWAY	hsa00260: Glycine, serine and threonine metabolism	5	0.011495
KEGG_PATHWAY	hsa04115: p53 signaling pathway	6	0.045851

DEGs, differentially expressed genes; KEGG, Kyoto Encyclopedia of Genes and Genomes.

Table VII.

Enriched KEGG pathways of DEG-tumor.

Category	Pathway name	Gene number	P-value
KEGG_PATHWAY	hsa00071: Fatty acid metabolism	15	1.37×10−09
KEGG_PATHWAY	hsa00280: Valine, leucine and isoleucine degradation	14	5.58×10−08
KEGG_PATHWAY	hsa04110: Cell cycle	21	1.08×10−06
KEGG_PATHWAY	hsa00830: Retinol metabolism	12	3.12×10−05
KEGG_PATHWAY	hsa00380: Tryptophan metabolism	10	7.32×10−05
KEGG_PATHWAY	hsa00650: Butanoate metabolism	9	1.33×10−04
KEGG_PATHWAY	hsa00250: Alanine, aspartate and glutamate metabolism	8	4.63×10−04
KEGG_PATHWAY	hsa04114: Oocyte meiosis	15	5.67×10−04
KEGG_PATHWAY	hsa00640: Propanoate metabolism	8	5.69×10−04
KEGG_PATHWAY	hsa03320: PPAR signaling pathway	11	0.001281
KEGG_PATHWAY	hsa00980: Metabolism of xenobiotics by cytochrome P450	10	0.00174
KEGG_PATHWAY	hsa00910: Nitrogen metabolism	6	0.003688
KEGG_PATHWAY	hsa00590: Arachidonic acid metabolism	9	0.004249
KEGG_PATHWAY	hsa00140: Steroid hormone biosynthesis	8	0.005165
KEGG_PATHWAY	hsa00982: Drug metabolism	9	0.007941
KEGG_PATHWAY	hsa00591: Linoleic acid metabolism	6	0.008896
KEGG_PATHWAY	hsa00340: Histidine metabolism	6	0.010346
KEGG_PATHWAY	hsa04115: p53 signaling pathway	9	0.013645
KEGG_PATHWAY	hsa00260: Glycine, serine and threonine metabolism	6	0.013718
KEGG_PATHWAY	hsa00410: β-Alanine metabolism	5	0.017838
KEGG_PATHWAY	hsa04920: Adipocytokine signaling pathway	8	0.036604
KEGG_PATHWAY	hsa00620: Pyruvate metabolism	6	0.037809
KEGG_PATHWAY	hsa00232: Caffeine metabolism	3	0.039507
KEGG_PATHWAY	hsa05222: Small cell lung cancer	9	0.042544
KEGG_PATHWAY	hsa04512: ECM-receptor interaction	9	0.042544

DEGs, differentially expressed genes; KEGG, Kyoto Encyclopedia of Genes and Genomes.

miRNA-mRNA regulated pairs

There were 50 overlapped DEGs, with 32 [including N-acetyltransferase (NAT2) and F-box only protein 5 (FBXO5)] exhibiting opposite expression between DEG-Ribavirin, and DEG-Tumor. A heat map of the 32 overlapped DEGs is presented in Fig. 3. Furthermore, three DEMs (miR-96, miR-145 and miR-183) were revealed to correspond to three DEGs (CCNB1, DEPDC1 and NTN4) which were included in the aforementioned 32 overlapped DEGs. Finally, five miRNA-mRNA regulated pairs were selected between the above three DEGs and the three DEMs, namely miR-183→CCNB1, miR-96→DEPDC1, miR-96→NTN4, miR-183→NTN4 and miR-145→NTN4.

Figure 3.

Heatmap of the overlapped DEGs with opposite expression between DEG-Ribavirin and DEG-Tumor. Adj, adjacent; DEG, differentially expressed gene.

Discussion

In the present study, the DEGs in HCC cells treated with ribavirin compared with PBS treated HCC tissue, HCC tissues and carcinoma adjacent tissues, were firstly identified, and 32 overlapped DEGs with opposite expression between DEG-Ribavirin and DEG-Tumor were selected. It was notable that NAT2 and FBXO5 were two mRNAs of them with opposite expression between DEG-Ribavirin and DEG-Tumor. NAT2 serves a function in the metabolic activation and detoxification of aromatic amines, which in turn serves a function in the metabolism of aromatic and heterocyclic amines, and hydrazines via N-acetylation and O-acetylation (28). As early as in 1996, Agúndez et al (29) reported that the slow acetylation was associated with an increased risk of HCC. Furthermore, it has been demonstrated that NAT2 activity is associated with smoking-associated HCC (30–32). A number of previous studies that have investigated the association between NAT2 genotypes and HCC risk have been published (32–36). FBXO5, also known as early mitotic inhibitor-1, is a key cell-cycle regulator that promotes S-phase and M-phase entry by inhibiting anaphasx10-promoting complex/cyclosome activity (37). Zhao et al (38) revealed that FBXO5 was overexpressed in HCC, which is in agreement with the results of the present study, and also reported that FBXO5 may control tumor cell proliferation in HCC. In the present study, it was identified that the expression of NAT2 was lower in HCC cells and HCC tissues. However, expression was increased following treatment with ribavirin. However, FBXO5 was overexpressed in HCC cells and HCC tissues, and decreased following treatment with ribavirin. Therefore, it is suspected that NAT2 and FBXO5 may be biomarkers of ribavirin in the treatment of HCC.

The cell cycle has been demonstrated to be associated with the progression and migration of HCC (39–41), and regulation of the cell cycle is considered an effective strategy for HCC treatment (42–45). The p53 signaling pathway has been heavily studied and is reported to serve a function in the occurrence and development of HCC (45–49). The association between the glycine, serine and threonine metabolism pathway and HCC has been less studied, and the glycine, serine and threonine metabolism pathway was also enriched in DEG-Tumor tissues. In this study, only three KEGG pathways of DEG-Ribavirin were obtained, namely cell cycle, p53 signaling pathway and glycine, serine and threonine metabolism. Cell cycle was the most significantly enriched function in this study, which was identified from the enriched GO terms of DEG-Ribavirin (e.g. cell cycle phase, cell cycle and M phase) and DEG-Tumor (e.g. M phase and cell cycle phase), as well as the enriched KEGG pathways of DEG-Tumor. The results of the present study suggest that these three KEGG pathways may be associated with the pathogenesis and treatment of HCC; however, more in-depth research is required.

In the present study, three DEMs (miR-96, miR-145 and miR-183) were identified to correspond to three overlapped DEGs (CCNB1, DEPDC1 and NTN4) with opposite expression in DEG-Ribavirin and DEG-Tumor, and 5 miRNA-mRNA regulated pairs were selected, namely miR-183→CCNB1, miR-96→DEPDC1, miR-96→NTN4, miR-183→NTN4 and miR-145→NTN4. It has been demonstrated that miR-96 downregulation may suppress the growth of HCC (50), and miR-96 may promote cell proliferation and invasion through targeting ephrinA5 in HCC (51). Chen et al (52) considered serum miR-96 as a promising biomarker for HCC with chronic hepatitis B virus infection. Previous studies have demonstrated that miR-145 may inhibit proliferation, migration and invasion, as well as promote apoptosis in HCC (53–56). MiR-183 may also regulate the growth, invasion and apoptosis of HCC (57–59). It has been identified that these miRNA and mRNA are possible biomarkers of ribavirin in HCC, and they may regulate HCC through the 5 miRNA-mRNA pairs.

In conclusion, a number of miRNAs (e.g. miR-96, miR-145 and miR-183) and mRNAs (e.g. NAT2, FBXO5, CCNB1, DEPDC1 and NTN4) may be associated with the effects of ribavirin on HCC. Furthermore, they may provide novel therapeutic targets for drugs of HCC.