lncRNA-AC079061.1/VIPR1 axis may suppress the development of hepatocellular carcinoma: a bioinformatics analysis and experimental validation
Journal of Translational Medicine volume 20, Article number: 379 (2022)
Hepatocellular carcinoma (HCC) is one of the most malignant tumors to threaten human life, and the survival rate remains low due to delayed diagnosis. Meanwhile, lncRNAs have great potential for application in tumor prognosis, therefore relevant research in hepatocellular carcinoma is indispensable.
Based on the EZH2 expression, the differentially expressed lncRNAs DElncRNAs), miRNAs (DEmiRNAs), and mRNAs (DEmRNAs) were identified in hepatocellular carcinoma by using the TCGA database. Bioinformatics technology was utilized to determine the effect of key genes in HCC progression. The methylation and immune infiltration analyses were performed to explore the underlying function of hub genes. Finally, cellular function experiments were performed to investigate the association between identified genes and biological phenotypes in HCC.
lncRNA-AC079061.1, hsa-miR-765, and VIPR1 were identified as independent factors that affect the prognosis of hepatocellular carcinoma. The immune infiltration analyses revealed that lncRNA-AC079061.1 can alter the immune microenvironment and thus inhibit the development of HCC by regulating the expression of an immune-related gene (VIPR1). Methylation analyses demonstrated that VIPR1 expression is negatively related to the methylation level in HCC. Experimental results suggested that lncRNA-AC079061.1 and VIPR1 were frequently downregulated in HCC cells, while hsa-miR-765 was significantly upregulated. Moreover, the lncRNA-AC079061.1/VIPR1 axis suppressed the proliferation and invasion of HCC cells.
The present study identified the lncRNA-AC079061.1/VIPR1 axis as a novel biomarker that inhibited the proliferation and invasion of hepatocellular carcinoma, affecting the ultimate disease outcome.
Hepatocellular carcinoma (HCC) is the second leading cause of cancer-related mortality worldwide and is commonly seen in patients with chronic liver inflammation associated with viral infection or metabolic syndrome . It accounts for about 90% of non-metastatic liver tumors, which may result in systemic manifestations [2, 3]. Over the past decades, patients have benefited from advancements in new treatment strategies. However, the overall survival rate remains low due to delayed diagnosis of the diseases.
In recent years, the use of statistical models and bioinformatics in identifying novel biomarkers has been increasing . Therefore, such research provided great potential for clinical application. Various studies showed that lncRNAs play a crucial role in various diseases including cancer . For instance, Ma et al. found that EGFR1-mediated lnc01503 can promote gastric cancer progression . While Fan et al. reported that MKL1-induced lncRNA SNHG18 drives the growth and metastasis of non-small cell lung cancer .
Enhancer of zeste homolog 2 (EZH2) is a member of the polycomb group genes (PcGs) family, which is an important epigenetic regulator . It has been identified as a common oncogene in various types of tumors and is involved in many biological processes including cell proliferation , invasion [10, 11], metabolism [12, 13], and apoptosis . Besides, EZH2 has also been identified as an oncogene in HCC. Chen et al. verified that EZH2 can promote the development of HCC through modulating the miR-22/galectin-9 axis. Therefore, we speculate that the EZH2-related gene can affect the progression of HCC.
In the present study, the differentially expressed lncRNAs/miRNAs/mRNAs are screened based on the EHZ2 expression. Then, the prognostic genes were selected according to the overall survival curve. Multiple databases were used to analyze the possible interaction between lncRNA, miRNA, and mRNA. The corresponding statistical models (cox regression analysis and nomogram) were established to evaluate the survival prediction function of key genes. Methylation and immune infiltration analyses were used to determine the potential function of VIPR1. Finally, PCR and a series of cellular function experiments were performed to confirm the comprehensive analysis. The prediction of hub gene-related drugs can help identify potential treatment targets for HCC.
Materials and methods
The whole workflow chart of the present study was depicted in Fig. 1. The RNA-seq in TPM (transcripts per million reads) format, miRNA-seq in RPM (reads per million mapped reads) format, and clinical data of hepatocellular carcinoma cases were downloaded from The Cancer Genome Atlas (TCGA) portal (http://cancergenome.nih.gov). The whole cohort consisted of 374 HCC samples and 50 normal samples. The GTEX data were obtained from UCSC XENA (https://xenabrowser.net/datapages/) database. The gene annotation files were downloaded from the GENCODE database and used to extract the lncRNA and mRNA expression matrix, respectively. The immune gene list was obtained from the ImmPort Portal database (https://www.immport.org/home). The targets of lncRNA were predicted from the LncBase v.2 (http://carolina.imis.athena-innovation.gr/diana_tools/web/index.php?r=lncbasev2%2Findex), and the targets of miRNAs were predicted from the miRDB (http://mirdb.org/). The copy-number alteration data was downloaded from the cBioPortal database (http://www.cbioportal.org/), and the immunohistochemistry data were obtained from the Human Protein Atlas (HPA) database (https://www.proteinatlas.org/).
Differential expression analysis
The RNA-seq were divided into high and low groups according to the EZH2 expression levels. The difference analysis was performed using the “DESeq2” . For better and effective screening of interaction networks, lncRNAs with |logFC|> 2, adjusted p-value < 0.05, miRNAs with |logFC|> 0.5, adjusted p-value < 0.05, and mRNAs with |logFC|> 1, adjusted p-value < 0.05 were considered as differentially expressed genes. All the expression data were log2 transformed.
Since the expression data of RNA-seq and miRNA-seq did not conform to a normal distribution, we utilized the Spearman correlation to analyze the association between lncRNAs, miRNAs, and mRNAs. The results of the Spearman correlation were visualized using the “ggplot2” R package.
Functional enrichment analysis
To assess the potential mechanisms of the lncRNA-miRNA-mRNA axis, the KEGG and GO enrichment analyses were performed by using the “clusterProfiler” R package . Biological process terms, cellular component terms, molecular function terms, and pathways with adjusted p-value < 0.05 were considered as potential mechanisms for the onset of HCC.
Immune infiltration analysis
The markers of 24 immune cells were derived from an immunity article, and the classification and description of specific cells can be found in this article . The immune infiltration analysis was performed through the ssGSEA algorithm in the “GSVA” R package . The Pearson correlation analysis was used to evaluate the correlation between immune cells and differentially expressed genes.
The methylation of hub gene in HCC and normal tissues was evaluated through the MEXPRESS (https://mexpress.be/) and MethSurv database (https://biit.cs.ut.ee/methsurv/). The search terms were VIPR1 and liver hepatocellular carcinoma.
According to the clinical data of HCC, the effects of hub genes on patients’ prognoses were assessed through the Kaplan–Meier curve. The univariate Cox regression analysis was performed to screen for factors to construct the nomogram model. Nomograms were used to evaluate the role of hub genes in HCC through the “rms” and “survival” R packages. P-value < 0.05 was considered significantly.
Cell lines and cell culture
Human liver cells L0-2, HCC cells (PLC/PRF/5, HepG2, Hep3B) (Cell Bank of the Chinese Academy of Sciences, Shanghai, China), MHCC97H and HCCLM3 (Liver Cancer Institute, Fudan University, Shanghai, China), and Huh7 (Japanese Cancer Research Resources Bank) were cultured in Dulbecco's modified Eagle medium (Gibco) with 10% fetal bovine serum (Gibco) and 1% penicillin–streptomycin (Invitrogen). Cell cultures were done in a thermostatic incubator at 37 °C with a humidified atmosphere of 95% air and 5% CO2.
The lncRNA-AC079061.1 small interfering RNA (siRNA) and negative control were purchased from RiboBio Company (Guangzhou, China). The transfection of the lncRNA-AC079061.1 siRNA and negative control was performed with Lipofectamine™ 3000 (#L3000008, ThermoFisher, China) according to the manufacturer’s instructions.
RNA extraction and quantitative reverse-transcription polymerase chain reaction (qRT-PCR)
The total RNA including miRNA was extracted from cells or tissues using the TRIzol Reagent (Invitrogen), while cDNA was synthesized from RNA using the Reverse Transcription Kit (Takara). Subsequently, cDNA was amplified using the Maxinma SYBR Green qPCR Master Mix (Thermo Scientific). The quantification of target genes was done with the 2−ΔΔCt method using glyceraldehyde-3-phosphate dehydrogenase (GAPDH) for normalization. The primers for hsa-miR-765 and U6 small nuclear RNA were obtained from RiboBio Company (Guangzhou, China). The sequences were covered by a patent. Analyses of miRNA expression were normalized to the expression of internal control U6 using the 2 − ΔΔCT method. Melting curve analysis was carried out to assess the specificity of PCR products. The primers have been used for real time PCR, lncRNA-AC079061.1 Sequence 5′-3′: Forward, CCCTCAGGCATCCACTCTACCC; Reverse, TCCAACGCCACCCACTTCAAAC; VIPR1 Sequence 5′-3′: Forward, ACAAGGCAGCGAGTTTGGATGAG; Reverse, GTGCAGTGGAGCTTCCTGAACAG.
Luciferase Reporter Assay
MHCC97H cells and HCCLM3 cells were seeded into 96-well plates and co-transfected with a mixture of 60 ng luciferase, 6 ng pRL-CMV Renilla luciferase reporter, and miR-765 mimic or the negative control using the Lipofectamine 2000 transfection reagent (RiboBio, Guangzhou, China), according to the manufacturer’s instructions.). After 48 h of incubation, the firefly and Renilla luciferase activities were measured using a dual-luciferase reporter assay (Promega, Madison, WI, USA).
As described, the protein was extracted from cells or tissues using RIPA cell lysis with Protease Inhibitor Cocktail. Mitochondrial protein extraction was performed with the use of a mitochondrial isolation kit (Beyotime Biotechnology), in accordance with the manufacturer’s directions. The proteins were quantified by using the BCA kit, subjected to 12% SDS-PAGE for separation, then transferred to 0.45 μM PVDF membranes (Millipore, USA). The membranes were blocked with skimmed milk and incubated with primary antibodies at 4 ℃ overnight, followed by incubation with the corresponding HRP-conjugated secondary antibody (PeproTech). The bands were visualized by enhanced chemiluminescence. The intensity of protein expression was measured using ImageJ software. The antibodies used for western blot were listed as follows: VIPR1, Bioss, # bs-2982R, 1:1000; GAPDH, Cell Signaling Technology, #5174, 1:1000.
CCK-8 and colony formation assays
The cell proliferation assays were performed using a CCK-8 Kit (Yeasen, Shanghai, China) and colony formation. Three thousand cells were seeded into each well in a 96-well plate. The CCK-8 solution (10 µl) was added to 100 µl of culture media, and the optical density was measured at 450 nm. For colony formation assay, one thousand cells were seeded into each well in a 6-well plate for one week, then washed twice with PBS, fixed with 4% paraformaldehyde, stained with crystal violet, and the numbers of foci were counted for each well. Three independent experiments were performed.
Transwell migration and invasion assays
For migration assay, 5 × 104 cells were suspended in 200 µl of DMEM without serum and placed in the cell culture insert (8 µm pore size; BD Falcon, San Jose, CA) of a companion plate (BD Falcon) with a prewarmed culture medium containing 10% fetal bovine serum. For invasion assay, 1 × 105 cells suspended in serum-free medium were seeded into the upper chamber coated with 1 µg/µl Matrigel (BD Biosciences, USA) of a 24-well transwell plate (8-μm pore size, Corning, NY, USA), then 600 μl DMEM with 10% FBS was added into the lower chamber. After incubation for an indicated period of time at 37 °C in 5% CO2, the migrating and invading cells on the outside of the upper chamber membrane were then fixed with 4% paraformaldehyde, stained with crystal violet, and counted under a light microscope (100× magnification) in eight randomly selected areas. Three independent experiments were performed.
The expression of hub genes was assessed using the Wilcoxon Signed rank test and Mann–Whitney U test. Experimental data were expressed as mean ± standard deviations from three independent experiments and were analyzed using SPSS software (21.0; SPSS, Inc, Chicago, IL). Continuous variables between two groups or among three groups were compared using the unpaired Student’s t-test or one-way analysis of variance (Bonferroni post hoc test) as appropriate. Categorical variables were compared using the Chi-square test or Fisher’s exact test as appropriate. All statistical tests were two-sided and a P-value < 0.05 was considered statistically significant.
The role of EZH2 in HCC
Previous studies have reported that EZH2 plays a crucial role in HCC. We hereby comprehensively discussed the important function of EZH2 in HCC through bioinformatics analysis [19,20,21]. The expression of EZH2 is frequently upregulated in HCC tissues compared to normal tissues (Fig. 2A). The related survival analysis of EZH2 showed that patients with high EZH2 expression have a worse overall survival (OS: HR = 1.99; 95% CI: 1.40–2.84; p < 0.001), disease-specific survival (DFS: HR = 2.44; 95% CI: 1.53–3.89; p < 0.001), and progression-free interval (PFI: HR = 1.81; 95% CI: 1.35–2.42; p < 0.001) (Fig. 2B). The association between EZH2 mRNA expression and copy number is shown in Fig. 2C and D), which suggested that the mRNA expression level of EZH2 is positively related to the copy number of EZH2 in HCC (R = 0.15, p = 4.227e−3). The protein expression of EZH2 in HCC and normal tissues were further evaluated through the HPA database. The immunohistochemistry of EZH2 showed that EZH2 is mainly located in the nuclear and highly expressed in tumor tissues compared with normal tissues (Fig. 2E). These results identified EZH2 as a prognostic factor, which promoted the progression of HCC. Future studies based on EZH2 may help explore more novel biomarkers and treatment strategies in HCC.
Screening of DElncRNAs, DEmiRNAs, and DEmRNAs
Based on the expression of EZH2, the whole HCC cohort was divided into EZH2high and EZH2Low groups using median partitioning. The difference analysis was performed to obtain the differentially expressed lncRNAs (DElncRNAs). lncRNAs with |logFC|> 2 and p-value < 0.05 was considered as DElncRNAs. Then, 163 upregulated and 31 downregulated lncRNAs were visualized through the volcano map (Fig. 3A). The top 10 significant downregulated lncRNAs in order of p-value were selected for gene co-expression analysis, and the result showed that there is a strong negative correlation between 9 lncRNAs (LINC00570, AC079061.1, HTR2A-AS1, LDLRAD4-AS1, AC022601.1, PGM5-AS1, FAM83A-AS1, AL391261.2, and TMEM252-DT) and EZH2 (Fig. 3B). We assessed the expression level of 10 lncRNAs between HCC tissues and normal tissues, the result suggested that 5 lncRNAs (AC079061.1, HTR2A-AS1, LDLRAD4-AS1, FAM83A-AS1, and AL391261.2) were highly expressed in normal tissues, while 1 lncRNAs (AC005165.1) was highly expressed in tumor tissues (Fig. 3C). Then, ROC analyses were performed to assess the diagnostic value of these 10 lncRNAs. The results showed that only FAM83A-AS1, HTR2A-AS1, AC079061.1, and LDLRAD4-AS1 have a better diagnostic value in HCC, the area under curve (AUC) of these four lncRNAs were 0.709, 0.853, 0.804, and 0.712, respectively (Fig. 3D). Furthermore, survival analyses were used to evaluate the impact of lncRNAs on the outcome of HCC. Finally, the results suggested that AC079061.1 is the only lncRNA that can affect the prognosis of HCC. Patients with high AC079061.1 expression have a better outcome compared with patients with low AC079061.1 expression (HR = 0.63; 95% CI: 0.44–0.89; p = 0.01) (Fig. 3E).
Based on the ceRNA mechanism, potential miRNAs that may bind to DElncRNAs were first predicted using the LncBase v.2 database (http://carolina.imis.athena-innovation.gr/diana_tools/web/index.php?r=lncbasev2%2Findex-predicted). Meanwhile, the differentially expressed miRNAs between HCC tissues and normal tissues were obtained by the “limma” R packages. Finally, 8 miRNAs that may bind to DElncRNAs were identified in HCC by intersecting the LncBase predicted miRNAs with the DEmiRNAs in TCGA (Fig. 4A and B). The expression levels of these 8 miRNAs were confirmed in Fig. 4C, which showed that 2 miRNAs (hsa-miR-4725-3p and hsa-miR-6514-3p) were downregulated in tumor tissues, and 6 miRNAs (hsa-miR-939-5p, hsa-miR-3690, hsa-miR-765, hsa-miR-362-3p, hsa-miR-766-5p, and hsa-miR-4742-3p) were elevated in tumor tissues. Next, ROC analyses were performed to evaluate the diagnostic value of DEmiRNAs. The results suggested that hsa-miR-362-3p (AUC = 0.755), hsa-miR-765 (AUC = 0.720), hsa-miR-939-5p (AUC = 0.713), and hsa-miR-4742-3p (AUC = 0.717) have a better diagnostic value than hsa-miR-766-5p (AUC = 0.596), hsa-miR-3690 (AUC = 0.604), hsa-miR-4725-3p (AUC = 0.588), and hsa-miR-6514-3p (AUC = 0.567) (Fig. 4D). Furthermore, survival analyses were utilized to assess the prognostic value of DEmiRNAs. Results showed that only 2 miRNAs (hsa-miR-3690 (HR = 1.47, p = 0.03) and hsa-miR-765 (HR = 1.43, p = 0.046) affected the prognosis of HCC (Fig. 4E). Collectively, AC079061.1 may regulate the progression of HCC by binding to hsa-miR-765.
In recent years, tumor immune microenvironment has been the focus of research and a hot topic. Many studies have suggested an inseparable relationship between tumorigenesis and the level of immune infiltration [22,23,24,25]. Hence, we report 6 immune-related genes that may bind to hsa-miR-765 by intersecting the miRDB-predicted miRNAs and immune-related gene list with DEmRNA (Fig. 5A and B). The expression level of differentially expressed immune-related genes (DEIRGs) was confirmed in Fig. 5C. The ROC and survival analysis results of 6 DEIRGs suggested that only VIPR1 (AUC = 0.975, HR = 0.67, p = 0.023) has both diagnostic and prognostic values in HCC (Fig. 5D and E).
Correlation between lncRNA, miRNA, and mRNA
To further explore the association between lncRNA- AC079061.1, hsa-miR-765, and VIPR1, we first analyzed the subcellular localization of the lncRNA- AC079061.1. The sequence of lncRNA-AC079061.1 was obtained from the LNCipedia database (https://lncipedia.org/db/transcript/lnc-EBAG9-9:1), and the sequence was imported into the lncLocator database to predict the subcellular location of lncRNA- AC079061.1. The results suggested that lncRNA- AC079061.1 was mainly located in the ribosome and cytosol (Fig. 6A). The Spearman correlation analysis was used to assess the correlation between lncRNA- AC079061.1, hsa-miR-765, and VIPR1 since the expression of these genes did not conform to normal distribution. Then, results showed a negative correlation between hsa-miR-765 and lncRNA- AC079061.1 or VIPR1. A positive correlation was also noted between lncRNA- AC079061.1 and VIPR1 (Fig. 6B–D).
Immune infiltration analysis of lncRNA, miRNA, and mRNA
Previous studies confirmed that immune cells were involved in the progression of HCC [26,27,28]. Thus, we explored the correlation between immune infiltration and these key genes. Initial results suggested that lncRNA-AC079061.1 was negatively correlated with NK CD56 bright cells, TFH, Th2 cells, macrophages, B cells, and Th1 cells, while positively correlated with neutrophils, Tcm, CD8 T cells (Fig. 7A). Subseqent results showed that hsa-miR-765 was negatively correlated with DC cells, cytotoxic cells, neutrophils, CD8 T cells, NK cells, NK CD56dim cells, and Treg, while positively correlated with NK CD56bright cells (Fig. 7B). VIPR1, an immune-related gene, was positively correlated with many immune cells such as NK cells, neutrophils, mast cells, CD8 T cells, and eosinophils, while negatively correlated with Th2 cells and pDC cells (Fig. 7C). Furthermore, we also analyzed the correlation between VIP or ADCYAP1 (VIPR1’s ligand) and the infiltration level of immune cells. Results showed that the VIP and ADCYAP1 were similar to VIPR1, which displayed a positive correlation with the infiltration level of NK cells, iDC, and TEM (Additional file 1: Fig. S1A–D), suggesting the possibility that altered expression of ligands and receptors may lead to the change in immune cell infiltration level and has a dramatic impact on HCC progression. Interestingly, we found that lncRNA AC079061.1, hsa-miR-765, and VIPR1 were all related to neutrophils and CD8 T cells. The correlation analysis showed that there is a positive relationship between lncRNA-AC079061.1/VIPR1 and CD8 T cells/Neutrophils (Fig. 7D and E). These results implied that the lncRNA-AC079061.1/hsa-miR-765/VIPR1 axis may influence the progression of HCC by regulating the immune infiltration level of neutrophils and CD8 T cells.
Functional enrichment analysis of VIPR1
The potential mechanism of VIPR1 remains unclear and was therefore uploaded to the String database. Then, the interaction network between VIPR1 and related genes was obtained (Fig. 8A). The KEGG and GO functional enrichment analyses were performed to explore the potential VIPR1-related mechanism. The results showed that VIPR1 and its related genes were highly enriched in neuroactive ligand-receptor interaction (adjusted p-value = 1.38e−04), cAMP signaling pathway(adjusted p-value = 0.004), and vascular smooth muscle contraction (adjusted p-value = 0.011) (Fig. 8B). Subsequently, the most highly enriched GO BP terms were cyclic-nucleotide-mediated signaling (adjusted p-value = 2.38e−09), G protein-coupled receptor signaling pathway, coupled to cyclic nucleotide second messenger (adjusted p-value = 5.42e−09), and adenylate cyclase-activating G protein-coupled receptor signaling pathway (adjusted p-value = 1.05e−08). The GO CC terms showed that these genes may be located in the perikaryon (adjusted p-value = 0.068), endocytic vesicle membrane (adjusted p-value = 0.068), and Mitochondria-associated ER Membrane (adjusted p-value = 0.068). Furthermore, the highly enriched GO MF terms were peptide hormone binding (adjusted p-value = 4.89e−07), G protein-coupled peptide receptor activity (adjusted p-value = 9.55e−06), and peptide receptor activity (adjusted p-value = 9.55e−06) (Fig. 8C, Table 1).
To further explore the possible mechanisms of VIPR1, the top 100 VIPR1-related genes were downloaded from the GEPIA database (http://gepia.cancer-pku.cn/detail.php?gene=VIPR1) for another functional enrichment analysis (Additional file 1: Table S1.). More physiological processes such as complement activation, lectin pathway (adjusted p-value = 2.33e−04), positive regulation of vasculature development (adjusted p-value = 2.33e−04), cytokine receptor binding (adjusted p-value = 1.99e−04), receptor-ligand activity (adjusted p-value = 0.027) were obtained, suggesting that VIPR1 and its related genes play a key role in multiple mechanisms including the immune process. In summary, the results of these functional enrichment analyses suggested that VIPR1, as an immune-related gene, is involved in a variety of physiological processes (such as cytokine binding, inflammation-stimulation, and stress response) by regulating relevant genes, indicating the tremendous research value of VIPR1 in HCC.
Methylation analysis of VIPR1 in HCC
Methylation is known to be aberrantly expressed in various forms of tumors. Meanwhile, it is capable of affecting not only the genetic phenotypes of cells but also the life span of humankind. There is a great potential for methylation in tumor diagnosis and treatments [29,30,31]. Therefore, multiple methods, as well as databases, were used to assess the correlation between VIPR1 and the progression of HCC. We first analyzed the correlation between VIPR1 expression and methyltransferases in HCC. As shown in Additional file 1: Fig. S2A, the results showed that there is a negative correlation between VIPR1 and DNMT1 (R = −0.150, p = 0.003) /DNMT3A (R = −0.250, p < 0.001). But no significant correlation was noted between VIPR1 and DNMT3B (R = −0.099, p = 0.057). Similarly, DiseaseMeth 2.0 database demonstrated that the methylation of VIPR1 was significantly lower in disease tissues compared with HCC normal tissues (p = 0.021, Additional file 1: Fig. S2B). Additionally, the CpG expression level in Additional file 1: Fig. S2C was obtained from the MethSurv platform, in which 8 methylation sites (cg27510539, cg12828331, cg23003500, cg10970409, cg04322483, cg22694480, cg27309542, cg23517013) in the DNA sequences of VIPR1 were found to be negatively associated with VIPR1 expression levels and 11 methylation sites (cg06783423, cg03791150, cg00328740, cg18559858, cg02061253, cg09475741, cg11484444, cg09996971, cg01919863, cg15066503, cg06290884) were found to be positively associated with VIPR1 expression level, analyzed through the MEXPRESS tool (Additional file 1: Fig. S2D). According to these results, altered methylation of VIPR1 may be an inducer of HCC.
Prognostic analysis in HCC
To assess the potential prognostic value of lncRNA AC079061.1, hsa-miR-765, and VIPR1 in HCC, univariate and multivariate cox regression analyses were performed. The results showed that lncRNA AC079061.1, hsa-miR-765, and VIPR1 as independent factors can affect the prognosis of HCC (Fig. 9A–C; Table 2, 3, 4). Then, we constructed prediction models to evaluate the relationship between lncRNA AC079061.1/hsa-miR-765/VIPR1 and the overall survival of HCC. The models were visualized by nomogram (Fig. 9D, E and F). All the results implied that the lncRNA AC079061.1/hsa-miR-765/VIPR1 axis acts as a prognostic indicator that may suppress the development of HCC.
Validation of lncRNA-AC079061.1/VIPR1 axis
To validate the feasibility of our comprehensive analysis, the expressions of these three genes were detected in the HCC cell lines and normal liver cells. We found that lncRNA and mRNA were frequently downregulated in HCC cells (Figs. 10A and B), while hsa-miR-765 was upregulated significantly in HCC cells (Fig. 10C). Then, we transfected the lncRNA- AC079061.1 siRNA into HCC cells, and the qRT-PCR was performed to confirm the stable downregulation of lncRNA- AC079061.1 (Additional file 1: Fig. S3). Moreover, to evaluate the relationship between lncRNA- AC079061.1 and the malignant phenotype of the tumor, a series of cellular function experiments was performed. The CCK-8 and colony formation assays suggested that the knockdown of lncRNA-AC079061.1 promoted the proliferation of HCC cells (Fig. 10D and E). Then, the transwell assays showed that the invasion and migration ability of si-lncRNA-AC079061.1 HCC cells were enhanced (Fig. 10F). These results indicated that lncRNA- AC079061.1 functions as an antioncogenic lncRNA to inhibit the proliferation, migration, and invasion of HCC cells in vitro. Additionally, to further verify the lncRNA-AC079061.1/VIPR1 axis, we first detected the expression level of hsa-miR-765 and VIPR1 in si- lncRNA- AC079061.1 HCC cells. Results showed that hsa-miR-765 is elevated in si- lncRNA- AC079061.1 cells compared with control cells (Fig. 10G), while the mRNA (Fig. 10G) and protein expression (Fig. 10H) of VIPR1 was decreased as the downregulation of lncRNA- AC079061.1 or upregulation of miR-765. Then, we performed dual-luciferase reporter assays in HCC cells and found that the luciferase activity was significantly downregulated after increasing the hsa-miR-765 expression in MHCC97H and HCCLM3 cells transfected with lncRNA-AC079061.1 or VIPR1 wildtype vectors. These results suggested that lncRNA-AC079061.1 may regulate the expression of downstream gene VIPR1 through sponging hsa-miR-765 (Fig. 10I). Furthermore, we investigated the effects of miR-765 on lncRNA-AC079061.1-mediated antioncogenic processes in vitro. Upregulation of miR-765 inhibited si-lncRNA-AC079061.1 induced proliferation (Fig. 10J), migration, and invasion (Fig. 10K) in MHCC97H cells. Collectively, these results indicated that the lncRNA- AC079061.1/hsa-miR-765/VIPR1 axis may regulate the progression of HCC.
Prediction of chemotherapeutic drugs
Clinically, chemotherapy remains one of the mainstream treatments for HCC. However, in recent years, HCC patients tend to develop tolerance to common chemotherapy drugs. Therefore, the exploration of novel drugs is currently an essential trend in research. By utilizing the GDSC and CTRP databases, we obtained the correlation between VIPR1 expression and drug sensitivity in pan-cancer. These results provided good reference for future research direction. The results showed that the key genes investigated in the present study did not only correlate with common drugs such as Gefitinib and Afatinib but also with some unreported drugs. These correlations may even help identify new therapeutic drugs for HCC (Additional file 1: Fig. S4).
Cancer-related mortality remains the leading cause of death and a major health burden in Asia. Approximately, 49.3% of new cancer cases are located in Asia. Interestingly, more than half of these cases were reported in China, indicating that tumors have always been a major health concern. According to the global cancer statistics, the number of new cases and deaths of liver cancer ranks 7th among all cancers in the United States . Furthermore, regardless of gender, liver cancer has always been ranked top 10 leading causes of death . The understanding of underlying molecular mechanisms in HCC is of great importance, especially in the early diagnosis and treatment of the disease. In the present study, bioinformatics analyses combined with experimental validation identified novel biomarkers for the diagnosis and treatment of HCC.
In recent years, the continuous improvements in bioinformatics have been well applied to the discovery of new research directions, new research targets, and therapeutic drugs [34,35,36]. Xu et al. identified several immune-related lncRNA signatures in HCC through bioinformatics technology . Liang et al. found that methylation-related genes (CTF1, FZD8, PDK4, and ZNF334) affected the progression of HCC . The present study identified the lncRNA- AC079061.1/hsa-miR-765/VIPR1 axis, providing great research value in HCC by employing bioinformatics technology.
In the present study, an EZH2-related ceRNA network that is associated with the prognosis of HCC was constructed. Previous studies confirmed EZH2 as an oncogene that accelerates the development of HCC and negatively regulates the expression of immune checkpoint inhibitor PD-L1 in HCC [19, 20, 39]. Based on the function of EZH2, we screened for relevant differentially expressed lncRNAs, miRNAs, and mRNAs in HCC. Among all the DEGs, lncRNA- AC079061.1, hsa-miR-765, and VIPR1 potentially regulate the progression of HCC. hsa-miR-765 may bind to the 3’-UTR pf lncRNA- AC079061.1 and VIPR1, resulting in mutual regulation.
Interestingly, lncRNA- AC079061.1 has not been reported or studied before in HCC. Therefore, our study is the first to identify that this lncRNA has a function that suppresses HCC progression and the malignant phenotype of HCC. We also found that lncRNA- AC079061.1 is closely related to the prognosis of HCC.
In addition, it was reported that hsa-miR-765 is involved in various types of malignant tumors. Xie et al. verified that miR-765 promoted cell proliferation in HCC . While Zhu et al. validated that LINC00994 repressed the malignant behaviors of pancreatic cancer cells through regulating the miR-765-3p/RUNX2 axis (which means that miR-765 accelerated the development of HCC) . The oncogenic effect of miR-765 is not only limited to these two tumors, but also plays a similar role in esophageal squamous cell carcinoma , gastric cancer , and osteosarcoma .
For VIPR1, there are few relevant studies in HCC. However, low expression of VIPR1 has been shown to have an adverse prognostic impact on HCC . However, it is not enough to study the role of VIPR1 in HCC alone, Zhao et al. also pointed out that VIPR1 has been confirmed to play the same role in lung adenocarcinoma . These studies revealed the good theoretical feasibility of VIPR1 in HCC.Then, DNA methylation is known to regulate gene transcription and silence tumor suppressors . A previous study reported that the H3K27 deacetylation and promoter methylation affect the expression of VIPR1 . Altogether, these findings suggested the antitumor effect of VIPR1 in HCC.
Combining the results of previous studies and corresponding prediction databases, lncRNA- AC079061.1, hsa-miR-765, and VIPR1 were collectively analyzed to investigate the influence on the outcome of HCC either as an independent prognostic factors or as a lncRNA- AC079061.1/hsa-miR-765/VIPR1 axis.
In summary, the lncRNA- AC079061.1/hsa-miR-765/VIPR1 axis may be a novel biomarker in HCC through comprehensive analyses and experimental validation. The prognostic analyses also suggest that it may suppress the development of HCC. Conclusively, the novel axis has a great potential to explore the pathogenesis of HCC and provide future directions for new therapeutic targets and drugs.
Availability of data and materials
All data generated or analyzed during this study are included in this article.
The Cancer Genome Atlas
- DElncRNAs, DEmiRNAs, DEmRNAs:
Differentially expressed lncRNAs/miRNAs/mRNAs
Kyoto Encyclopedia of Genes and Genomes
Tang W, Chen Z, Zhang W, Cheng Y, Zhang B, Wu F, Wang Q, Wang S, Rong D, Reiter FP, et al. The mechanisms of sorafenib resistance in hepatocellular carcinoma: theoretical basis and therapeutic aspects. Signal Transduct Target Ther. 2020;5:87.
Huang A, Yang XR, Chung WY, Dennison AR, Zhou J. Targeted therapy for hepatocellular carcinoma. Signal Transduct Target Ther. 2020;5:146.
Sangro B, Chan SL, Meyer T, Reig M, El-Khoueiry A, Galle PR. Diagnosis and management of toxicities of immune checkpoint inhibitors in hepatocellular carcinoma. J Hepatol. 2020;72:320–41.
Iasonos A, Chapman PB, Satagopan JM. Quantifying treatment benefit in molecular subgroups to assess a predictive biomarker. Clin Cancer Res. 2016;22:2114–20.
Bhan A, Soleimani M, Mandal SS. Long noncoding RNA and cancer: a new paradigm. Cancer Res. 2017;77:3965–81.
Ma Z, Gao X, Shuai Y, Wu X, Yan Y, Xing X, Ji J. EGR1-mediated linc01503 promotes cell cycle progression and tumorigenesis in gastric cancer. Cell Prolif. 2021;54: e12922.
Fan H, Yuan J, Li Y, Jia Y, Li J, Wang X, Li X. MKL1-induced lncRNA SNHG18 drives the growth and metastasis of non-small cell lung cancer via the miR-211-5p/BRD4 axis. Cell Death Dis. 2021;12:128.
Duan R, Du W, Guo W. EZH2: a novel target for cancer treatment. J Hematol Oncol. 2020;13:104.
Varambally S, Dhanasekaran SM, Zhou M, Barrette TR, Kumar-Sinha C, Sanda MG, Ghosh D, Pienta KJ, Sewalt RG, Otte AP, et al. The polycomb group protein EZH2 is involved in progression of prostate cancer. Nature. 2002;419:624–9.
Kleer CG, Cao Q, Varambally S, Shen R, Ota I, Tomlins SA, Ghosh D, Sewalt RG, Otte AP, Hayes DF, et al. EZH2 is a marker of aggressive breast cancer and promotes neoplastic transformation of breast epithelial cells. Proc Natl Acad Sci USA. 2003;100:11606–11.
Xia L, Zhu X, Zhang L, Xu Y, Chen G, Luo J. EZH2 enhances expression of CCL5 to promote recruitment of macrophages and invasion in lung cancer. Biotechnol Appl Biochem. 2020;67:1011–9.
Yiew NKH, Greenway C, Zarzour A, Ahmadieh S, Goo B, Kim D, Benson TW, Ogbi M, Tang YL, Chen W, et al. Enhancer of zeste homolog 2 (EZH2) regulates adipocyte lipid metabolism independent of adipogenic differentiation: role of apolipoprotein E. J Biol Chem. 2019;294:8577–91.
Ahmad F, Patrick S, Sheikh T, Sharma V, Pathak P, Malgulwar PB, Kumar A, Joshi SD, Sarkar C, Sen E. Telomerase reverse transcriptase (TERT) - enhancer of zeste homolog 2 (EZH2) network regulates lipid metabolism and DNA damage responses in glioblastoma. J Neurochem. 2017;143:671–83.
Wu ZL, Zheng SS, Li ZM, Qiao YY, Aau MY, Yu Q. Polycomb protein EZH2 regulates E2F1-dependent apoptosis through epigenetically modulating Bim expression. Cell Death Differ. 2010;17:801–10.
Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15:550.
Yu G, Wang LG, Han Y, He QY. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS. 2012;16:284–7.
Bindea G, Mlecnik B, Tosolini M, Kirilovsky A, Waldner M, Obenauf AC, Angell H, Fredriksen T, Lafontaine L, Berger A, et al. Spatiotemporal dynamics of intratumoral immune cells reveal the immune landscape in human cancer. Immunity. 2013;39:782–95.
Hänzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinformatics. 2013;14:7.
Xiao G, Jin LL, Liu CQ, Wang YC, Meng YM, Zhou ZG, Chen J, Yu XJ, Zhang YJ, Xu J, Zheng L. EZH2 negatively regulates PD-L1 expression in hepatocellular carcinoma. J Immunother Cancer. 2019;7:300.
Chen S, Pu J, Bai J, Yin Y, Wu K, Wang J, Shuai X, Gao J, Tao K, Wang G, Li H. EZH2 promotes hepatocellular carcinoma progression through modulating miR-22/galectin-9 axis. J Exp Clin Cancer Res. 2018;37:3.
Zhao H, Xu Y, Mao Y, Zhang Y. Effects of EZH2 gene on the growth and migration of hepatocellular carcinoma HepG2 cells. Hepatobiliary Surg Nutr. 2013;2:78–83.
Xu Q, Chen S, Hu Y, Huang W. Landscape of immune microenvironment under immune cell infiltration pattern in breast cancer. Front Immunol. 2021;12: 711433.
Dumauthioz N, Labiano S, Romero P. Tumor resident memory T Cells: new players in immune surveillance and therapy. Front Immunol. 2018;9:2076.
Efimova I, Catanzaro E, Van der Meeren L, Turubanova VD, Hammad H, Mishchenko TA, Vedunova MV, Fimognari C, Bachert C, Coppieters F, et al. Vaccination with early ferroptotic cancer cells induces efficient antitumor immunity. J Immunother Cancer. 2020;8:e001369.
Fujii SI, Shimizu K. Immune networks and therapeutic targeting of iNKT cells in cancer. Trends Immunol. 2019;40:984–97.
Ma J, Zheng B, Goswami S, Meng L, Zhang D, Cao C, Li T, Zhu F, Ma L, Zhang Z, et al. PD1(Hi) CD8(+) T cells correlate with exhausted signature and poor clinical outcome in hepatocellular carcinoma. J Immunother Cancer. 2019;7:331.
Wu X, Luo H, Shi B, Di S, Sun R, Su J, Liu Y, Li H, Jiang H, Li Z. Combined antitumor effects of sorafenib and GPC3-CAR T cells in mouse models of hepatocellular carcinoma. Mol Ther. 2019;27:1483–94.
Di Blasi D, Boldanova T, Mori L, Terracciano L, Heim MH, De Libero G. Unique T-cell populations define immune-inflamed hepatocellular carcinoma. Cell Mol Gastroenterol Hepatol. 2020;9:195–218.
Klutstein M, Nejman D, Greenfield R, Cedar H. DNA methylation in cancer and aging. Cancer Res. 2016;76:3446–50.
Morgan AE, Davies TJ, Mc Auley MT. The role of DNA methylation in ageing and cancer. Proc Nutr Soc. 2018;77:412–22.
Pan Y, Liu G, Zhou F, Su B, Li Y. DNA methylation profiles in cancer diagnosis and therapeutics. Clin Exp Med. 2018;18:1–14.
Sung H, Ferlay J, Siegel RL, Laversanne M, Soerjomataram I, Jemal A, Bray F. Global cancer statistics 2020: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J Clin. 2021;71:209–49.
Siegel RL, Miller KD, Fuchs HE, Jemal A. Cancer statistics, 2022. CA Cancer J Clin. 2022;72:7–33.
Mahar M, Cavalli V. Intrinsic mechanisms of neuronal axon regeneration. Nat Rev Neurosci. 2018;19:323–37.
Zhang Z, Zhou L, Xie N, Nice EC, Zhang T, Cui Y, Huang C. Overcoming cancer therapeutic bottleneck by drug repurposing. Signal Transduct Target Ther. 2020;5:113.
Browne CJ, Godino A, Salery M, Nestler EJ. Epigenetic mechanisms of opioid addiction. Biol Psychiatry. 2020;87:22–33.
Xu Q, Wang Y, Huang W. Identification of immune-related lncRNA signature for predicting immune checkpoint blockade and prognosis in hepatocellular carcinoma. Int Immunopharmacol. 2021;92: 107333.
Liang Y, Ma B, Jiang P, Yang HM. Identification of methylation-regulated differentially expressed genes and related pathways in hepatocellular carcinoma: a study based on TCGA database and bioinformatics analysis. Front Oncol. 2021;11: 636093.
Gao SB, Xu B, Ding LH, Zheng QL, Zhang L, Zheng QF, Li SH, Feng ZJ, Wei J, Yin ZY, et al. The functional and mechanistic relatedness of EZH2 and menin in hepatocellular carcinoma. J Hepatol. 2014;61:832–9.
Xie BH, He X, Hua RX, Zhang B, Tan GS, Xiong SQ, Liu LS, Chen W, Yang JY, Wang XN, Li HP. Mir-765 promotes cell proliferation by downregulating INPP4B expression in human hepatocellular carcinoma. Cancer Biomark. 2016;16:405–13.
Zhu X, Niu X, Ge C. Inhibition of LINC00994 represses malignant behaviors of pancreatic cancer cells: interacting with miR-765-3p/RUNX2 axis. Cancer Biol Ther. 2019;20:799–811.
Jiang B, Xu G, Lv HQ, Huang M, Li Z. Up-regulation of miR-765 predicts a poor prognosis in patients with esophageal squamous cell carcinoma. Eur Rev Med Pharmacol Sci. 2018;22:3789–94.
Lin W, Miao Y, Meng X, Huang Y, Zhao W, Ruan J. miRNA-765 mediates multidrug resistance via targeting BATF2 in gastric cancer cells. FEBS Open Bio. 2020;10:1021–30.
Lv DB, Zhang JY, Gao K, Yu ZH, Sheng WC, Yang G, Gao YZ. MicroRNA-765 targets MTUS1 to promote the progression of osteosarcoma via mediating ERK/EMT pathway. Eur Rev Med Pharmacol Sci. 2019;23:4618–28.
Lu S, Lu H, Jin R, Mo Z. Promoter methylation and H3K27 deacetylation regulate the transcription of VIPR1 in hepatocellular carcinoma. Biochem Biophys Res Commun. 2019;509:301–5.
Zhao L, Yu Z, Zhao B. Mechanism of VIPR1 gene regulating human lung adenocarcinoma H1299 cells. Med Oncol. 2019;36:91.
Dor Y, Cedar H. Principles of DNA methylation and their implications for biology and medicine. Lancet. 2018;392:777–86.
We acknowledge TCGA and other public databases for providing their platforms and contributors for uploading their meaningful datasets.
This work was supported by the National Natural Science Foundation of China (grant no. 8197100979) and Natural Science Foundation of Autonomous Region (grant no. XZ 2019 ZR G-117).
Ethics approval and consent to participate
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
: Table S1. The functional enrichment analysis of 100 VIPR1-related genes. Figure S1. Immune infiltration analysis of VIP and ADCYAP1 A. C. Explored the correlation between immune infiltration and VIP and ADCYAP1 through the ssGSEA algorithm in the “GSVA” R package. B. D. Evaluated the correlation between Treg/Tem/NK/iDC cells and VIP, Tem/NK/iDC/Mast cells and ADCYAP by the Pearson correlation analysis. *P<0.05, **P<0.01, ***P<0.001. Figure S2. The methylation analysis of VIPR1 in HCC. A. The correlation between VIPR1 expression and methyltransferases in HCC. B. The methylation level of VIPR1 in disease tissues and normal tissues by DiseaseMeth 2.0 database. C. The CpG expression level of VIPR1 in HCC. D. The methylation site of VIPR1 DNA sequence association with gene expression was visualized using MEXPRESS. The expression of VIPR1 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. Figure S3. Confirmation of the stable downregulation of lncRNAAC079061.1 transfected with lncRNA- AC079061.1 siRNA into MHCC97H and HCCLM3 HCC cells by qRT-PCR. *P<0.05, **P<0.01, ***P<0.001. Figure S4. The correlation between drugs and the VIPR1 mRNA expression of HCC in CTRP and GDSC database
About this article
Cite this article
Lin, XH., Zhang, DY., Liu, ZY. et al. lncRNA-AC079061.1/VIPR1 axis may suppress the development of hepatocellular carcinoma: a bioinformatics analysis and experimental validation. J Transl Med 20, 379 (2022). https://doi.org/10.1186/s12967-022-03573-7