Investigation of Potential Molecular Biomarkers for Diagnosis and Prognosis of AFP-Negative HCC

Introduction

Hepatocellular carcinoma (HCC) is one of the most common malignant cancers and the second leading cause of cancer-associated deaths worldwide with the 5-year survival rate of <5%.^1,2 Most of the HCC cases are diagnosed at an advanced stage due to the lack of specific biomarkers and poor clinical symptoms, and over 80% of them have a poor prognosis.³ Alpha-fetoprotein (AFP) is the most important and commonly used diagnostic indicator for HCC.^4,5 A previous study has indicated that AFP can promote cell growth and escape from the host immune surveillance of liver cancer cells in vivo, which is an independent risk factor associated with the prognosis.⁶ Only 33–65% of the tumors with a diameter of <3 cm in HCC patients had an elevated serum AFP level,^1,7 and nearly 30% of the patients had normal AFP levels.^8,9 Although AFP has its limitations as a biomarker for screening for liver cancer, it is still most widely used biomarker. Thus, it is crucial to identify new biomarkers for the diagnostic and prognostic evaluation of AFP-negative HCC and understand the underlying mechanisms of its tumorigenesis.

Non-coding RNA (ncRNA) is a group of noncoding RNA transcripts, which have many pivotal functions in various physiological and pathological processes.^10,11 Long non-coding RNAs (lncRNAs) can regulate gene expression indirectly although they have no protein-coding function.^12,13 As increasing number of studies have reported that the dysregulation of lncRNAs plays a vital role in the tumorigenesis of different types of cancer and are related to the poor prognosis in HCC,^14–16 it is necessary to understand the role of lncRNAs in cancer.

The competing endogenous RNA (ceRNA) hypothesis proposed that lncRNAs can competitively bind to microRNA (miRNA) and indirectly regulate its expression levels by acting as miRNAs sponges.^17–19 The ceRNA network has been shown to play a critical role in HCC pathogenesis and progression.^15,20,21 However, the underlying molecular mechanisms of AFP-negative HCC and the role of ceRNA in its tumorigenesis are still unknown.

In this study, the lncRNA, miRNA, and mRNA gene expression data of HCC were downloaded from the TCGA database and screened for AFP-negative HCC. After identifying the differentially expressed RNAs (DERNAs) between AFP-negative HCC samples and adjacent non-cancerous samples, the lncRNA-miRNA-mRNA regulatory network and protein-protein interaction (PPI) network for AFP-negative HCC were constructed. Receiver operating characteristic (ROC) curves of the top 10 hub genes were used to assess the diagnostic efficiency. Meanwhile, RNAs which played a key role in the prognosis of AFP-negative HCC among the ceRNA were identified. This study will contribute to assist clinical diagnosis of liver cancer and evaluate its prognosis when the level of AFP is negative.

Materials and Methods

Data Source and Pre-Treatment

