Construction of CDKN2A-related competitive endogenous RNA network and identification of GAS5 as a prognostic indicator for hepatocellular carcinoma

Abstract

BACKGROUND

Competitive endogenous RNA (ceRNA) is an innovative way of gene expression modulation, which plays a crucial part in neoplasia. However, the intricacy and behavioral characteristics of the ceRNA network in hepatocellular carcinoma (HCC) remain dismal.

AIM

To establish a cyclin dependent kinase inhibitor 2A (CDKN2A)-related ceRNA network and recognize potential prognostic indicators for HCC.

METHODS

The mutation landscape of CDKN2A in HCC was first explored using the cBioPortal database. Differential expression analysis was implemented between CDKN2A^high and CDKN2A^low expression HCC samples. The targeted microRNAs were predicted by lncBasev3.0, and the targeted mRNAs were predicted by miRDB, and Targetscan database. The univariate and multivariate analysis were utilized to identify independent prognostic indicators.

RESULTS

CDKN2A was frequently mutated and deleted in HCC. The single-cell RNA-sequencing analysis revealed that CDKN2A participated in cell cycle pathways. The CDKN2A-related ceRNA network-growth arrest specific 5 (GAS5)/miR-25-3p/SRY-box transcription factor 11 (SOX11) was successfully established. GAS5 was recognized as an independent prognostic biomarker, whose overexpression was correlated with a poor prognosis in HCC patients. The association between GAS5 expression and methylation, immune infiltration was explored. Besides, traditional Chinese medicine effective components targeting GAS5 were obtained.

CONCLUSION

This CDKN2A-related ceRNA network provides innovative insights into the molecular mechanism of HCC formation and progression. Moreover, GAS5 might be a significant prognostic biomarker and therapeutic target in HCC.

Key Words: Hepatocellular carcinoma, Competitive endogenous RNA, Molecular mechanism, Prognosis, Biomarker

Core Tip: In this study, the mutation landscape, and molecular biological mechanisms of cyclin dependent kinase inhibitor 2A (CDKN2A) in hepatocellular carcinoma (HCC) was first explored, and a CDKN2A-related competitive endogenous RNA axis was constructed. We identified growth arrest specific 5 as an independent prognostic biomarker and established a prognostic nomogram model for HCC. Moreover, we analyzed its methylation level, immune infiltration, and targeted agents, which might be an independent prognostic biomarker and therapeutic target in HCC.

INTRODUCTION

Hepatocellular carcinoma (HCC) takes up 75%-85% of primary liver cancer and is also the third leading cause of cancer deaths globally[1,2]. The development of HCC involves many genetic alterations and interactions, which is the molecular basis for its recurrence, metastasis, and drug resistance[3-5]. Although there have been some advances in HCC treatment, the 5-year survival rate for HCC patients is less than 17%[6]. The tumor recurrence rate is still up to 70% within five years after surgery or radiofrequency ablation for early HCC[7]. Therefore, seeking valuable prognostic predictors and therapeutic targets is imperative to ameliorate defects in the prevention, diagnosis, and treatment of HCC.

Competitive endogenous RNA (ceRNA) has received much scholarly attention in recent years, representing a novel model of gene expression regulation. ceRNA regulatory networks are more elaborate and complex than microRNA (miRNA) regulatory networks, involving many RNA molecules, including mRNAs, long noncoding RNAs (lncRNAs), and pseudogenes encoding genes[8-10]. ceRNA molecules increase the expression of target genes by competitively linking the same miRNA with mRNA, thereby diluting the concentration of free miRNA in the cell and reducing the inhibition of mRNA by miRNA[11,12]. Besides, Salmena et al[8] proposed that lncRNAs mostly functioned as endogenous competitors in the cytoplasm. Regarding HCC, the DLEU2L/TAOK1 axis was recognized as a valuable prognostic factor engaging in HCC[13]. Li et al[14] proposed that TRHDE-AS1 influenced the prognosis of hepatitis B virus-infected HCC via miR-23b/PKIA axis. According to Hao et al[15], the transcript produced from the pseudogene UBE2MP1 acted as a molecular sponge for miR-145-5p, promoting the proliferation and resistance to apoptosis of HCC cells in vitro.

Cyclin dependent kinase inhibitor 2A (CDKN2A) is frequently mutated or deleted in various tumors, including HCC, and is known to be an essential tumor suppressor gene[16]. Inactivation of CDKN2A may result in uncontrolled cell growth and proliferation, leading to cancer development. Notably, CDKN2A was one of the cuproptosis (copper-induced cell death) genes discovered by Tsvetkov et al[17]. Prior research has revealed that deleted CDKN2A was a common activity in HCC and was connected to invasion, metastasis, recurrence, and prognosis[18-20]. For example, Desjonqueres et al[21] mentioned that aberrant methylation changes of CDKN2A at early stages of malignant transformation in HCC correlated with increased cell proliferation.

