Comprehensive analysis to identify noncoding RNAs mediated upregulation of maternal embryonic leucine zipper kinase (MELK) correlated with poor prognosis in hepatocellular carcinoma

Expression and survival analysis of MELK in the pan-cancer dataset

Data on MELK mRNA expression in 18 cancers were obtained from TCGA (https://genome-cancer.ucsc.edu/). The data were normalized and the package “limma” in R was used to analyze differential expression [13]. A p-value < 0.05 was taken as statistically significant. The GEPIA database (http://gepia.cancer-pku.cn/) [14] was used to assess MELK and lncRNA levels in various cancers, again using p < 0.05 to represent significance.

Prediction of MELK-binding miRNAs

We used PITA, RNA22, miRmap, microT, miRanda, PicTar, and TargetScan to predict the potential binding of miRNAs. Candidate miRNAs predicted by two or more programs were used for further analysis. StarBase (http://starbase.sysu.edu.cn/), a database for analyzing miRNAs and their interactions, [15] was used to investigate relationships between miRNA-MELK, lncRNA- MIR3142HG, or lncRNA- LINC00265 in HCC. The levels of hsa-miR-101-3p in tumor and control tissue were determined, as were potential lncRNA candidates for binding to has-miR-101-3p.

MELK expression in HCC

Gene expression patterns and clinical data pertaining to HCC patients were downloaded from the TCGA database (https://portal.gdc.cancer.gov). RNA-seq data from 375 HCC patients were obtained. MELK expression patterns in HCC were also investigated in the Oncomine dataset and International Cancer Genome Consortium LIRI-JP cohort (n = 229) (https://dcc.icgc.org/projects/LIRI-JP). The Human Protein Atlas (HPA, http://www.proteinatlas.org/) was searched for information on MELK protein expression, and information on mutations was obtained from the cBioPortal for Cancer Genomics (http://www.cbioportal.org/). The institutional ethics committee of the First Affiliated Hospital of Jinzhou Medical University approved this study.

Survival analyses and prognostic model development

Kaplan-Meier curves and log-rank tests were used to measure survival and to assess the prognostic relevance of MELK. Univariate Cox regression was conducted to explore the associations between MELK expression and OS to identify prognostic biomarkers. Multivariate Cox regression was subsequently used to identify factors that were independently associated with outcomes.

DNA methylation analyses

DNA methylation is controlled by DNA methyltransferases and influences cancer cell behavior. We investigated the expression of DNA methyltransferases in relation to MELK expression in data from the TCGA database. The UALCAN (http://ualcan.path.uab.edu/) and DiseaseMeth v 2.0 (http://bio-bigdata.hrbmu.edu.cn/diseasemeth/) databases were then used to examine MELK expression in HCC tumors and paracancerous tissues and MEXPRESS (https://mexpress.be) [16] was used to determine relationships between MELK and DNA methylation.

Immune cell infiltration in relation to MELK levels

TIMER (https://cistrome.shinyapps.io/timer/) [17] was used to investigate tumor immune infiltration, specifically, to determine whether there was an association between MELK expression and the infiltration by different types of immune cells. Correlations between MELK levels and the expression of immune checkpoint genes, specifically, CD274, CTLA4, HAVCR2, PDCD1, PDCD1LG2, TIGIT, LAG3, and SIGLEC15, were explored. In addition, the relationships between these variables and MELK copy numbers were assessed in HCC patients, as was their prognostic relevance. Moreover, the correlations between MELK expression and markers associated with 16 different tumor-infiltrating lymphocytes (TILs) were assessed; these included B cells, monocytes, T cells, CD8+ T cells, neutrophils, M1/M2 macrophages, natural killer (NK) cells, DCs, exhausted T cells, Tregs, tumor-associated macrophages (TAMs), and Th1, Th2, Th17, and Tfh cells. Outcome modules were used to investigate relationships between TILs, genomic alterations, and clinical outcomes in the TCGA-LIHC dataset.

Statistical analysis

Clinicopathological variables associated with MELK expression were analyzed using Pearson chi-squared tests and Fisher’s exact test, as appropriate. Disease-free survival (DFS) was defined as the length of time between surgery and disease recurrence, whereas OS was defined as the length of time from diagnosis to death or the most recent follow-up. Patients for whom these data were not available were excluded from the analysis. Survival outcomes were compared via Kaplan-Meier curves with the log-rank test. Hazard ratios (HRs) with 95% confidence intervals (CIs) were calculated for DFS and OS using univariate Cox proportional hazards regression analysis, with variables found to be significant in the univariate analysis (p < 0.05) being incorporated into a multivariate analysis. A two-sided p < 0.05 was used as the significance threshold. Analyses were conducted in R (v 3.6) and GraphPad Prism 8.3.

MELK expression in pan-cancer

We initially examined MELK expression in 18 different cancer types, finding that MELK levels were significantly elevated in all 18 cancers in comparison with normal tissue. The cancers investigated included bladder carcinoma (BLCA), breast invasive carcinoma (BRCA), cholangiocarcinoma (CHOL), colon adenocarcinoma (COAD), esophageal carcinoma (ESCA), glioblastoma multiforme (GBM), head-neck squamous cell carcinoma (HNSC), kidney chromophobe (KICH) and kidney renal clear cell carcinoma (KIRC), kidney renal papillary cell carcinoma (KIRP), liver hepatocellular carcinoma (LIHC), lung adenocarcinoma (LUAD), lung squamous cell carcinoma (LUSC), prostate adenocarcinoma (PRAD), rectum adenocarcinoma (READ), stomach adenocarcinoma (STAD), thyroid cancer (THCA), and uterine carcinosarcoma (UCEC). The findings were verified in the GEPIA database, which confirmed the significantly elevated MELK levels (Figure 1). Taken together, MELK was upregulated in all the above cancers, suggesting that it may act as a carcinogenic modulator.

MELK and pan-cancer prognosis

The relationship between MELK and survival in the 18 cancer types was evaluated using the GEPIA database, using both overall survival (OS) and disease-free survival (DFS). Elevated levels of MELK in KIRC, KIRP, LIHC, LUAD, and PAAD were linked to poor OS (Figure 2A–2E) while for DFS, the raised expression of MELK in KIRP, LIHC, PAAD, PRAD, STAD, and THCA was found to be associated with poor prognosis (Figure 2F–2K). There were no significant associations between MELK levels and prognosis in the other cancer types. The results suggest the potential of MELK as a biomarker for unfavorable outcomes in HCC.

Overexpression of MELK in HCC

Having established the abnormal overexpression of MELK in HCC, we then investigated its clinical significance. Kaplan-Meier survival curves (Figure 2D and 2G) showed a significant correlation between MELK overexpression and reduced OS. Data from the Human Protein Atlas (HPA) database indicated overexpression of MELK in HCC tissues compared with normal liver tissue (Figure 3A). In addition, we examined MELK genomic expression and copy number to examine the possible reason for its overexpression. Using cBioPortal, we observed MELK amplification in approximately 1.4% of all HCC samples (Figure 3B), while there was no relationship between the copy number and mRNA levels (Figure 3C and 3D). These results indicate that copy number amplification is not the major mechanism responsible for MELK overexpression in HCC.

Overexpression of MELK indicated poor prognosis in HCC

We next examined the relationships between MELK and various clinical parameters in the TCGA-LIHC cohort. The basic characteristics of the patients are shown in Table 1. This showed that MELK levels were significantly linked to both tumor diameter (p = 0.029) and tumor-node-metastasis (TNM) stage (p = 0.032). However, no relationship was observed between MELK and age, sex, BMI, and metastases (both lymph node and distant) (p > 0.05; Table 1).

The median OS for patients with high MELK expression (48.95 ± 8.56 months, 95% CI 32.17–65.73) was significantly shorter than those with low expression (80.68 ± 11.86, 95% CI (57.43–103.93), log-rank p < 0.001). Of the 291 cases, 158 patients (54.3%) experienced tumor recurrence during the follow-up period (median time 21 months.55 ± 2.68, 95% CI 16.29–26.81). Recurrence was greater in the group with high MELK expression (106/187) than in the low-expression group (77/187). The median DFS in the high-expression group (13 months 07 ± 2.22, 95% CI 8.71–17.43) was significantly shorter than in the low-expression group (33 months 90 ± 5.13, 95% CI (23.83–43.97), log-rank p < 0.001).

miRNAs and MELK expression

It is well documented that ncRNAs modulate gene expression. To evaluate the possibility that ncRNAs regulate MELK expression, we predicted ncRNAs that could bind to MELK. This yielded 35 miRNAs; Figure 5A shows a visualization of their interactions with MELK determined by Cytoscape (Figure 5A). According to the established mechanism of miRNA-gene interactions, a negative relationship should exist between the miRNA and the MELK level. Investigation of this issue showed a significant negative association between MELK and has-miR-101-3p (Figure 5B and 5C) with a correlation coefficient R = −0.28 (p = 6.3e-08), but not between MELK and the remaining 34 miRNAs. We then examined the levels and prognostic ability of hsa-miR-101-3p, observing that miR-101-3p expression was reduced in HCC and that higher levels were associated with better prognosis (Figure 5D). This indicates that miR-101-3p negatively regulates MELK in HCC.

lncRNAs upstream of has-miR-101-3p

We next used StarBase to identify lncRNAs upstream of has-miR-101-3p, and a regulatory network of the 25 predicted lncRNAs and has-miR-101-3p was created using Cytoscape (Figure 6A). Measurement of the levels of these lncRNAs in HCC tissue showed significant negative regulation of only LINC00265 and MIR3142HG (Figure 6B and 6C) while LINC00265 and MIR3142HG positively regulated MELK (Figure 6D and 6E). Raised expression of MIR3142HG and LINC00265 was also found to be linked with reduced OS (Figure 6F and 6G). The competing endogenous RNA (ceRNA) hypothesis proposes that lncRNAs modulate mRNA levels by competitive binding to shared miRNAs. This suggests either a negative relationship between lncRNAs and miRNAs or a positive relationship between lncRNAs and mRNAs. These findings indicate that MIR3142HG and LINC00265 are likely candidates for lncRNAs operating upstream of the hsa-miR-101-3p/MELK axis, as shown in Figure 6H.

Association between DNA methylation and MELK expression

To further investigate the mechanisms controlling MELK expression in HCC, we assessed the relationships between MELK levels and methylation. Firstly, we began by comparing the expression of the DNMT1, DNMT3A, and DNMT3B methyltransferases in relation to MELK tumor expression, revealing that all three were up-regulated in the context of higher MELK expression (Figure 7A). Secondly, a UALCAN analysis further revealed a possible link between DNMT1 and increased methylation in normal liver tissue relative to HCC tissue (p = 0.101, Figure 7B), while DiseaseMeth v 2.0 analysis showed that MELK methylation was significantly reduced in HCC tumors compared with paired normal tissues (p < 0.001; Figure 7C). Hypermethylated regions were present within the 3′- and 5′-UTR regions, whereas the TSS1500 and TSS200 regions tended to be hypomethylated. Thirdly, key methylation sites (eg.cg14339556) in the MELK DNA sequences were negatively associated with MELK expression and poor OS (HR = 2.198, p = 0.023, Figure 7D). The differential patterns of MELK expression are represented using heatmaps (Figure 7E).

Immune cell infiltration and MELK

Tumor-infiltrating lymphocytes (TILs) are generally considered to be key independent predictors for lymph node metastasis as well as survival outcomes in many cancer types [17]. Accordingly, we utilized the TIMER database to determine potential links between MELK levels and TILs in HCC. An initial ‘SCNA’ module analysis revealed several infiltrating immune cell populations that were not associated with changes in MELK gene copy number in HCC, including CD4+ T cells, dendritic cells (DCs), B cells, and macrophages (Figure 8A). A subsequent ‘Gene’ module analysis revealed no correlations between MELK and tumor purity, although it was positively linked to B cell, CD8+T, CD4+ T cell, neutrophil, macrophage, and dendritic cell infiltration in HCC (Figure 8B). When we examined the association between MELK levels and patient OS at 60 months, we found that MELK and increased neutrophil (HR = 4.636, 95% CI 1.170, 18.369, p = 0.029) and macrophage (HR = 6.202, 95% CI 1.578, 24.378, p = 0.009) levels were linked with poorer HCC patient survival outcomes, as shown in Figure 8C and 8D. These findings suggest that MELK may impact HCC patient prognosis and clinical outcomes in part by modulating intra-tumoral immune cell infiltration.

To confirm the relationship between MELK expression and TILs, we examined the levels of immunological marker genes associated with six cell types. This revealed that five of these marker genes, namely, CD19, IRF5, ITGAM, NRP1, and ITGAX which are associated with B cells, M1 macrophages, neutrophils, and dendritic cells, respectively, correlated with MELK levels (Table 4). As such, the interplay between MELK and these immune cell populations may shape HCC patient prognosis.

MELK and HCC immune checkpoints

The potential links between MELK and the immune checkpoints PD1/PD-L1 (CD274), CTLA-4, HAVCR2, PDCD1, PDCD1LG2, TIGIT (CADM4), LAG3, and SIGLEC15 were assessed using TIMER. After adjustment for tumor purity, a significant positive correlation was observed between MELK and the expression of these immune checkpoint genes (Figure 9A–9I). This indicated a positive relationship between MELK levels and all the immune checkpoints, implicating MELK in HCC immune escape.