The original gene expression data and clinical information including the AFP value of HCC patients were obtained from the TCGA database (https://cancergenome.nih.gov/). The miRNA-seq data from 150 AFP-negative (AFP <20 ng/mL) HCC samples and 50 adjacent non-tumorous samples, as well as mRNA-seq data from 147 AFP-negative HCC samples and 50 adjacent non-tumorous samples were used. As data from TCGA are publicly available, no further ethical approval was required for the conduct of this study.

Differential Expression Analysis

The RNA and miRNA sequencing data were downloaded to distinguish differentially expressed lncRNAs, miRNAs, and mRNAs between the AFP-negative HCC and normal samples. After the gene expression matrix was merged, the ensemble IDs were converted to gene names. The data was calibrated and standardized, and the differential expression was identified by applying the ‘edgeR’ package of R software. The False Discovery Rate (FDR) was used to correct the P value for determining the statistical significance of multiple testing. The dysregulated RNAs were identified based on the cut-off criteria FDR <0.05 and |log2 (fold change) |>1.

CeRNA Network Construction

The lncRNA-miRNA interactions were predicted by using the miRcode (http://www.mircode.org/).²² In addition, miRDB (http://www.mirdb.org/), miRTarBase (http://mirtarbase.mbc.nctu.edu.tw), and TargetScan (http://www.targetscan.org) were applied to predict the DEmiRNA target genes,^23–25 and only miRNA-mRNA pairs found in all the three databases were utilized to construct the ceRNA networks. To further increase the reliability of the ceRNA network, the overlap between DEmiRNA target genes and the DEmRNA in AFP-negative HCC obtained by Venn diagram was used as DEmRNA. Finally, the ceRNA network based on the matched DElncRNA–DEmiRNA and DEmiRNA–DEmRNA pairs was constructed to clarify the gene interactions in AFP-negative HCC. Cytoscape software²⁶ was used to visualize lncRNA-miRNA-mRNA network.

Functional Enrichment Analysis

To explore the potential biological processes and pathways of overlapping differentially expressed genes (DEGs), Gene Ontology (GO) enrichment at the significant level (p-value < 0.01) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses at the significant level (p-value < 0.05) was performed using Database for Annotation, Visualization and Integrated Discovery (DAVID: available online: http://david.abcc.ncifcrf.gov/).

PPI Network and Hub Genes

To explore the interrelationship of DEmRNAs, STRING (https://string-db.org/),²⁷ a widely used interactive tool, was utilized to set up a PPI network, which was then visualized using the Cytoscape software. The top ten well connected DEmRNAs were identified as hub genes with CytoHubba. ROC curves were generated to further explore the diagnostic value of the 10 hub genes, and the area under the curve (AUC) was used to evaluate the diagnostic efficiency.

Survival Analysis

To explore the prognostic DERNA signature, differentially expressed lncRNAs, miRNAs, and mRNAs associated with survival were analyzed using the ‘Survival’ package of R software. Survival curves were estimated based on Kaplan–Meier survival curve analysis and P < 0.05 was considered as statistically significant.

Results

Differentially Expressed lncRNA, mRNA, and miRNA

A total of 131 lncRNAs (87 upregulated and 44 downregulated), 185 miRNA (141 upregulated and 44 downregulated), and 1309 mRNAs (786 upregulated and 523 downregulated) were observed to be differentially expressed in AFP-negative HCC on the criteria that |log2 (fold change)| >1, FDR <0.05. The visual volcano plots are shown in Figure 1.

Figure 1 Volcano plots of differentially expressed lncRNA (A), miRNA (B) and mRNA (C) in AFP-negative HCC. Red dots indicate genes significantly upregulated and green indicates genes significantly downregulated.

Construction and Analysis of ceRNA Regulatory Network

The DElncRNA-DEmiRNA and DEmiRNA-DEmRNA pairs were used to construct the ceRNA network to further understand how lncRNA mediates mRNA by binding miRNA in AFP-negative HCC. The result showed that 81 DElncRNA-DEmiRNA interactions were identified, with 12 DElncRNAs targeting 23 DEmiRNAs from the miRcode database (Table 1). The miRDB, miRTarBase, and TargetScan databases were combined to predict the candidate mRNA targets of the DEmiRNAs to improve the reliability of the results. A total of 1119 target genes that intersected with 1309 DEmRNA were obtained, resulting in 74 DEmRNA (Figure 2). Finally, the ceRNA network comprising 12 DElncRNAs, 23 DEmiRNAs, and 74 DEmRNAs was visualized by Cytoscape software (Figure 3).

Table 1 12 DElncRNAs Target 23 DEmiRNAs from the miRcode Database

Figure 2 Venn diagram analysis for the intersection set between differentially expressed mRNA in AFP-negative HCC and miRNA target genes.

Figure 3 The lncRNA-miRNA-mRNA ceRNA network in AFP-negative HCC. The diamond represents lncRNA, the square represents miRNA and the circle denotes mRNA. The red nodes indicate upregulated RNAs and the green nodes indicate downregulated RNAs.

Function Enrichment

To explore the function of the ceRNA network, the GO and KEGG enrichment analysis of DEmRNAs in the ceRNA network was performed using DAVID to explore the function of the ceRNA network. As shown in Figure 4, 11 GO biological process (BP) terms such as “negative regulation of apoptotic process,” “response to drug,” “positive regulation of transcription,” and so on were significantly (P < 0.01) enriched. Also, in the cellular component (CC), the significantly (P < 0.01) enriched GO terms were “nucleus,” “membrane,” “nucleoplasm,” and’ “nuclear chromatin,” and those in the molecular function (MF) were “protein kinase binding,” “chromatin binding,” “protein binding,” “ATP binding,” and “beta-catenin binding.” The KEGG results revealed that 74 DEGs were involved in 26 significantly (P < 0.05) enriched pathways (Figure 5). Among these pathways, “p53 signaling pathway”, “Hepatitis B” “PI3K-Akt” were related to hepatocarcinogenesis. And some other pathways such as “Bladder cancer,” “Colorectal cancer,” “Small cell lung cancer,” “Prostate cancer,” “MicroRNAs in cancer,” and “Pathways in cancer” were related to cancer-related pathways.

Figure 4 The GO enrichment analysis of DEmRNAs in the ceRNA network of AFP-negative HCC.

Figure 5 The enriched KEGG pathway of DEmRNAs in the ceRNA network of AFP-negative HCC.

PPI Network of DEmRNAs in the ceRNA Network

A PPI network was built to elaborate interaction relationship between the 74 DEGs using the STRING database (interaction score >0.4) and displayed using the Cytoscape software after discarding the unrelated nodes (Figure 6A). Moreover, the top ten DEmRNAs obtained by CytoHubba were identified as the hub genes including CCNB1, E2F1, RRM2, MCM4, EZH2, PBK, ATAD2, CCNE1, ESR1, CHAF1A (Figure 6B). To explore whether the top 10 hub genes have excellent diagnostic efficiency in AFP-negative HCC patients, the ROC curves based on the TCGA database were performed. The results (Figure 7) showed the AUC of CCNB1, E2F1, RRM2, MCM4, EZH2, PBK, ATAD2, CCNE1, ESR1, and CHAF1A to be 0.954, 0.949, 0.916, 0.848, 0.962, 0.929, 0.816, 0.873, 0.920, and 0.928, respectively.

Figure 6 (A) A PPI network of DEmRNAs involved in the ceRNA network of AFP-negative HCC. Red nodes represent the upregulated DEmRNAs and green nodes represent the downregulated DEmRNAs. (B) The top ten highly interacted hub genes were recognized by cytoHubba.

Figure 7 The ROC curves of the top ten hub genes. The AUC of CCNB1, E2F1, RRM2, MCM4, EZH2, PBK, ATAD2, CCNE1, ESR1, and CHAF1A were 0.954, 0.949, 0.916, 0.848, 0.962, 0.929, 0.816, 0.873, 0.920, and 0.928, respectively.

Survival Analysis of RNAs in the Network

To study the prognostic capability of genes in the ceRNA network, the overall survival of HCC patients was determined with Kaplan–Meier curve analysis between DERNAs (P < 0.05 as the screening criterion). The results (Figure 8) showed that 2 lncRNAs (LINC00261, LINC00482), 3 miRNAs (hsa-miR-93, hsa-miR-221, and hsa-miR-222), and 2 mRNAs (EGR2, LPCAT1) were significantly correlated with the overall survival. High expression of LINC00482, hsa-miR-221, hsa-miR-222, and LPCAT1 was associated with poor overall survival, suggesting their potential negative role in AFP-negative HCC patients. However, high expression of LINC00261, hsa-miR-93 EGR2 was associated with good overall survival, which may be designated as AFP-negative HCC suppressor genes.

Figure 8 Kaplan–Meier survival curves of lncRNAs (A and B), miRNA (C–E), mRNA (F and G) associated with overall survival in AFP-negative HCC.

Discussion

AFP, as a conventional biomarker widely used in clinical practice, has provided critical clues for the diagnosis of HCC.²⁸ However, there are still approximately a number of HCC patients are negative for AFP.⁹ In addition, the AFP level is an independent risk factor correlated with tumor differentiation, TNM stage, and survival of patients with HCC.²⁹ Therefore, it is vital to identify more potential biomarkers for the effective diagnostic and evaluation of AFP-negative HCC.

LncRNA and miRNA, as non-coding RNA, are closely associated with the pathogenesis and progression of tumors. The emergence of lncRNA-miRNA-mRNA ceRNA revealed the regulatory relationship between non-coding and coding RNA and provided a novel guiding theory to explain the mechanism of tumorigenesis.^17,30 Hence, to better understand the molecular pathogenesis involved in AFP-negative HCC, a ceRNA network was built based on the TCGA database. In this study, the DERNAs between AFP-negative HCC and normal samples were identified, and a ceRNA network comprising 12 DElncRNAs, 23 miRNAs, and 74 DEmRNAs was constructed. Among them, lncRNA FAM99B have been reported to suppress HCC progression, and the low expression of FAM99B was associated with advanced histologic grade, T-stage, and vascular invasion.³¹ Furthermore, lncRNA RUSC1-AS1 can promote HCC cell proliferation and inhibit HCC cell apoptosis through Notch signaling pathway, and high expression of RUSC1-AS1 indicated a poor prognosis.³² MiR-214 has been shown to be a suppressor gene in many types of cancers including HCC, and it can reduce EZH2 to suppress tumorigenesis in various types of cancer.^33,34 MiR-216a and miR-137 have been demonstrated to be associated with hepatocarcinogenesis and tumor recurrence in HCC.^35–37 The ceRNA network in this study identified many RNAs that play an important role in the occurrence and development of HCC.

The AFP level has been reported to be associated with the pathological grade, progression, and prognosis of HCC.^38,39 Several studies have identified many prognostic indicators for HCC. However, there are few studies on the prognostic indicators of AFP-negative HCC based on a ceRNA network. Therefore, in this study, the relationship between aberrantly expressed RNAs in the ceRNA and prognosis was assessed to determine the potential prognostic indicators for AFP-negative HCC. It was noted that 2 lncRNAs (LINC00261, LINC00482), 3 miRNAs (hsa-miR-93, hsa-miR-221, and hsa-miR-222) and 2 mRNAs (EGR2, LPCAT1) were correlated with the overall survival in AFP-negative HCC patients. Among these RNAs, LINC00261 has been reported to suppress cell proliferation and invasion in HCC, and associate with AFP levels and postoperative outcome of HCC.^40,41 Furthermore, the low expression of LINC00261 in AFP-negative HCC also showed a poor prognosis in the present study. Xuan L suggested LPCAT1 to be a potent target molecule for inhibiting HCC progression.⁴² The high expression of LPCAT1 demonstrated a poor prognosis in this study, indicative of a positive correlation with the pathologic stage of AFP-negative HCC.

To further understand the ceRNA network, GO and KEGG analyses were performed to reveal the functions and pathways of 74 DEmRNAs in the ceRNA. The GO analyses demonstrated that DEmRNAs may focus on positive regulation transcription, DNA-templated, negative regulation transcription from RNA polymerase promoter, protein binding, ATP binding, protein kinase binding, and chromatin binding. KEGG pathway analysis indicated that the DEmRNAs were significantly associated with PI3K-Akt, p53 signaling pathway, Hepatitis B, and many other tumor-related pathways such as miRNAs in cancer, pathways in cancer, colorectal cancer, small cell lung cancer, and prostate cancer.

After constructing a PPI network, 10 hub genes including CCNB1, E2F1, RRM2, and PBK were identified. Previous studies have reported that overexpression of CCNB1, PBK and RRM2 predicted poor survival in HCC, which may be potential therapeutic targets.^43–45 EZH2 can initiate molecular processes such as NF-κB activation, miRNA silencing, and tumor immune evasion, which are considered to be the basic characteristics of cancer,⁴⁶ and therefore, different types of EZH2 inhibitors have been developed for cancer therapy.⁴⁷ ATAD2, an emerging oncogene closely related to the tumorigenesis in multiple cancers, has been validated to be a promising drug target for treatment development.⁴⁸ As the hub genes are closely associated with tumorigenesis, we are interested in whether they can be used as potential good diagnostic indicators of AFP-negative HCC. The 10 hub genes identified in the present study demonstrated superior performance in the diagnosis of AFP-negative HCC, indicating their potential as diagnostic biomarkers for AFP-negative HCC.

In the present study, we identified the dysregulation RNAs and constructed a ceRNA network which may lead to the tumorigenesis in the AFP-negative HCC. Despite the results provide some candidate biomarkers for AFP-negative HCC, there are inevitably some limitations. More experimental validations are needed for clinical application. And whether there is a correlation between the ceRNA and HCC clinical progress needs more exploration.

Conclusion

In conclusion, a ceRNA was constructed and 2 lncRNAs (LINC00261, LINC00482), 3 miRNAs (hsa-miR-93, hsa-miR-221, hsa-miR-222), and 2 mRNAs (EGR2, LPCAT1) were associated with the prognosis of AFP-negative patients. And ten hub genes including EZH2, CCNB1, E2F1, PBK, CHAF1A, ESR1, RRM2, CCNE1, MCM4, and ATAD2 showed good diagnostic power in AFP-negative patients. Our research identified potential diagnostic and prognostic biomarkers for AFP-negative HCC, which provided novel insights into the tumorigenesis mechanism of AFP-negative HCC.