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

Object: Maternal embryonic leucine zipper kinase (MELK) is involved in the development and progression of various cancers. This work investigated the usefulness of MELK in the prediction of hepatocellular carcinoma (HCC) prognosis. Methods: Information on MELK expression was obtained by pan-cancer analysis using The Cancer Genome Atlas (TCGA) database. The TCGA-liver hepatic cancer (TCGA-LIHC), Oncomine datasets, International Cancer Genome Consortium (ICGC) datasets were used to investigate MELK expression in HCC. The prognostic roles of MELK in HCC were assessed by univariate and multivariate survival analyses. The underlying mechanism for noncoding RNAs (ncRNAs) involved in MELK expression was investigated by in silico studies, correlation, methylation, and survival analyses. The relationships between MELK expression and immune cells, immune markers, and checkpoint markers were also analyzed. Results: (1) MELK was identified as an independent predictor of overall survival (OS) in HCC patients (MELK high vs. low expression, HR 2.469; 95% CI 1.217–5.008; p = 0.012) in a multivariate Cox analysis, with a concordance index (C-index) value of 0.727 (95% CI 0.750–0.704). (2) The noncoding RNA miR3142HG and the LINC00265/has-miR-101-3p axis were found to regulate MELK expression in HCC tissue. (3) MELK levels were linked to various immune functions, including tumor infiltration and the expression of immune checkpoints and biomarkers in HCC. Conclusion: MELK may have an oncogenic function in HCC and was found to be up-regulated by ncRNAs and associated with immune cell infiltration and unfavorable prognosis.


INTRODUCTION
Hepatocellular carcinoma (HCC) is a common malignancy, particularly in East Asia, and is associated with high mortality [1,2]. The cancer has a high rate of recurrence, estimated to be between 77% and 100%, largely from the remnant liver, resulting in poor outcomes [3]. The incidence of HCC appears to be increasing, with risk factors including obesity, alcohol abuse, diabetes, chronic viral infection (hepatitis B and C viruses), and metabolic disease [4]. Thus, the determination of biomarkers and targets that will facilitate treatment is necessary.
Maternal embryonic leucine zipper kinase (MELK) is an AMPK serine/threonine kinase [5] and has been proposed as a potential therapeutic target in several cancers [6]. Tang et al. showed that MELK was necessary for proliferation, metastasis, and apoptosis in lung cancer [7], and its use as a diagnostic marker in colorectal [8], ovarian [9], and head and neck cancer has also been reported [10]. MELK promotes mTOR signaling in endometrial cancer [11], and its abnormal expression has been linked to unfavorable breast cancer prognosis [12]. MELK has been found to be targeted by the microRNA (miRNA) miRNA-214-3p in HCC to block proliferation [13]. However, there is as yet no detailed investigation of the role of MELK in HCC, nor is its relationship with immune infiltration of tumors known.
Here, we evaluated the relationship between MELK levels and patient outcomes in a variety of cancers. We also examined MELK regulation in HCC by noncoding RNAs (ncRNAs), as well as an investigation of its function in terms of methylation, tumor infiltration, biomarker expression, and immune checkpoint analysis. The results indicated that up-regulation of MELK by ncRNAs is linked to both tumor infiltration and poor outcome in HCC.

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 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.

Further evaluation of MELK expression in HCC using the Oncomine and ICGC datasets
To further verify the expression pattern and prognostic value of MELK in HCC, we performed meta-analysis using the Oncomine database and survival analysis ICGC database. Oncomine analysis of cancer vs. normal tissue showed that MELK was significantly overexpressed in HCC tissue in different datasets ( Figure 4A). In the Wurmbach Liver dataset, higher MELK mRNA levels were associated with both tumor grade and vascular invasion ( Figure 4B and 4C). In the ICGC liver cancer dataset, higher expression of MELK was associated with poor OS ( Figure 4D).

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 miRNAgene interactions, a negative relationship should exist between the miRNA and the MELK level. Investigation of this issue showed a significant negative association AGING 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.

