The transcriptome data of The Cancer Genome Atlas Liver Hepatocellular Carcinoma program ‘(TCGA)-LIHC’ and the corresponding clinical data updated, was downloaded on November 3, 2022 from TCGA (https://portal.gdc.cancer.gov/repository) (10), which contains expression profiles of mRNAs and lncRNAs of patients with HCC. After removing duplicate samples, the expression data of 369 patients with HCC and 50 non-tumor liver patients were obtained and converted into TPM format after preprocessing. The GSE144269 dataset was downloaded from Gene Expression Omnibus (GEO) (https://www.ncbi.nlm.nih.gov/geo/, accessed on November 5, 2022) (11) for validation of gene expression and prognosis. Mutation data of HCC were collected and visualized in R using ‘maftools’, and the characteristics of the mutation status were investigated.

The expression patterns of MK5-AS1 and MK5 in HCC tissues and normal tissues from TCGA were compared. GSE144269 was used as an external validation dataset to verify the expression status of the two genes. The 369 patients with HCC were divided into high expression group and low expression group using the medians of MK5-AS1 and MK5, and the relationships between gene expression and clinical parameters were analyzed.

The expression profiles were compared between the MK5-AS1 high and low expression groups; a total of two sets of differentially expressed genes (DEGs) were obtained from TCGA-LIHC and GSE144269. lncRNAs and mRNAs in both gene sets were screened and imported into Metascape (http://metascape.org, accessed on 19 December 2022) (12) for functional enrichment analysis. KEGG enrichment analysis was performed using the ‘clusterProfiler’ R package to screen out relevant pathways (13).

GSEA (14) was performed using the ‘clusterProfiler’ R package in the Molecular Signatures Database (MSigDB)_v7.0_GMTs with the reference dataset ‘c5.all.v7.0.entrez.gmt’. Significantly altered pathways were validated with 1,000 iterative calculations, and the expression level of MK5-AS1 was considered as a phenotypic marker. According to the reference information of GSEA official software, it is generally considered that the NES absolute value is greater than 1, NOM P<0.05, false discovery rate (FDR) q-value <0.25 of the enrichment pathways are significantly enriched between high and low MK5-AS1 expression groups in HCC (15–18).

HCC tissue microarrays used for analysis were commercial products purchased from Shanghai Xinchao Biological Technology Co. Ltd. (cat. no. hlivh180su15). All experiments were performed following the manufacturer's instructions. The corresponding clinical information was also provided by the company. The protocol was approved (approval no. SHYJS-CP-1901001 on 11th January 2019, and extended as approval no. SHYJS-BC-2310001 on 20th October, 2023) by the ethics committee of Shanghai Xinchao Biological Technology Co., Ltd. (Shanghai, China). The tissue microarray contained tumor samples from 90 patients with HCC who received surgical treatment from June 2007 to October 2008. Clinical data included HBsAg, anti-hepatitis C (HCV), alanine aminotransferase (ALT), AFP, and other liver indicators, as well as pathological characteristics such as TNM stage and histological grade. The follow-up period was 3 to 5 years, and data on the overall survival (OS) and disease-free survival (DFS) were collected.

Tissue sections were deparaffinized in xylene for 0.5 h twice. The sections were hydrated in ethanol, followed by high-pressure antigen recovery by heating the tissue with EnVision FLEX TARGET RETRIEVAL SOLUTION LOW pH (pH, 6.1, 3 min; cat. no. K8005; Dako; Agilent Technologies, Inc.). Endogenous peroxidase was blocked by incubating the sections with 3% hydrogen peroxide at room temperature for 15 min. In addition, the sections were blocked with 1X Antibody Diluent/Block (cat. no. ARD1001EA; Akoya Biosciences) at room temperature for 30 min. Subsequently, the slides were stained with a primary antibody against MAPKAPK5 (1:100; cat. no. HPA015515; Atlas Antibodies) at 4°C overnight. The slides were then incubated with an appropriate HRP-conjugated secondary antibody (EnVision FLEX⁺, Mouse, High pH; no dilution required; cat. no. SM802; Dako; Agilent Technologies, Inc.) at 37°C for 30 min. Furthermore, DAB (25°C, 5 min) and hematoxylin (25°C, 1 min) were used for visual antibody staining. An optical microscope was used for observation and to capture images. MAPKAPK5 cytoplasmic staining was scored using four grades (0: negative, 1: weakly positive, 2: moderately positive, 3: strongly positive). The percentage of positive cells was categorized into five grades (0: 0%, 1: 1–5%, 2: 6–25%, 3: 26–50%, 4: 51–100%). The final IHC score was calculated as follows: intensity grad × positive cell percentage grade. A score of 0–3 indicated low expression, and a score of 4–9 indicated high expression.

CIBERSORT is an accurate and robust algorithm that calculates the immune cell composition of tumor tissues by the gene expression profiles (13). The normalized mRNA expression matrix of patients with HCC in TCGA was analyzed using the ‘CIBERSORT’ source R package, which was obtained from the website https://rdrr.io/github/singha53/amritr/src/R/supportFunc_cibersort.R. Other required data were obtained from the supplementary data of a previously published article (13) (https://www.nature.com/articles/nmeth.3337#MOESM207), which included the gene expression matrix for 22 types of immune cells. The number of permutations was 1,000, and the components of various immune cells were compared in MK5-AS1 high and low expression groups. The single sample gene set enrichment analysis (ssGSEA) algorithm was used to assess the infiltration degree of 28 different immune cells. A matrix of immune cell (B cells, CD4⁺ T cells, CD8⁺ T cells, neutrophils, macrophages and dendritic cells) infiltration levels of TCGA-LIHC samples were also downloaded using the TIMER database (https://cistrome.shinyapps.io/timer/ (19), accessed on 21 December 2022) and the correlation between the expression status of MK5-AS1 and immune cell infiltration level in HCC was detected.

A number of clinical studies have demonstrated an association between TMB, MSI and the effect of immunotherapy (20–22). Therefore, such emerging biomarkers may be implicated in the regulation of TME, therefore the relationship between these biomarkers and the expression of MK5-AS1 were explored. Spearman correlation analysis was conducted using SangerBox (http://vip.sangerbox.com/home.html, accessed on 18 October 2022). The TMB value of patients with HCC and high and low expression of MK5-AS1 was calculated. A boxplot of TMB value was generated to visualize the results using the ‘ggbetweenstats’ function of ‘ggstatsplot’ package. Immunotherapy is a strategy for tumor treatment and includes immunotherapy targeting programmed death protein-1 (PD-1)/programmed death protein ligand-1 (PD-L1) and cytotoxic T lymphocyte-associated antigen (CTLA-4). These immunotherapies have demonstrated great efficacy and application prospect in the treatment of advanced HCC (23). Immune checkpoint molecules play an important role in maintaining immune homeostasis and can be exploited for the immune evasion of tumor cells. The relationship between MK5-AS1 and eight immune checkpoints were examined to identify its potential effect in immunotherapy.

The ‘pRRophetic’ R package, developed by Professors Paul Geeleher, Nancy Cox and R. Stephanie Huang at the University of Minnesota, was designed to predict phenotypes by gene expression data (24). It uses data of the CGP cell line from the Cancer Genome Project to predict clinical outcomes and drug sensitivity of external cell lines (Cancer Cell Line Encyclopedia) (25). The potential chemotherapeutic efficacy of 16 common drugs in HCC was examined on the basis of gene expression profiles in TCGA and the differences in the sensitivity of the drugs between MK5-AS1 high and low expression groups were compared.

Diseasemeth 2.0 (http://bio-bigdata.hrbmu.edu.cn/diseasemeth/, accessed on 18 November 2022) (26) was used to examine the methylation pattern of the promoter region of MK5-AS1 and its relation to different clinical stages and histological grades. The somatic mutation data of HCC was downloaded from TCGA and the overall mutation status was analyzed. The results were visualized using the cBioPortal database (http://www.cbioportal.org (27), accessed on 1 February 2023).

Based on high-throughput sequencing data from 15 cell lines, lncATLAS (http://lncatlas.crg.eu/, accessed on 5 September 2022) (28) was used to collect data of 6,768 lncRNAs and to evaluate specific subcellular localization using the ‘relative concentration index’. This database was used to predict the localization of MK5-AS1. Target miRNAs that interact with MK5-AS1 were predicted using the ENCORI platform (http://starbase.sysu.edu.cn/, 5 January 2023) (29), which contains networks of interactions among RNAs. The screening conditions were as follows: i) ‘miRNA-Target: miRNA-lncRNA’; ii) ‘Genome: human’; and iii) ‘Target Gene: lnc-MAPKAPK5-AS1’. The expression of the miRNAs in HCC and the correlation with MK5-AS1 expression using data from TCGA were explored.

WGCNA is an algorithm widely used to identify biomarkers through clustering sets of genes with similar expression patterns (30). It calculates the associations between distinct modules and specific clinical features. A total of 8,900 DEGs were identified from HCC tissues and normal tissues and 5,100 DEGs from the high MK5-AS1 expression group and low MK5-AS1 expression group. A total of 3,255 DEGs were obtained by intersection of the two gene sets. The ‘WGCNA’ R package was applied to construct the co-expression network of 3,255 genes and 369 HCC samples were clustered. β=3 was selected as the soft threshold power to construct the scale-free network. A hierarchical clustering tree was constructed using the dynamic hybrid cutting technology to gather genes with similar expression patterns. The STRING (http://string.embl.de/) (31) database was used to explore the interaction network between proteins, which helps to identify the key regulatory genes.

The cytolytic activity score (CYT) is a robust transcriptome-based immune signature across multiple cancer types and defined as the mean of GZMA and PRF1 expression (TPM format). It has been revealed that a higher CYT is associated with improved outcomes (32).

The IFNG6 score can reflect the overall immune activity and predict the therapeutic effect of pembrolizumab in patients with HCC. This score is calculated from the average expression of six genes (CXCL9, CXCL10, IDO1, IFNG, HLA-DRA and STAT1) (33).

R-4.1.2 and GraphPad Prism-8.00 were used for data cleaning, statistical analysis and graphing. A Kolmogorov Smirnov test was performed to determine the distribution and then based on that, an unpaired t-test was used to explore the difference in MK5-AS1 expression between the HCC and normal groups in the TCGA-LIHC dataset. As for the GSE144269 dataset, the tumor and normal groups were originated from the same patients and then a paired t-test was used to analyze the differences in MK5-AS1 expression between the aforementioned two groups. The histogram was visualized using the median gene expression with an interquartile range. Kruskal-Wallis test with post hoc test (Dunn's or Steel-Dwass) for multiple testing correction was used for multi-group comparison. The chi-square test and Fisher's exact test were employed in the analysis of the relationship between gene expression and clinicopathological characteristics of patients with HCC. Survival curves were drawn using Kaplan-Meier method and comparison among different groups was performed using the log-rank test. Cox regression models were utilized for univariate and multivariate analysis. The prognostic ability of the two genes was assessed by the area under the ROC curve (AUC). Correlations between genes and immune infiltrating cells were compared using Spearman correlation analysis. Two-sided P<0.05 were considered statistically significant in all analyses. Bonferroni correction was used for pairwise comparisons between multiple groups.

As depicted in Fig. 1A and B, the expression levels of MK5-AS1 and MK5 showed significant upregulation in HCC tissues compared with normal liver tissues in TCGA (both P<0.0001). The median with interquartile range of MK5-AS1 in HCC and normal liver tissues was 6.965 and 5.247–10.010 vs. 2.536 and 2.003–3.094, respectively; for MK5, they were 4.268 and 3.199–5.719 vs. 2.180 and 1.870–2.525, respectively. Similar results were observed in the GSE144269 dataset (both P<0.0001, Fig. 1C and D).

Antisense lncRNAs are often correlated with the expression of their sense strand genes, suggesting that they probably be widely involved in the expression regulation of protein-coding genes (5). MK5-AS1 is transcribed from the antisense strand of its protein-coding gene MK5, and the two genes have partially overlapping sequences (Fig. 1E). Spearman's correlation analysis revealed a positive correlation between the two genes (Fig. 1F).

The relationship between the expression levels of MK5-AS1 and MK5 and several widely recognized clinicopathological factors was explored. Analysis of UALCAN (http://ualcan.path.uab.edu/, accessed on 20 October 2022) (34) revealed that a higher expression level of MK5-AS1 was associated with higher clinical stage and histological grade of HCC tissues (Fig. 2A and B). The 369 patients with HCC were divided into the high expression groups and low expression groups using the median of gene expression. As shown in Table I, higher expression of MK5-AS1 was significantly associated with advanced clinical stage (χ²=5.372, P=0.020), T stage (χ²=5.280, P=0.022), histological grade (χ²=17.825, P<0.01) and higher AFP (χ²=29.950, P<0.01). Notably, MK5 exhibited similar tendencies: MK5 expression was negatively linked with the clinical stage (χ²=4.554, P=0.033), T stage (χ²=3.983, P=0.046), histological grade (χ²=17.825, P<0.01) and higher AFP (χ²=23.348, P<0.01) in patients with HCC (Table II). These findings indicated that high expression of MK5-AS1 and MK5 may be potential risk factors for HCC.

Time-dependent ROC curve analysis revealed that MK5-AS1 had high sensitivity and specificity for the five-year survival rate of patients with HCC (AUC=0.78, Fig. 2C). The AUC of MK5 for predicting five-year survival rate was 0.77 (Fig. 2D). To further examine the relationship between gene expression and survival, Kaplan-Meier survival curves were used to examine the effect of MK5-AS1 expression on the overall survival rate of patients with HCC in TCGA and GEO. The median survival time of patients with HCC with high MK5-AS1 expression in TCGA was only 39 months, while it was ~70 months for patients with low MK5-AS1 expression (P=0.01, Fig. 2E). Analysis of GEO data also suggested that the high expression group had improved survival outcomes (P=0.013, Fig. 2F). The results of the survival analysis of MK5 were consistent with these findings; its high expression was an adverse factor for the prognosis of patients with HCC (P=0.013, P=0.05; Fig. 2G and H).

Cox regression analysis can be used to examine whether the gene expression level is a risk factor that affects survival. Univariate Cox regression analysis showed that compared with patients with low MK5-AS1 and MK5 expression, patients with high expression of MK5-AS1 and MK5 indicated a substantially higher risk of mortality. The variables with a statistically significant effect on survival were further included in multivariate Cox regression analysis. The results revealed that MK5-AS1 and MK5 may be independent risk factors for poor survival when M stage was contained (Table III, Fig. 2I). These results suggested that high expression of MK5-AS1 in patients with HCC was associated with tumor progression and adverse prognosis.

Using ‘DESeq2’ R package, 4,179 and 3,706 differentially expressed lncRNAs and mRNAs (|logFC|≥0.6, FDR <0.25, P<0.05) were identified in the TCGA-LIHC and GSE144269 datasets, respectively. A Venn diagram was plotted to select the intersecting genes of the aforementioned two gene sets and 676 genes co-expressed with MK5-AS1 were finally obtained (Fig. 2J). GO analysis demonstrated that MK5-AS1-related genes were primarily enriched in biological regulation, metabolic progress, response to stimulus, multicellular organismal process and immune system process (Fig. 3A). Based on the functional correlation, a network of enriched terms colored by cluster ID was constructed in accordance with correlation and similarity, where nodes that share the same cluster ID were typically close to each other (Fig. 3B). The relative number of genes in each pathway is demonstrated in Fig. 3C; a darker color indicates a greater number of genes, as observed in the pathway ‘biological regulation’. KEGG analysis results suggested that DEGs may be involved in various metabolic-related pathways, such as the PRAK signaling pathway, cell adhesion molecules, carbon metabolism, bile secretion, synaptic vesicle cycle, biosynthesis of amino-acids, complement and coagulation cascades and fatty acid metabolism (Fig. 3D).

GSEA was performed on the basis of normalized enrichment score and the FDR. Regulation of immune system process, biological adhesion, collagen containing extracellular matrix, immune effector process, small molecule metabolic process and defense response were significantly enriched signaling pathways (P<0.05, Fig. 3E).

To make our results more credible at the histological level, an external validation of the gene expression pattern and prognostic significance from 90 HCC tissue chips was carried out. Immunohistochemical semi-quantitative evaluation of clinicopathological specimens showed that the expression score of MK5 in HCC was significantly higher than that in normal tissues (Fig. 4A and B); the median with interquartile range of MK5 in HCC and normal tissues were 1.067 and 1.000–1.500 vs. 0.811 and 0.500–1.000, respectively. Based on a threshold of P<0.05, it was found that the clinical stage (χ²=7.701, P=0.006), ALT (χ²=5.011, P=0.025), and PD-L1 expression (χ²=7.003, P=0.008) was strikingly associated with the expression level of MK5 (Table IV). Nevertheless, the remaining clinicopathological factors such as sex, age, pathology grade, tumor size, recurrence, HBsAg, HBcAb, AntiHCV, AFP and CTLA4 were not statistically significant. As detailed in Fig. 4C and D, patients with HCC with decreased MK5 expression had longer OS and DFS, which shed light on the probability that MK5 acted as a risk factor in the development of HCC.

To clarify the relationship between MK5-AS1 and tumor infiltrating immune cells, the association between MK5-AS1 expression and the infiltration levels of six immune cells in HCC was investigated using the immune cell infiltration data downloaded from the TIMER online database. All correlation analysis conducted with Spearman's test exhibited statistically significant positive correlations, including B cells (R=0.20, P=1.1×10⁻⁴), CD4⁺ T (R=0.31, P=8.9×10⁻¹⁰), CD8⁺ T (R=−0.10, P=0.07), neutrophils (R=0.31, P=2.9×10⁻⁹), macrophages (R=0.34, P=1.5×10⁻⁹) and dendritic cells (R=0.28, P=4.3×10⁻⁷) (Fig. 5A). MK5 was also found to be positively associated with six types of cells of the immune system (Fig. 5B). Furthermore, as suggested in Table V, MK5-AS1 expression is positively linked with multiple immune cell biomarkers in HCC. The aforementioned results demonstrated its effectiveness in regulating TME of HCC.

The immune cell infiltration ratio of TCGA-LIHC was evaluated based on CIBERSORT algorithm. Patients with HCC were divided according to the median of MK5-AS1 expression and the proportions of immune infiltrating cell subtypes in high and low expression groups were calculated (Fig. 5C). Using the R package ‘GSVA’, ssGSEA was used to calculate the abundance of immune cells based on the corresponding data set. For the majority of the 28 types of immune cells, including myeloid-derived suppressor cells, gamma delta T cells, effector memory CD4⁺ T cells, mast cells, memory B cells and natural killer T cells were more significantly enriched in MK5-AS1 high expression group (Fig. 5D). The aforementioned findings indicated that MK5-AS1 may regulate the progression of HCC by affecting cellular infiltration of the immune system.

Inhibiting immune checkpoint signaling pathways are key strategies for the treatment of a range of cancers. The association of MK5-AS1 with ten common immune checkpoints in human cancers, including CD274, CTLA4, HAVCR2, LAG3, PDCD1, PDCD1LG2, TIGIT, SIGLEC15, ITPRIPL1 and IGSF8, was explored. The expression levels of six immune checkpoints were significantly increased in the MK5-AS1 high expression group (Fig. 6A).

As emerging markers of immunotherapy, the predictive value of TMB and MSI in certain cancers has been validated in clinical trials. The effects of TMB and MSI status on the expression level of MK5-AS1 were assessed using the Sangerbox database. Radar charts showed that the expression of MK5-AS1 was positively correlated with TMB and MSI in HCC (R=0.138, P=0.04; R=0.17, P=0.04, respectively) (Fig. 6B and C). The expression of MK5-AS1 revealed a significant effect on the TMB of patients with HCC in the analysis of data from TCGA (Fig. 6E). Collectively, these findings suggested that the group with higher TMB may have a shorter survival time owing to the overexpression of MK5-AS1 in these patients with HCC.

To further explore the clinical significance and drug sensitivity of MK5-AS1, based on ‘pRRophetic’ R package, the potential relationships between MK5-AS1 expression and drug sensitivity of targeted therapeutic drugs that commonly used in patients with HCC were further probed. Obviously, patients with low expression of MK5-AS1 were more sensitive to Axitinib, Bosutinib, Cyclopamine, dasatinib, Docetaxel, Embelin, Gefitinib, Lapatinib, Metformin, Methotrexate and Vorinostat. By contrast, the half-maximal inhibitory concentration (IC₅₀) calculated utilizing ‘pRRophetic’ R package of Bexarotene, Bleomycin, Cisplatin, Doxorubicin and Gemcitabine was lower in MK5-AS1 high expression group (Fig. 6D). These findings suggested that MK5-AS1 might act as an effective biomarker for the efficacy of targeted treatment of patients with HCC.

To seek out the possible reasons for the upregulation of MK5-AS1 in HCC, the methylation level of promoter region near MK5-AS1 was analyzed based on sample type, clinical stage and histological grade and the mutation status of two genes in HCC was probed.

Firstly, the methylation data of the promoter region near MK5-AS1 was obtained using the Diseasemeth 2.0 database. The results revealed that methylation level in HCC tissues was significantly lower than that in normal tissues (Fig. 6F), and there were significant differences in the methylation expression level of MK5-AS1 in different clinical stages and histological grades of HCC (Fig. 6G and H), which provided evidence that the upregulation of MK5-AS1 in HCC might partly due to hypomethylation of its promoter. In addition, the overall result of the ‘MAF’ file was plotted and it was found that missense mutation accounted for the predominant part when the mutation types were classified according to different categories (Fig. 7A). Moreover, single nucleotide polymorphism appeared more frequently than insertions or deletions, with C>T being the most common mutation in single nucleotide variants. Then, the mutation status of MK5-AS1 and MK5 was explored using cBioPortal database, and the results showed that the incidence of MK5-AS1 mutation was 6% (21/360) in HCC, and only two of the 360 patients had missense mutation in MK5, demonstrating that this mutation is fairly unlikely to be the major reason of MK5-AS1 upregulation (Fig. S1A and B).

Evidence suggests that lncRNAs in the cytoplasm can affect the stability and translation regulation of mRNA mainly through the ceRNA regulatory mechanism by adsorbing miRNAs. Using lncATLAS database, the expression pattern of MK5-AS1 in different cell lines was uncovered, while most of the data exhibited that MK5-AS1 is primarily located in the cytoplasm (Fig. 7B). Therefore, it was hypothesized that MK5-AS1 might promote the progression of HCC through the sponge adsorption of miRNAs. ENCORI database unveiled that four miRNAs (hsa-miR-452-5p, hsa-miR-556-3p, hsa-miR-4676-3p and hsa-miR-892c-3p) could directly bind to the gene body of MK5-AS1 and MK5. It has been already noted that hsa-miR-452-5p was overexpressed in HCC tissues compared with normal tissues, which lead to a poor prognosis of HCC through modulating the RNA levels of downstream target genes (35). Unfortunately, based on data from the LIRI-JP dataset of the ICGC database (https://dcc.icgc.org/, accessed on 10 January 2021) (36), there was no discernible variation in the expression pattern of these miRNAs between HCC tissues and normal tissues (Fig. S2A-D). Furthermore, there was no concomitant negative association between these miRNAs and the expression level of MK5-AS1 or MK5 owing to the limitation of sample size and data source (Fig. S2E-L); thus, further molecular experiments are required to confirm the specific mechanism.

A total of eight modules were generated in the hierarchical clustering tree (Fig. 7C). The correlations between all feature genes of these modules and CYT and IFNG6 scores reflecting immune activity are shown in Fig. 7D, in which the blue module showed the strongest association with the aforementioned two scores (R=0.75, P<0.001; R=0.61, P<0.001, Fig. 7F and G). Using MS>0.8 and GS>0.3, 15 genes (NCAPD2, TACC3, PRR11, TPX2, MFSD10, DBN1, ECT2, CENPF, NCAPH, CENPO, TICRR, CHEK1, SPINDOC, H2AX and TRIM59) were screened from the blue module as candidate hub genes. To further explore the biological function of these genes, a protein-protein interaction (PPI) network was constructed using the STRING database (Fig. 7E). Spearman's correlation analysis revealed strong correlation between genes and six types of cells of the immune system (Fig. 7H and I).