In our study, we established a CDKN2A-related ceRNA network-GAS5/miR-25-3p/SOX11 axis in HCC (Figure 1). First, we sought differentially expressed RNAs (DEGs) between CDKN2A^high and CDKN2A^low expression HCC patients and obtained CDKN2A-related lncRNA-miRNA-mRNA triple regulatory networks. Subsequently, a hub ceRNA network was selected through differential expression and survival analysis, and GAS5 was considered an independent prognostic factor. Finally, we carried out methylation analysis, immune infiltration analysis and targeted drug analysis to explore potential biological roles of GAS5 in HCC. In conclusion, our study constructed GAS5/miR-25-3p/SOX11 axis and identified GAS5 as a significant prognostic indicator and therapeutic target in HCC.

Figure 1 Flowchart of construction of the cyclin dependent kinase inhibitor 2A-related competitive endogenous RNA network. HCC: Hepatocellular carcinoma; TCGA: The Cancer Genome Atlas; CDKN2A: Cyclin dependent kinase inhibitor 2A; GAS5: Growth arrest specific 5; SOX11: SRY-box transcription factor 11; DElncRNAs: Differentially expressed long non-coding RNAs; DEmiRNAs: Differentially expressed microRNAs; DEmRNAs: Differentially expressed messenger RNAs.

MATERIALS AND METHODS

Data gathering and handling