DISCUSSION
HCC is notable for its high mortality. It is hoped that clarification of its etiology may assist in the development of treatment targets as well as identifying biomarkers to assist diagnosis and prognosis. The involvement of MELK in many human cancers, including HCC, has been documented [18]. However, our understanding of its role in HCC remains largely undetermined. Here, we first analyzed MELK levels in a variety of cancers using data from the TCGA, followed by verification using GEPIA. Further evaluation demonstrated a link between raised MELK levels and poor outcomes in HCC. We also found that the noncoding RNAs LINC00265/has-miR-101-3p mediated upregulation of MELK.
Xia et al. proposed that MELK promotes endometrial carcinoma progression through the E2F1/MELK/ mTORC1/2 axis [19]. Overexpression of MELK was found to correlate with early tumor recurrence and poor patient survival in HCC. The overexpression of MELK in HCC samples strongly correlated with the expression of cell cycle-and mitosis-related genes, while silencing MELK inhibited the cell growth, invasion, stemness, and tumorigenicity of HCC cells by inducing apoptosis and mitosis [19]. These findings, together with our own, suggest that MELK functions as an oncogene in HCC.
Furthermore, the crosstalk between ncRNAs in modulating gene expression through the ceRNA process has been well documented [20,21].
We used seven algorithms, namely, PITA, RNA22, miRmap, microT, miRanda, PicTar, and TargetScan, to identify potential MELK-binding miRNAs. The miRNA has-miR-101-3p was found to be the most likely candidate and may thus be a likely HCC biomarker. miR-101-3p target molecules have been implicated in HCC carcinogenesis [22]. Further investigation indicated that miR-101-3p has tumor-suppressing functions upstream of MELK. This miRNA has been previously found to block HCC proliferation and migration. The ceRNA hypothesis [23] suggests that lncRNAs interacting with the miR-101-3p/MELK axis may have oncogenic functions in HCC. We next identified 25 lncRNAs likely to function upstream of the miR-101-3p/MELK axis, with further examination revealing MIR3142HG and LINC00265 to be the most probable. LINC00265 has been linked to the progression of various cancers, including HCC [24,25]. Taken together, MIR3142HG and the LINC00265/has-miR-101-3p/MELK axis were determined to be potential modulatory factors in HCC. The involvement of this axis explains both MELK's elevated expression and links to poor survival outcomes in HCC, together with the reduction in miR-101-3p.   We did not observe any relationship between MELK levels and copy number variations in HCC. The role of DNA methylation in gene expression is well documented [26,27] and it is possible that this accounts for the dysfunctional expression of MELK in HCC. Hypomethylation of MELK was seen in HCC tissue, which supports the observation of elevated MELK levels together with elevated levels of the DNA methyltransferases (DNMT1, DNMT3A, and DMNT3B). Certain methylation sites were also observed to correlate with HCC patient outcomes.
Immune cell infiltration is known to influence the effectiveness of cancer treatments, including immunotherapy and adjuvant therapy, as well as patient outcomes [28,29]. Here, we found a positive  relationship between MELK and TILs, including B cells, CD8+ T cells, CD4+ T cells, macrophages, neutrophils, and dendritic cells, in HCC. In addition, there was an association between MELK and biomarkers of these cells. MELK was associated with poorer HCC patient survival and neutrophil and macrophage levels. These observations suggest that immune infiltration may be at least partially responsible for MELK's role in HCC carcinogenesis.
Successful immunotherapy relies not only on adequate numbers of TILs but also on the level of immune checkpoint expression [30,31]. For this reason, we examined the link between MELK and immune checkpoints, observing that raised levels of MELK were linked to all the examined checkpoints, implying that targeting MELK could enhance the effectiveness of immunotherapy.
Here, we identified a ceRNA-based MIR3142HG and LINC00265/has-miR-101-3p/MELK axis that may be a potential prognostic biomarker in clinical applications. Nevertheless, several limitations must also be noted. First, the binding affinities of the lncRNAs, miRNAs, and mRNAs obtained from the database require further experimental investigation. Second, further experimental verification of the function and mechanism of the MIR3142HG and LINC00265/has-miR-101-3p/MELK axis in HCC is required.
In conclusion, we demonstrated elevated expression of MELK in various human cancers, including HCC, which was linked to poor survival outcomes. We further showed that MIR3142HG and the LINC00265/ has-miR-101-3p axis acted as upstream regulators of MELK. MELK may thus be a significant prognostic indicator of HCC outcome. Nevertheless, these findings require verification by further research and future clinical trials.