The RNA-Seq expression profiles of 374 HCC samples and 50 normal samples were provided by the The Cancer Genome Atlas (TCGA) database (http://portal.gdc.com), together with corresponding clinical information of HCC patients. The raw read counts of gene expression were converted into transcripts per million. The clinical information of patients in the TCGA cohort is illustrated in Supplementary Table 1.

Besides, the mutation landscape of CDKN2A in HCC was investigated via the cBioPortal for Cancer Genomics (http://www.cbioportal.org/).

Single-cell RNA-seq analysis

The Single-cell RNA-seq (scRNA-seq) dataset (GSE112271) of HCC was acquired from the GEO database (https://ncbi.nlm.nih.gov/geo/). The “Seurat” package was used to preprocess before dimensionality reduction. Then the principal component analysis (PCA) was performed to identify the top 20 PC based on the top 2000 highly variable genes. Cluster visualization was performed by uniform manifold approximation and projection reduction, and the marker genes for each cluster were selected using the FindAllMarkers function with the adjusted P value < 0.01 and absolute log2 (fold change) value > 1. Finally, each cluster was annotated based on the corresponding canonical marker gene.

For metabolism analysis, the mean expression levels of the cells contained in various clusters were initially determined by the Average Expression function. For every cluster, the scores of the relevant pathways were determined using the “GSVA” R package[22], which were finally visualized by the “pheatmap” R package.

DEGs screening

For DEGs in CDKN2A^high and CDKN2A^low HCC samples, we selected differentially expressed long non-coding RNAs (DElncRNAs), differentially expressed microRNAs (DEmiRNAs), and differentially expressed messenger RNAs (DEmRNAs) with limits of all P < 0.05 and |logFC| > 0.5, |logFC| > 0.3, |logFC| > 0.5, respectively. For DEGs between HCC samples and adjacent non-tumor samples, we set P < 0.05 and |logFC| > 0.5 as the cut-off value of DElncRNAs, P < 0.05 and |logFC|> 0.3 as the cut-off value of DEmiRNAs, P < 0.05 and |logFC| > 0.7 as the cut-off value of DEmRNAs.

Construction and validation of the ceRNA network in HCC

The following traits of ceRNA are based on its notion and mechanism: (1) The presence of same miRNA binding sites for ceRNAs; (2) ceRNA are modulated by miRNAs; (3) the mutual regulatory relationship between ceRNAs and the consistent trend in expression levels; and (4) lncRNAs mainly function as ceRNAs in the cytoplasm. Hence, the ceRNA network was constructed based on the following procedures: (1) LncBasev.3 (https://diana.ece.uth.gr/Lncbase-v3/interactions) was utilized to forecast potential miRNAs targeted by DElncRNAs and lncRNA-miRNA matched bindings; (2) miRDB (http://www.mirdb.org/) and Targetscan (version 7.2, http://www.targetscan.org/) were applied to predict target genes of the DEmiRNAs and search for miRNA-mRNA matched bindings; (3) “VennDiagram” package was used to screen the target genes that intersected with DEmRNAs in this study; (4) LncATLAS (https://Lncatlas.crg.eu/) database[23] was utilized to select cytoplasmic lncRNAs for the subsequent analysis; (5) Cytoscape software (version 3.9.1) was used to construct and draw the lnRNA-miRNA-mRNA triple regulatory network; and (6) The Cytoscape plug-in CytoHubba was applied to determine the hub triple regulatory network which ranked top 15 by maximum clique centrality (MCC) algorithm.

Functional enrichment research and gene network analysis

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) research of DEmRNAs was performed using “clusterProfiler” package[24]. To investigate the potential interplay of 15 hub genes, we applied an interactive network analysis through the GENEMANIA website[25].

Clinical correlation analysis

The univariate and multivariate Cox regression was utilized to recognize independent prognostic indicators for HCC. Using “rms” package, the nomogram was established to incorporate every predictive variable in the multivariate Cox regression model. The model’s predictive ability was assessed using the calibration curve.

Validation of differential expression of GAS5

For validation, four independent datasets (GSE45267, GSE45436, GSE55092, and GSE64041) were collected from the GEO database. The data of GSE45267, GSE45436, GSE55092 were generated using the platform GPL570, while the data of GSE64041 was produced by the platform GPL6244. We adopted a column bar graph to present differential expression of GAS5 in a variety of datasets using GraphPad Prism v9.3.1.

To enhance the statistical credibility, we implemented a meta-analysis of differential expression results of diverse datasets by inputting mean and standard deviation. We took Q test (I2 statistics) result to measure the heterogeneity between different data. In addition, a forest plot was created to display the mean difference and correlated 95% confidence interval of GAS5. All procedures were accomplished by the software Review Manager 5.4.

Methylation and immune infiltration analysis of GAS5

DNA methylation, a common epigenetic mechanism, could modulate gene expression. We first used the GSCALite database[26] to investigate the correlation between methylation and GAS5 expression in LIHC. Then we separately utilized UALCAN and DiseaseMeth version 2.0 to evaluate the methylation degree of GAS5. Besides, the methylation site of GAS5 DNA sequence correlation with gene expression was visualized through MEXPRESS.

The TIMER database[27] was utilized to evaluate the abundance degree of different immune cell types in HCC tumor microenvironment. Additionally, we explored the association between three significant immune checkpoints (ICKs), including PDCD1, CD274, HAVCR2, and GAS5 in HCC using the Spearman correlation coefficient, since the expression of ICK-related genes influences the immunotherapy effect of ICK blockades[28].

Screening traditional Chinese medicine effective agents targeting GAS5

The CTD database (http://ctdbase.org/) was utilized to recognize traditional Chinese medicine (TCM) effective components targeting GAS5, and PubChem database (pubchem.ncbi.nlm.nih.gov) was applied to visualize their structures.

Statistical analysis

As for statistical description, numerical variables were depicted as mean ± SD, while categorical variables were depicted as frequency and proportion. The student’s t test and nonparametric tests were utilized to make a comparison between two groups of data. All analysis methods and R packages were conducted by R (3.6.3) software. A difference of P < 0.05 was defined as statistically significant.

RESULTS

The mutation landscape of CDKN2A in HCC

We first explored genomic alterations, copy number, and mRNA expression of CDKN2A using the cBioPortal database. As shown in Figure 2A, CDKN2A had an 8% genomic alteration rate, among which deep depletion was the most common type. Furthermore, HCC samples harboring CDKN2A deletion exhibited the highest mRNA expression, while HCC samples containing CDKN2A deep deletion showed the lowest mRNA expression (Figure 2B). There was a positive association between mRNA expression and CDKN2A copy number values in HCC samples in terms of both Spearman correlation coefficient and Pearson correlation coefficient (Figure 2C). As shown in Figure 2D, there were nine mutations occurring at different positions in the CDKN2A protein for HCC patients, including eight putative drivers (five missense mutations, two truncating mutations, and one splice mutation) and one variant of uncertain significance (one missense mutation). In addition, we found that CDKN2A participated in cell cycle pathways in HCC using pathway sources from the cBioPortal database (Figure 2E).

Figure 2 The mutation landscape of cyclin dependent kinase inhibitor 2A in hepatocellular carcinoma based on cBioPortal database. A: OncoPrint plot indicated different types and proportions of cyclin dependent kinase inhibitor 2A (CDKN2A) mutations; B and C: The association between CDKN2A copy number and mRNA expression are shown in the dot plot (B) and correlation plot (C); D: Demonstration of mutation sites, mutation types, and mutation cases on the two-dimensional structure of the CDKN2A protein; E: The cell cycle pathway mapper of CDKN2A in hepatocellular carcinoma. CDKN2A: Cyclin dependent kinase inhibitor 2A.

Given that genetic mutation can drive tumorigenesis and cause aberrant cell proliferation, CDKN2A was very likely to participate in HCC progression.

Single-cell analysis re-confirmed that CDKN2A was involved in cell cycle pathways in HCC

One scRNA-seq dataset (GSE112271) was collected to further explore CDKN2A expression in the cells from HCC tissues. As shown in Figure 3A, three quality measures were applied. After the top 2000 highly variable genes were screened for clustering (Figure 3B), we utilized PCA to identify 20 PCs with significant difference (Figure 3C). Thirteen clusters were identified (Figure 3D), and three immune cell types (NK cell, B cell, and T cell) were annotated (Figure 3E). Notably, almost only T cells expressed CDKN2A (Figure 3F), More interestingly, we found that the many pathways enriched in T cells were associated with cell cycle, including “tumor proliferation signature”, “DNA repair”, “DNA replication”, and “G2M checkpoint” (Figure 3G), which re-indicated that CDKN2A may participate in cell cycle modulation.

Figure 3 Single cell analysis of cyclin dependent kinase inhibitor 2A. A: Quality control assessment and data filtering of single cell sequence data; B: Highly variable genes used for clustering and cell identification; C: Top 20 principal components were identified based on P value < 0.05; D: Visualization of different cell clusters; E: Corresponding annotation of the cell clusters; F: cyclin dependent kinase inhibitor 2A expression between different immune cell clusters; G: The pathways highly enriched in different immune cell clusters.

DEGs screening

To construct a ceRNA network related to CDKN2A, we first screened out the DElncRNAs, DEmiRNAs, and DEmRNAs in HCC samples with CDKN2A^high and CDKN2A^low expression groups as well as in HCC and adjacent non-tumor tissues using TCGA database. The P value and each |logFC| threshold was listed in the material and method part.

Totally, 2032 DElncRNAs, 23 DEmiRNAs, and 2055 DEmRNAs were selected from HCC samples with CDKN2A^high and CDKN2A^low expression groups. Meanwhile, 4544 DElncRNAs, 69 DEmiRNAs, and 6560 DEmRNAs were determined between HCC samples and normal liver tissue samples.

We employed volcano plots to visualize the distribution of DElncRNAs, DEmiRNAs, and DEmRNAs (Figure 4A-C; Supplementary Figure 1A-C), and the heatmaps to exhibit the expression of 15 significant variable genes in CDKN2A^high expression and CDKN2A^low expression HCC samples, as well as in HCC and normal samples (Figure 4D-F;Supplementary Figure 1D-F).

Figure 4 Screening of differentially expressed long non-coding RNAs, differentially expressed microRNAs, and differentially expressed messenger RNAs. A-C: The volcano plots describe 2032 differentially expressed long non-coding RNAs (A; |log2fold change| > 0.5 and adjusted P value < 0.05), 23 differentially expressed microRNAs (B; |log2fold change| > 0.3 and adjusted P value < 0.05), and 2055 differentially expressed messenger RNAs (C; |log2fold change| > 0.5 and adjusted P value < 0.05); D-F: The horizontal axis of the heatmap indicates the samples, and the vertical axis of the heatmap indicates 15 significant differentially expressed genes. ^aP < 0.05, ^bP < 0.01, ^cP < 0.001. lncRNA: Long non-coding RNAs.

Establishment of lncRNA-miRNA-mRNA triple regulatory networks

We then put 4544 DElncRNAs into lncBasev3 to find targeted miRNAs that overlapped with DEmiRNAs between CDKN2A^high HCC and CDKN2A^low HCC group and 3 targeted miRNAs (miR-25-3p, miR-210-3p, and miR-148a-3p) and their corresponding lncRNAs (MALAT1, NEAT1, SNHG7, and GAS5) were sorted out. We analyzed the localization of 4 LncRNAs through the lncATLAS database and found that only GAS5 was mainly located in the cytoplasm (Figure 5A). So, we put GAS5 targeting miRNA-miR-25-3p into the miRDB database and ENCORI database, respectively and acquired 63 shared targeted mRNAs. Then we utilized cytoscape software to construct the lncRNA-miRNA-mRNA triple regulatory network in HCC (Figure 5B) and the MCC algorithm to select top 15 score genes as hub genes and visualized them by Cytoscape plug-in CytoHubba (Figure 5C). Additionally, we implemented GO and KEGG analysis of the DEmRNAs in the network. We discovered that they were enriched in the “inclusion body”, “apical dendrite”, “actin-based cell projection”, and “positive regulation of neuron differentiation” (Figure 5D), indicating that they possibly played a crucial role in the growth and differentiation of cells. Finally, we built a gene network to analyze the interaction of hub genes by GENEMANIA website and revealed that these genes took part in positive regulation of cell development (Figure 5E).

Figure 5 Construction and functional enrichment analysis of long noncoding RNA-miRNA-mRNA triple regulatory networks. The diamonds denote long noncoding RNAs (lncRNAs), ellipses denote miRNAs, rectangles denote mRNAs. A: Cytoplasmic/nuclear localization of growth arrest-specific 5; B: The lncRNA-miRNA-mRNA triple regulatory network in hepatocellular carcinoma. Red indicates upregulated, and blue represents downregulated; C: Top 15 score genes identified by maximum clique centrality algorithm; D: Functional enrichment analysis of the differentially expressed messenger RNAs (DEmRNAs) in the network; E: Interactive network of the DEmRNAs. RCI: Relative concentration index; GAS5: Growth arrest-specific 5.

Construction and validation of the ceRNA regulatory network

To establish a ceRNA regulatory network that truly influenced the progression of HCC, we analyzed whether the 15 hub genes were differentially expressed and the relationship between 15 hub genes expression and OS in HCC patients. The results showed that 3 mRNAs (CACNA1, SOX11, and PDZD2) were differentially expressed, and only SOX11 was significantly linked to OS in HCC patients. Therefore, a CDKN2A-related ceRNA network-GAS5/miR-25-3p/SOX11 axis was constructed (Figure 6A). To validate whether the regulatory relationship is truly established, we explored the expression level of GAS5 and SOX11 in HCC patients and discovered that both GAS5 and SOX11 were remarkably upregulated in HCC patients (Figure 6B and C). Then, we examined whether SOX11 expression influenced the prognosis of HCC patients, and Kaplan-Meier analysis revealed that SOX11 overexpression was associated with poor OS in HCC patients (Figure 6D). Moreover, spearman correlation analysis revealed that CDKN2A expression was positively related to SOX11 expression (Figure 6E).

Figure 6 Construction and validation of the competitive endogenous RNA regulatory network. A: The model of the competitive endogenous RNA axis in the progression of hepatocellular carcinoma (HCC); B: The expression of growth arrest-specific 5 was significantly upregulated in HCC; C: the expression of SRY-box transcription factor 11 (SOX11) was significantly upregulated in HCC; D: Kaplan-Meier curves of SOX11 in patients with HCC in the low and high expression groups; E: Correlation analysis between SOX11 and cyclin dependent kinase inhibitor 2A in HCC. HCC: Hepatocellular carcinoma; CDKN2A: Cyclin dependent kinase inhibitor 2A; GAS5: Growth arrest-specific 5; SOX11: SRY-box transcription factor 11.

Clinical correlation analysis of the GAS5/miR-25-3p/SOX11 axis in HCC patients

GAS5 expression, SOX11 expression, and clinical characteristics were subjected to univariate and multivariate Cox analyses to identify independent prognostic factors. The univariate Cox regression analysis (Figure 7A) revealed that tumor diameter, distant metastasis, GAS5 expression, and SOX11 expression were tightly connected with a poor OS in HCC patients. The multivariate Cox regression analysis (Figure 7B) suggested that the only factors still associated with OS were tumor diameter and GAS5 expression. All variables in the multivariate Cox regression model were incorporated to establish a nomogram (Figure 7C). The calibration curve (Figure 7D) demonstrated the nomogram’s superior prediction ability. In addition, the hazard proportional curve (Figure 7E) showed Cox regression coefficients for SNHG3 at various prognostic times.

Figure 7 Clinical significance of growth arrest-specific 5 in hepatocellular carcinoma patients. A and B: Univariate Cox analysis (A) and multivariate Cox analysis (B) of clinical characteristics, growth arrest-specific 5 (GAS5) expression and SRY-box transcription factor 11 (SOX11) expression; C: The predictive nomogram model; D: The corresponding calibration curve; E: The hazard proportional curve of GAS5. 95%CI: 95% confidence interval; HR: Hazard ratio; GAS5: Growth arrest-specific 5; SOX11: SRY-box transcription factor 11.

Therefore, GAS5 might develop into an independent prognostic factor for HCC patients.

External validation of GAS5 overexpression

To verify that GAS5 was overexpressed in HCC tumor, we collected four independent GEO datasets (GSE64041, GSE45267, GSE55092, and GSE45436) and carried out a meta-analysis to aggregate research results. In GEO datasets, GAS5 exhibited a significant overexpression trend in HCC tumor tissues (GSE64041: t = 3.443, P = 0.0008; GSE45267: t = 6.596, P < 0.0001; GSE55092: t = 7.257, P < 0.0001; GSE45436: t = 4.220, P < 0.0001; Figure 8A-D). Because several datasets (including the TCGA-LIHC cohort) differ from one another, we used a random effects model in our meta-analysis. After combining all outcomes, we concluded that GAS5 was significantly overexpressed in HCC tissues (Figure 8E), which revealed the role of GAS5 as an oncogene.

Figure 8 External validation of growth arrest-specific 5 overexpression. A-E: Differential expression analysis of growth arrest-specific 5 in GSE64041 (A), GSE45267 (B), GSE55092 (C), GSE45436 (D), and meta-analysis outcome (E). ^aP < 0.05, ^bP < 0.01, ^cP < 0.001, ^dP < 0.0001. 95%CI: 95% confidence interval.

Correlation between methylation and GAS5 expression in HCC

To better understand the methylation level of GAS5 in HCC, we assessed the methylation of all hub genes through the GSCALite database. It is revealed that GAS5 expression was most correlated with methylation level among all hub genes (Figure 9A). We also discovered that the methylation level of GAS5 in normal samples was higher than that in HCC tumor samples both in the UALCAN database (Figure 9B) and the DiseaseMeth database (Figure 9C), which proved differential expression of GAS5 in HCC and normal samples from another point of view since methylation could inhibit gene expression. Furthermore, we detected ten methylation sites (cg08947084, cg06644515, cg03044573, cg16290996, cg17025683, cg07177756, cg05396044, cg03584783, cg06704773, and cg17868533) in the DNA sequences of GAS5 had a negative correlation with their expression at a standard of P value < 0.001 (Figure 9D).

Figure 9 Correlation between methylation and growth arrest-specific 5 expression in hepatocellular carcinoma. A: Spearman correlation analysis of methylation and expression of growth arrest-specific 5 (GAS5); B: Methylation was evaluated using UALCAN; C: Methylation was assessed using DiseaseMeth version 2.0; D: The methylation site of GAS5 DNA sequence association with gene expression was visualized using MEXPRESS. The expression of GAS5 is illustrated by the blue line in the center of the plot. Pearson’s correlation coefficients and P values for methylation sites and query gene expression are shown on the right side. GAS5: Growth arrest-specific 5.

Association between immune infiltration and GAS5 expression in HCC

As shown in Figure 10A, infiltration abundance of CD8⁺ T cells, macrophages, neutrophils, and dendritic cells were possibly associated with GAS5 gene copy number. In addition, we detected that GAS5 expression was remarkably positively associated with tumor purity and infiltration abundance of B cells, CD4⁺ T cells, CD8⁺ T cells, macrophages, neutrophils, and dendritic cells (Figure 10B).

Figure 10  Correlation analysis of growth arrest-specific 5 expression and immune infiltration in hepatocellular carcinoma. A: Association between growth arrest-specific 5 (GAS5) gene copy number and immune cell infiltration levels; B: Correlation of GAS5 expression with immune infiltration levels; C-E: Association between GAS5 expression and the expression of PDCD1, CD274, and HAVCR, respectively. ^aP < 0.05, ^bP < 0.01, ^cP < 0.001.

Besides, we explored the relationship between GAS5 expression and the expression of PDCD1, CD274, and HAVCR2 in HCC patients, respectively. The GAS5 expression showed a positive correlation with PDCD1 expression and a negative correlation with CD274 expression, while GAS5 expression bore no notable association with HAVCR2 expression (Figure 10C).

Identification of TCM effective agents targeting GAS5

To find out potential TCM active ingredients targeting GAS5, we acquired chemical components of TCM that may correlate with GAS5 using the CTD database. The top 10 chemicals were listed as follow: Valproic acid, bisphenol A, carbon tetrachloride, 1-methyl-4-phenylpyridinium, acetaminophen, tetrachlorodibenzodioxin, trichostatin A, 4-(5-benzo (1,3) dioxol-5-yl-4-pyridin-2-yl-1H-imidazol-2-yl) benzamide, aflatoxin B, and dorsomorphin. The structure of Valproic Acid, bisphenol A, 1-Methyl-4-phenylpyridinium, and Acetaminophen was visualized by the PubChem database (Supplementary Figure 2).

DISCUSSION

HCC pathophysiology is a multistep, intricate process, which entails gene mutation and signaling pathway abnormalities. Although substantial advancements have been made in diagnosing and managing HCC, additional studies are still required to understand the disease’s potential molecular mechanism fully and to enhance the prognosis for HCC patients. Herein, considering that CDKN2A was frequently mutated or deleted in HCC, we established a CKND2A-related ceRNA network-GAS5/miR-25-3p/SOX11, which might help to illuminate the molecular process of HCC carcinogenesis and identified GAS5 as an independent prognostic biomarker for HCC.

We first discovered that CDKN2A exhibited a relatively high mutation frequency in HCC patients. Cancer is featured in aberrant and limitless cell proliferation caused by genetic mutations, which are called “drivers” for they drive tumorigenesis, and their mutation form affects the homeostasis of a collection of crucial cell functions[29]. The single-cell analysis re-indicated that CDKN2A participated in cell cycle pathways in HCC. Therefore, CDKN2A was most likely to engage in the progression of HCC and was chosen for further study. Interestingly, CDKN2A, also called multiple tumor suppressor 1, has been proven to be the most common somatic mutated gene in metastatic tumors. If the CDKN2A gene is mutated or missing, the inhibition of the cyclinD-CDK4 complex is lifted, allowing it to function as an inhibitor of retinoblastoma protein and removing barriers in the cell cycle, which provides the cell unlimited proliferative capacity.

The significance of ceRNA network in the development and spread of malignancies has been increasingly clear in recent years. For example, several studies have found that the expression of key oncogene PTEN was downregulated by ceRNAs (e.g., PTEN pseudogene, VAPA, BCL11B, etc.), thus leading to the development of lung adenocarcinoma[30], renal clear cell carcinoma[31], and glioma[32]. Nevertheless, the predictive significance of ceRNA for HCC remains incompletely understood, and scant research has pinpointed ceRNA regulatory networks associated with CDKN2A that influence patient survival in HCC patients. Accordingly, we intended to construct a CDKN2A-related ceRNA network that likely predicts the prognosis of HCC and provides a potential target for the therapy of HCC.

First, after excluding lncRNAs whose subcellular localization was the nucleus, the lncRNA-miRNA-mRNA triple regulatory network containing one lncRNA (GAS5), one miRNA (miR-25-3p), and 63 mRNAs was acquired by bioinformatic analysis. Next, we identified the hub ceRNA network using the cytoHubba plug-in based on the MCC score. Functional analysis and gene network revealed that 15 hub genes might participate in the positive regulation of cell development. After differential expression analysis and survival analysis of 15 hub genes in HCC, we discovered that only differential expression of SOX11 was significantly linked to the prognosis of HCC patients. Furthermore, we demonstrated consistent expression trends of GAS5 and SOX11 in HCC patients, which lived up to the features of ceRNAs. Thus, CDKN2A-related ceRNA regulatory network- GAS5/miR-25-3p/SOX11 axis was established.

Based on published data, we observed that studies into the roles of GAS5, miR-25a-3p, and SOX11 in cancer or their connections to cancer have taken place. Hsieh et al[33] evaluated the influences of lncRNA GAS5 genetic variants on the clinicopathological characteristics of lung adenocarcinoma. Similarly, miR-25-3p was also identified as a stable tumor biomarker for its knockdown suppressed breast cancer cell proliferation and invasion[34]. As for SOX11, many studies have demonstrated that SOX11 participated in various cancer progression, such as lung adenocarcinoma[35], breast cancer[36], bladder cancer[37], neuroblastoma[38]. However, few studies have explored the role of ceRNA regulatory networks containing GAS5, miR-25-3p, or SOX11 in HCC.

In this study, we performed a clinical correlation analysis of the GAS5/miR-25-3p/SOX11 axis in HCC patients. GAS5 and SOX11 were markedly correlated with the OS of HCC patients in univariate Cox regression, and only GAS5 was remarkably linked to the OS of HCC patients in multivariate Cox regression. Hence, GAS5 was more likely to be an independent prognostic indicator for HCC patients. Besides, we established a nomogram containing all variables in the multivariate Cox regression model, exhibiting an excellent predictive capacity. We also demonstrated GAS5 abnormally high expression in four independent GEO datasets to improve the credibility of the conclusion. To further investigate the potential clinical significance of GAS5 for HCC treatment, we carried out methylation analysis, immune infiltration analysis and identified TCM effective drugs targeting GAS5. We validated the abnormal methylation level of GAS5 in HCC and detected ten methylation sites in the DNA sequences of GAS5. Both copy number and expression of GAS5 were tightly associated with the abundance of some immune cells, such as CD8+ cells, macrophages, neutrophils, and dendritic cells. Moreover, we discovered that GAS5 expression was positively linked to PDCD1 expression and negatively correlated with CD274 expression. These findings indicated that these differences induced by the GAS5/miR-25-3p/SOX11 axis might influence the tumor immune microenvironment, thus promoting the progression of HCC. Since Chinese herbal medicine was considered an essential source for discovering innovational medications[39,40], we identified ten potential TCM chemical components, which promisingly become potent methods for the treatment of HCC patients.

GAS5 has been widely considered as a potential cancer suppressor gene, with its downregulation associated with cancer in multiple different tissues[41-43]. In our study, however, we discovered that GAS5 was overexpressed in HCC tumor tissues, and its overexpression led to a poor prognosis of HCC patients. There was also some contradictory information about GAS5’s correlation with HCC. Two previous studies pointed out that downregulation of GAS5 promoted tumor progression in HCC[44,45]. In contrast, a study using bioinformatic tools revealed that GAS5 overexpression was significantly linked to shorter relapse free survival in the TCGA-LIHC cohort[46]. Similarly, Li et al[47] suggested that gain of GAS5 was an independent indicator of poor prognosis in HCC from the TCGA-LIHC dataset. Hence, we predicted that GAS5 might have its distinctive function and involve in the occurrence and progression of LIHC. Besides, GAS5 might exhibit different roles based on different genetic backgrounds.

Of course, these findings need to be demonstrated by detailed experiments in vivo and in vitro in the future, and there are still some defects in our study. First, although we validated the robustness of the predictive effect of GAS5 in the TCGA cohort, more prospective studies with a suitably high sample size are needed to confirm each other as this study was solely retrospective. Second, we are supposed to detect the function and mechanism of the GAS5/miR-25-3p/SOX11 axis in HCC by in vivo and in vitro experiments.

CONCLUSION

In conclusion, we constructed a CDKN2A-related ceRNA network-GAS5/miR-25-3p/SOX11 for HCC patients, which could possibly serve as a prognostic indicator for clinical outcome and therapeutic responsiveness, provide innovative insights into the molecular basis of HCC initiation and development. Besides, regulating GAS5 expression might be an effective therapeutic strategy in inducing immune infiltration and immunotherapy for HCC.

ARTICLE HIGHLIGHTS

Research background

Hepatocellular carcinoma (HCC) is a tumor with multi-etiology and multi-pathway involvement, and it is characterized by a low 5- year survival rate. Competitive endogenous RNA (ceRNA) is an innovative way of gene expression modulation, which plays a crucial part in HCC. Nevertheless, the complexity and behavioral features of the ceRNA network in HCC are not fully elucidated.

Research motivation

To construct a cyclin dependent kinase inhibitor 2A (CDKN2A)-related ceRNA network in HCC.

Research objectives

To establish a CDKN2A-related ceRNA network in HCC and investigate the effect of the ceRNA on the prognosis of HCC.

Research methods

TCGA database, GEO database, cBioPortal database, lncBasev.3 database, miRDB database, Targetscan database, lncATLAS database, principal component analysis, The Gene Ontology enrichment, Kyoto Encyclopedia of Genes and Genomes pathway, meta-analysis, Kaplan-Meier plotter, nomogram, calibration curve, univariate and multivariate Cox regression, Spearman correlation coefficient, GSCALite database, UALCAN database, DiseaseMeth version 2.0 database, MEXPRESS database, TIMER database, CTD database, and PubChem database were used in this study.

Research results

CDKN2A was frequently mutated and deleted in and participated in cell cycle pathways in HCC. The CDKN2A-related ceRNA network- Growth arrest specific 5 (GAS5)/miR-25-3p/SRY-box transcription factor 11 (SOX11) was successfully established. GAS5 was recognized as an independent prognostic biomarker for HCC patients. GAS5 expression was correlated with methylation level, immune infiltration, and the expression of three immune checkpoint genes in HCC.

Research conclusions

The CDKN2A-related ceRNA network-GAS5/miR-25-3p/SOX11 was successfully established and GAS5 was an independent prognostic indicator for HCC patients.

Research perspectives

The CDKN2A-related ceRNA network provides innovative insights into the molecular mechanism of HCC formation and progression. Moreover, regulating GAS5 expression might be an effective therapeutic strategy in inducing immune infiltration and immunotherapy for HCC.