Online ISSN: 1945-4589
Learn about our FREE Post-Publication Promotion Services
Research Paper Volume 15, Issue 21 pp 11918—11939

The development and experimental validation of hypoxia-related long noncoding RNAs prognostic signature in predicting prognosis and immunotherapy of cutaneous melanoma

Gang Wang1, *, , Yuliang Sun1, *, , Qingjia Xu1, ,

  • 1 Department of Orthopedics, Qilu Hospital of Shandong University, Jinan, China
* Equal contribution and share first authorship

Received: May 23, 2023       Accepted: September 26, 2023       Published: November 2, 2023      

https://doi.org/10.18632/aging.205157
How to Cite

Copyright: © 2023 Wang et al.. This is an open access article distributed under the terms of the Creative Commons Attribution License (CC BY 4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Abstract

Cutaneous melanoma (CM) is widely acknowledged as a highly aggressive form of malignancy that is associated with a considerable degree of morbidity and poor prognosis. Despite this recognition, the precise role of hypoxia-related long noncoding RNAs (HRLs) in the pathogenesis of CM remains an area of active research. This study sought to elucidate the contribution of HRLs in CM by conducting a thorough screening and extraction of hypoxia-related genes (HRGs). In particular, we conducted univariate and multivariate Cox regression analyses to assess the independence of the prognostic signature of HRLs. Our results demonstrated that a novel risk model could be established based on five prognostic HRLs. Remarkably, patients with low-risk scores exhibited significantly higher overall survival rates compared to their high-risk counterparts, as confirmed by Kaplan-Meier survival analysis. Furthermore, we utilized consensus clustering analysis to categorize CM patients into two distinct subtypes, which revealed marked differences in their prognosis and immune infiltration landscapes. Our nomogram results confirmed that the HRLs prognostic signature served as an independent prognostic indicator, offering an accurate evaluation of the survival probability of CM patients. Notably, our findings from ESTIMATE and ssGSEA analyses highlighted significant disparities in the immune infiltration landscape between low- and high-risk groups of CM patients. Additionally, IPS and TIDE results suggested that CM patients in different risk subtypes may exhibit favorable responses to immunotherapy. Enrichment analysis and GSVA results indicated that immune-related signaling pathways may mediate the role of HRLs in CM. Finally, our tumor mutation burden (TMB) results indicated that patients with low-risk scores had a higher TMB status. In summary, the establishment of a risk model based on HRLs in this study provided an accurate prognostic prediction and correlated with the immune infiltration landscape of CM, thereby providing novel insights for the future clinical management of this disease.

Introduction

Cutaneous melanoma (CM) is a highly metastatic and aggressive malignancy with an increasing incidence worldwide [1]. Accounting for 5% of skin malignancies and causing more than 70% of deaths, CM poses a significant challenge to human health [2]. The main cause of death in melanoma patients is widespread, and the 5-year survival rate is less than 23% [3]. Despite the promising attention received by immunotherapies such as PD-1/PD-L1 and CTLA4 in CM treatment, the prognosis of CM remains unsatisfactory [4]. Therefore, there is an urgent need to investigate the mechanisms underlying the tumorigenesis and progression of CM to identify novel markers for its diagnosis and treatment.

Hypoxia is the most important and prevalent characteristic of the microenvironment in tumors and is closely associated with tumor proliferation and metastasis, which typically indicates poor prognosis [5]. Several reports have indicated that the majority of malignant tumors are associated with hypoxia, including prostate cancer, glioblastoma multiforme, malignant melanoma, breast cancer, and metastatic liver cancer [6]. In tumor biology, the vasculature of tumors is often unable to keep pace with the rapidly proliferating tumor cells, resulting in highly hypoxic regions. This aggressive phenotype in tumors is further conferred by upregulating angiogenic, survival, proliferative, and metastatic pathways [7]. Previous studies have demonstrated that hypoxia is a well-accepted aggravating factor in tumor development that facilitates metastasis, such as promoting lymph node metastasis in melanoma [8]. It has been reported that hypoxia can promote uveal melanoma cell angiogenesis and metastasis by upregulating the expression of glycosylate - secreted protein ANGPTL4 and VEGF [9]. Notably, as a characteristic feature of virtually all solid tumors, hypoxia can also directly modulate the tumor immune microenvironment [10]. Hypoxia can lead to tumor immunosuppression and immune escape [11]. However, few studies have described the underlying mechanisms of hypoxia in CM.

Long non-coding RNAs (lncRNAs) are a novel class of non-protein-coding transcripts longer than 200 nucleotides that lack apparent protein coding potential [12, 13]. Previous studies have reported that lncRNAs are expressed abnormally in multiple diseases and participate in various biological and physiological processes, including cell proliferation, apoptosis, and migration [14, 15]. With the in-depth study of lncRNAs, dysregulation of lncRNAs has been found to be involved in the tumorigenesis and development of human tumors, including melanoma [16]. Furthermore, lncRNAs have close relationships with the development of CM, including cell cycle arrest, inhibition of tumor microenvironment formation, activation of tumor cell signal pathway, and poor prognosis [17, 18]. Nevertheless, the role of hypoxia-related long noncoding RNAs (HRLs) in CM remains elusive.

Currently, bioinformatics approaches are widely utilized to characterize diseases, and lncRNA-based prognostic models have been developed to evaluate patient prognosis [15]. In this study, an HRL-based model was established using The Cancer Genome Atlas (TCGA) database, and 5 HRL signatures were used to predict the prognosis and evaluate the immune infiltration landscape of CM patients. Moreover, the response to drug sensitivity and immunotherapy of patients in different risk subtypes was comprehensively investigated, and the possible molecular mechanisms were illustrated in detail. Abnormal expressions of HRLs in CM were further verified in cell lines using qRT-PCR. Collectively, the findings of this study provide novel insights and perspectives for the clinical management of CM patients.

Materials and Methods

TCGA dataset collection

For this study, we obtained the transcriptome matrix in RNA-seq FPKM format and clinical information from The Cancer Genome Atlas database (TCGA) (https://portal.gdc.cancer.gov/). We excluded patients with melanoma who had no survival time or a survival time of less than zero from the analysis, resulting in 454 melanoma samples for further analysis. To annotate the transcriptome matrix symbols (mRNA and lncRNA), we used the ensembles human genome browser GRCh38.p13 (http://asia.ensembl.org/index.html) with the assistance of Perl scripts. We retrieved the patients’ clinical characteristics, such as age, gender, stage, and TN stage, from the TCGA database using Perl scripts. We excluded samples in M stage due to the significant difference in sample size. All clinical information and data were obtained from public databases, and therefore, written informed consent from patients and approval from the ethics committee were not required for this study.

Identification of HRLs and risk model construction

We obtained 200 hypoxia-related genes (HRGs) from the Molecular Signature Database (MSigDB) (https://www.gsea-msigdb.org/gsea/) (details provided in Supplementary Table 1). Subsequently, we performed Pearson correlation analysis to identify a set of 186 lncRNAs that were significantly associated with HRGs, termed HRLs, based on a threshold of |correlation coefficient| > 0.5 and P < 0.001 (|r| > 0.5, P < 0.001) (details provided in Supplementary Table 2). Using Perl scripts, we extracted the expression levels of the HRLs from the transcriptome matrix and merged them with the corresponding clinical characteristics information for further analysis. We performed univariate Cox regression analysis using the R package “survival” to identify HRLs significantly associated with the overall survival (OS) rate of CM. From this, we selected 51 HRLs that met the predefined criteria (details provided in Supplementary Table 3) for further analysis. The least absolute shrinkage and selection operator (LASSO) regression analysis was employed using the R package “glmnet” to identify the most characteristic variables among the prognostic HRLs. The candidate HRLs were selected using multivariate Cox regression analysis, which was performed using the R package “survival” to establish the risk model. The risk score for each CM patient was calculated using the following formula: risk score = (-0.357 x expression of LINC00324) + (0.221 x expression of EBLN3P) + (0.218 x expression of MIR205HG) + (-0.118 x expression of THCAT158) + (-0.439 x expression of USP30-AS1). Based on the median risk score, patients were divided into low- and high-risk groups. We used the Kaplan-Meier survival curve to evaluate the OS rate of patients in the two risk groups using the log-rank algorithm with the R package “survival”. We used the R package “ggplot2” to perform principal component analysis (PCA) to investigate the distribution pattern between the low- and high-risk groups. We visualized the expression of HRLs in the low- and high-risk groups using the R package “pheatmap”.

Function enrichment analysis and tumor mutational burden landscape

The tumor mutation data in maf format of CM samples were retrieved from the TCGA database. Perl scripts were used to extract the mutation data from the raw data, and the “Maftools” package in R software was employed to create a waterfall diagram. To identify differentially expressed genes (DEGs) between patients in the low- and high-risk group, the R package “limma” was utilized, with a threshold of |fold change| ≥ 2 and P < 0.05. The KEGG terms of CM patients in the low- and high-risk group were calculated using Gene Set Variation Analysis (GSVA) with P < 0.05 considered significantly different. Furthermore, the “clusterProfiler” R package was used to conduct Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis to enrich the DEGs into pathways [19].

Independent prognosis analysis

Univariate and multivariate Cox regression analyses are commonly used statistical methods to examine the association between survival outcomes and various factors, such as the HRLs prognostic signature and clinicopathological characteristics. The pROC package is utilized to assess the accuracy of diagnostic tests and predict AUC values. The timeROC package is employed to compute the AUC values for different time points. Nomogram models can integrate diverse prognostic factors to predict the survival probability, and the rms package is used to construct nomograms. The C-index plot is a tool to evaluate the predictive accuracy of survival models. The regplot package can be used to generate calibration diagrams that compare the predicted survival probability with the actual survival probability. The ggDCA package is employed to perform decision curve analysis, which assesses the clinical utility of prediction models.

Validation of risk model and consensus clustering analysis

The CM samples were divided into training and test cohorts in a random manner using the R package “caret,” at a ratio of 7:3. A total of 318 CM samples were assigned to the training cohort, while 136 samples were assigned to the test cohort. Utilizing the prognostic HRLs, the risk score for each CM sample was calculated according to the risk formula in both cohorts. The patients in both cohorts were subsequently categorized into low- and high-risk groups based on the median risk score. To classify the CM samples into different molecular subtypes, the R package “ConsensusClusterPlus” was employed. The clustering was based on partitioning around medoids with “euclidean” distances, using 1000 iterations, and with the maximum K value set at 9. Based on the optimal classification of K between 2 and 9, the CM samples were assigned to different molecular subtypes for further analysis.

Cell culture and qRT-PCR analysis

The human fibroblasts cell line HFB4 and human melanoma cell line A375 were obtained from the American Type Culture Collection (ATCC). To initiate the culture, the cryopreserved A375 and HFB4 cells were thawed in a 37° C water bath and transferred to sterile 15 mL centrifuge tubes containing 10 mL of DMEM/F12 culture medium under aseptic conditions. The tubes were then incubated in a humidified incubator at 37° C with 5% CO2 to promote cell growth and maintenance of cell viability. Subsequently, RNA extraction from both cell lines was performed using Trizol reagent (Catalogue number: 15596018, Thermo Fisher Scientific), followed by cDNA synthesis using a reverse transcription kit with gDNA Eraser (Perfect Real Time, Takara Bio). Real-time quantitative qRT-PCR (Catalogue number: RR047A, Takara Bio) was carried out to perform further analysis. Finally, the lncRNA expression levels were measured using SYBR Pre-mix Ex Taq II (TliRNaseH Plus) (Catalogue number: RR820B, Takara Bio).

Immune infiltration landscape, immunotherapy response and drug sensitivity analysis

ESTIMATE algorithm was used to estimate the proportion of stromal and immune cells in CM samples, and the stromal, immune, ESTIMATE scores, and tumor purity were estimated using R package “estimate”. Then, single sample gene set enrichment analysis (ssGSEA) algorithm was used to evaluate the proportion of 23 types of immune cells and the immune function score of each CM sample via the “GSVA” R package. Spearman-ranked correlation analysis was used to investigate the association of prognostic IHRLs and immune cells, which were visualized in a heatmap using the “ggplot2” R package. In addition, the Immunophenoscore (IPS) results of CM patients were downloaded from the TCIA database, and the TIDE score of each sample was calculated using the TIDE database. Finally, the Genomics of Drug Sensitivity in Cancer (GDSC) database was utilized to evaluate drug sensitivity (IC50), and the response to antineoplastic drugs for each CM sample was predicted using the “pRRophetic” R package. Correlation analysis was used to investigate the correlation between risk score and drug sensitivity (IC50), and all statistical analysis were visualized using the “ggplot2” R package.

Statistical analysis

In this study, all statistical analyses were conducted using R software version 4.1.0 (http://www.R-project.org) and Perl scripts. The correlation between two variables was calculated using the Spearman’s rank correlation algorithm, and statistical significance was set at a threshold of P < 0.05. Differential functions were analyzed using the Wilcoxon rank-sum test between the two groups, and statistical significance was set at a threshold of p < 0.05.

Availability of data and materials

The data used to support the findings of this study are included within the article. The data and materials in the current study are available from the corresponding author on reasonable request.

Results

Risk model development based on HRLs prognostic signature in CM

In this study, we developed a novel risk model to investigate the prognostic value of HRLs in predicting the prognosis for CM. To identify the HRLs, we performed Pearson correlation analysis and identified 186 lncRNAs that were associated with HRGs (Supplementary Figure 1A). Then, we used the least absolute shrinkage and selection operator (LASSO) analysis and identified 8 HRLs that were associated with OS rate (Figure 1A and Supplementary Figure 1B) based on univariate Cox regression analysis. Multivariate Cox regression analysis showed that 5 HRLs could independently predict the OS rate of CM, and they were used to construct the risk model. Based on the median of risk score, we ranked and divided the patients with CM into low- and high-risk groups. The scatter dot plot showed that the HRLs prognostic signature was inversely associated with survival time in CM (Figure 1B). The Kaplan-Meier survival curve analysis indicated that patients with low-risk scores had a significantly higher OS rate compared to those with high-risk scores (Figure 1C). Moreover, principal component analysis (PCA) illustrated a clear separation between the low- and high-risk groups based on the HRLs prognostic signature (Figure 1D). The heatmap visualizable diagram result showed that the expression of LINC00324, USP30-AS1, EBLN3P, and THCAT158 were significantly higher in the low-risk group, whereas the expression of MIR205HG was higher in the high-risk group (Figure 1E). These findings indicate that the construction of the risk model for the HRLs prognostic signature is closely associated with the prognosis of patients with CM.

Risk model construction based on the HRLs prognostic signature in CM. (A) LASSO regression analysis shows the optimal coefficients and minimum lambda. (B) Risk score distribution, and the scatter dot plot shows the correlation between risk score and survival time in CM. (C) Kaplan-Meier survival curve analysis suggests that the OS rate of patients with low-risk score is longer than those with high-risk score. (D) Principal component analysis illustrates a significant separation between low- and high-risk group. (E) Heatmap diagram shows the expression of 5 prognostic HRLs of patients in the low- and high-risk group.

Figure 1. Risk model construction based on the HRLs prognostic signature in CM. (A) LASSO regression analysis shows the optimal coefficients and minimum lambda. (B) Risk score distribution, and the scatter dot plot shows the correlation between risk score and survival time in CM. (C) Kaplan-Meier survival curve analysis suggests that the OS rate of patients with low-risk score is longer than those with high-risk score. (D) Principal component analysis illustrates a significant separation between low- and high-risk group. (E) Heatmap diagram shows the expression of 5 prognostic HRLs of patients in the low- and high-risk group.

Validation of risk model based on HRLs prognostic signature

An internal validation was conducted to assess the accuracy and independence of the HRLs prognostic signature in predicting the prognosis of CM patients. The patients were randomly split into a training cohort and a test cohort in a 7:3 ratio, resulting in 318 samples in the training cohort and 136 samples in the test cohort. Using the 5 prognostic HRLs, the risk score of each CM sample was calculated and categorized into low- and high-risk groups in both cohorts. The patients in the training cohort were ranked based on their median risk score, and the scatter dot plot showed an inverse correlation between the risk score and survival time (Figure 2A). Similarly, the patients in the test cohort were ranked according to the median risk score, and the scatter dot plot indicated an inverse correlation between the risk score and survival time (Figure 2C). The Kaplan-Meier survival curve analysis showed that patients in the low-risk group had a significantly higher OS rate than those in the high-risk group (Figure 2D). The PCA score plot revealed a clear separation between the low- and high-risk groups in both cohorts (Figure 2E–2H). Furthermore, the heatmap diagram showed that the low-risk group had higher expression levels of LINC00324, USP30-AS1, EBLN3P, and THCAT158, while the high-risk group had a higher expression level of MIR205HG, in both cohorts. These findings suggest that the HRLs prognostic signature accurately predicts the prognosis of CM patients.

Risk model construction in training cohort and test cohort based on the prognostic HRLs. (A) Distribution of risk score in training cohort. (B) Kaplan-Meier survival curve analysis of patients with CM in training cohort. (C) Distribution of risk score in test cohort. (D) Kaplan-Meier survival curve analysis of patients with CM in test cohort. (E) PCA score plot shows a clear separation between low- and high-risk group in training cohort. (F) Heatmap diagram displays the expression of 5 prognostic HRLs in training cohort. (G) PCA score plot shows a clear separation between low- and high-risk group in test cohort. (H) Heatmap diagram displays the expression of 5 prognostic HRLs in test cohort.

Figure 2. Risk model construction in training cohort and test cohort based on the prognostic HRLs. (A) Distribution of risk score in training cohort. (B) Kaplan-Meier survival curve analysis of patients with CM in training cohort. (C) Distribution of risk score in test cohort. (D) Kaplan-Meier survival curve analysis of patients with CM in test cohort. (E) PCA score plot shows a clear separation between low- and high-risk group in training cohort. (F) Heatmap diagram displays the expression of 5 prognostic HRLs in training cohort. (G) PCA score plot shows a clear separation between low- and high-risk group in test cohort. (H) Heatmap diagram displays the expression of 5 prognostic HRLs in test cohort.

The prognostic signature based on HRLs was an independent prognosis indicator for CM

Univariate and multivariate Cox regression analyses were employed to assess the ability of the risk score based on HRLs prognostic signature to serve as an independent prognostic indicator for CM. The results of the univariate analysis showed a close association between age (HR = 1.020, P < 0.001), stage (HR = 1.473, P < 0.001), T (HR = 1.445, P < 0.001), N (HR = 1.443, P < 0.001), risk score (HR = 1.798, P < 0.001) and OS rate in patients with CM (Figure 3A). On the other hand, the multivariate Cox regression analysis indicated that age (HR = 1.013, P = 0.025), T (HR = 1.315, P = 0.002), N (HR = 1.550, P < 0.001), and risk score (HR = 1.546, P < 0.001) were all independent prognostic indicators for CM patients (Figure 3B). The ROC curve demonstrated that the risk model had an AUC of 0.729, indicating a satisfactory predictive ability for CM (Figure 3C). Furthermore, a stratified subgroup analysis was conducted to evaluate the prognostic value of the HRLs prognostic signature in different clinicopathological characteristics. The patients with CM were categorized into low- and high-risk groups based on the median risk score, and the analysis was performed across various clinicopathological characteristics, including age (age ≤ 65 vs age > 65), gender (female vs male), N (N 0-1- vs N 2-3), stage (stage 0-1 vs stage 2-4), T (T 0-1 vs T 2-4). The Kaplan-Meier survival curve analysis revealed that the OS rate of patients with a low-risk score was significantly higher than those with a high-risk score across different clinicopathological characteristics (Figure 3D–3M). These results collectively demonstrate that the risk score based on the HRLs prognostic signature is an independent prognostic indicator that can effectively predict the prognosis of CM patients in comparison to other clinicopathological characteristics.

Independent prognosis analysis of HRLs prognostic signature and clinicopathological characteristics. (A) Univariate Cox regression analysis shows that age, stage, T, N, and risk score are closely associated with OS rate in CM. (B) Multivariate Cox regression analysis reveals that risk score is an independent prognostic indicator of patients with CM. (C) ROC curve shows the AUC of HRLs prognostic signature and different clinicopathological characteristics. (D–M) The Kaplan-Meier survival curve shows the OS rate of patients with low- and high-risk score in different clinicopathological characteristics.

Figure 3. Independent prognosis analysis of HRLs prognostic signature and clinicopathological characteristics. (A) Univariate Cox regression analysis shows that age, stage, T, N, and risk score are closely associated with OS rate in CM. (B) Multivariate Cox regression analysis reveals that risk score is an independent prognostic indicator of patients with CM. (C) ROC curve shows the AUC of HRLs prognostic signature and different clinicopathological characteristics. (D–M) The Kaplan-Meier survival curve shows the OS rate of patients with low- and high-risk score in different clinicopathological characteristics.

Nomogram construction based on the HRLs prognostic signature and clinicopathological characteristics

A nomogram has been developed to accurately assess the probability of survival at 1, 3, and 5 years for patients with CM, based on their HRLs prognostic signature and clinicopathological characteristics (Figure 4A). The concordance index (C-index) curve showed that the predictive capability of the HRLs signature in predicting prognosis of CM patients was better than other clinicopathological characteristics (Figure 4B). The calibration curve indicated that the predicted OS rates for 1, 3, and 5 years using the nomogram were consistent with the actual OS rates (Figure 4C). The decision curve analysis (DCA) and ROC curve results also suggested a satisfactory predictive ability of the nomogram for predicting the survival probability of CM patients (Figure 4D, 4E). The time-dependent ROC curve indicated that the AUC for 1, 3, and 5 years was 0.689, 0.658, and 0.699, respectively (Figure 4F). In summary, these findings demonstrate that the nomogram construction based on HRLs prognostic signature accurately predicts the prognosis of CM patients and is highly reliable.

Nomogram construction based on HRLs prognostic signature and clinicopathological characteristics. (A) Construction of nomogram model to predict 1-, 3-, and 5-years survival probability of patients with CM. (B) Concordance index curve of risk score and clinicopathological characteristics. (C) Calibration curve reveals the consistence of nomogram-predict OS and actual OS. (D) Decision curve analysis (DCA). (E) ROC curve shows the AUC of nomogram, risk score, and clinicopathological characteristics. (F) Time-dependent ROC curve shows the AUC of 1-, 3-, and 5-years.

Figure 4. Nomogram construction based on HRLs prognostic signature and clinicopathological characteristics. (A) Construction of nomogram model to predict 1-, 3-, and 5-years survival probability of patients with CM. (B) Concordance index curve of risk score and clinicopathological characteristics. (C) Calibration curve reveals the consistence of nomogram-predict OS and actual OS. (D) Decision curve analysis (DCA). (E) ROC curve shows the AUC of nomogram, risk score, and clinicopathological characteristics. (F) Time-dependent ROC curve shows the AUC of 1-, 3-, and 5-years.

Functional enrichment analysis of differential expressed genes (DEGs)

To investigate the potential molecular mechanisms of CM patients in the low- and high-risk groups, GSVA and enrichment analysis were employed. The volcano diagram (Figure 5A) illustrates the differentially expressed genes (DEGs) in the low- and high-risk groups. GSVA was used to calculate the activity of KEGG pathways, and the results suggested a remarkable down-regulation in immune-related signaling pathways of CM patients in the high-risk group (Figure 5B). KEGG enrichment analysis indicated that the DEGs were significantly enriched in cytokine-cytokine receptor interaction, cell adhesion molecules, and chemokine signaling pathways (Figure 5C). GO enrichment analysis illustrated that DEGs were enriched in immune-related biological processes, such as leukocyte-mediated immunity, positive regulation of cell activation, and positive regulation of leukocyte activation (Figure 5D). These results suggest that immune-related processes may play a role in mediating the effects of HRLs in CM patients.

Functional enrichment analysis of differential expression genes (DEGs) in low- and high-risk group. (A) Volcano diagram shows the DEGs in the low- and high-risk group with the threshold set at |Fold Change| ≥ 2 and P B) GSVA reveals the activity of KEGG signal pathways of each CM patient in the low- and high-risk group. (C) KEGG enrichment analysis of DEGs. (D) GO enrichment analysis of DEGs.

Figure 5. Functional enrichment analysis of differential expression genes (DEGs) in low- and high-risk group. (A) Volcano diagram shows the DEGs in the low- and high-risk group with the threshold set at |Fold Change| ≥ 2 and P < 0.05. (B) GSVA reveals the activity of KEGG signal pathways of each CM patient in the low- and high-risk group. (C) KEGG enrichment analysis of DEGs. (D) GO enrichment analysis of DEGs.

Consensus clustering analysis and immune infiltration landscape

The molecular subtypes of cutaneous melanoma (CM) were further investigated using a set of five prognostic HRLs. Consensus clustering was employed to classify the CM patients into different molecular subtypes, resulting in an optimal classification of K=2, with 309 samples in Cluster A and 145 samples in Cluster B, as illustrated in the heatmap (Figure 6A). Kaplan-Meier survival analysis revealed a significant difference in the overall survival (OS) rate between patients in Cluster A and Cluster B, with Cluster A exhibiting a higher OS rate than Cluster B (Figure 6B). Additionally, the results of principal component analysis (PCA) demonstrated a clear separation of patients in Cluster A and Cluster B based on the five prognostic HRLs (Figure 6C). Furthermore, the ESTIMATE assessment algorithm indicated that patients in Cluster B had higher stromal, immune, and ESTIMATE scores and lower tumor purity than those in Cluster A (Figure 6D–6G). Similarly, the results of the ssGSEA algorithm revealed that the proportion of most immune cells was higher in patients in Cluster B (Figure 6H). Moreover, immune function analysis demonstrated that patients in Cluster B had a higher immune function score than those in Cluster A, as evidenced by higher cytolytic activity and T cell co-inhibition (Figure 6I). These findings demonstrate that the five prognostic HRLs can accurately classify CM samples into different molecular subtypes that are associated with prognosis and immune infiltration landscape.

Consensus clustering of CM samples and immune infiltration landscape evaluation. (A) Consensus clustering heatmap shows the optimal classification under K= 2-9. (B) The Kaplan-Meier survival curve analysis of patients with CM in Cluster A and Cluster B. (C) PCA score plot illustrates a clear distribution between Cluster A and Cluster B. (D–G) Stromal, immune and ESTIMATE scores, and tumor purity. (H) The proportion of 23-type immune cells in Cluster A and Cluster B. (I) Immune function score of patients in Cluster A and Cluster B.

Figure 6. Consensus clustering of CM samples and immune infiltration landscape evaluation. (A) Consensus clustering heatmap shows the optimal classification under K= 2-9. (B) The Kaplan-Meier survival curve analysis of patients with CM in Cluster A and Cluster B. (C) PCA score plot illustrates a clear distribution between Cluster A and Cluster B. (D–G) Stromal, immune and ESTIMATE scores, and tumor purity. (H) The proportion of 23-type immune cells in Cluster A and Cluster B. (I) Immune function score of patients in Cluster A and Cluster B.

Association of HRLs prognostic signature and immune infiltration landscape

The present study investigated the association between the HRLs prognostic signature and the immune infiltration landscape. The ESTIMATE analysis showed that patients in the low-risk group had higher stromal, immune, and ESTIMATE scores, and lower tumor purity (Figure 7A–7D). Additionally, the immunotherapy prediction analysis indicated that patients with high-risk scores had lower TIDE scores, suggesting a better response to immunotherapy for patients in the high-risk group (Figure 7E). Furthermore, the response to anti-CTLA4 and anti-PD-1 immunotherapies was evaluated for CM patients in the low- and high-risk groups. The IPS analysis suggested that the low-risk group showed a promising response to anti-CTLA4, antiPD-1, and anti-CTLA4/PD-1 (Figure 7F–7H). ssGSEA and immune function score analysis revealed that patients in the low-risk group had a higher proportion of immune cells and immune function scores than those in the high-risk group (Figure 7I, 7J). Moreover, the expression of ICI analysis showed higher expression of LAG3, CTLA4, PD−1, PDCD1LG2, and PD−L1 in the low-risk group (Figure 7K). Finally, a correlation analysis was conducted to investigate the association between the 5 prognostic HRLs and immune infiltration landscape. The results showed that LINC00324 and USP30−AS1 were positively correlated with 23 types of immune cells, while EBLN3P and THCAT158 were negatively associated with most of the 23 types of immune cells (Figure 7L). These findings suggest that the risk model based on HRLs prognostic signature is closely associated with the immune infiltration landscape and may aid in the evaluation of the immunotherapy response of CM patients in different risk subgroups.

Association of risk score and immune infiltration landscape in CM. (A–D) Stromal, immune, ESTIMATE scores, and tumor purity. (E) TIDE score. (F–H) IPS score. (I) The proportion of 23-type immune cells of patients in the low- and high-risk group. (J) Immune function score. (K) Expression of immune checkpoints inhibitor (ICI) of in low- and high-risk group. The expression of ICI was transformed by log2 (expression + 1). (L) Correlation analysis of 5 prognostic HRLs and 23-type immune cells.

Figure 7. Association of risk score and immune infiltration landscape in CM. (A–D) Stromal, immune, ESTIMATE scores, and tumor purity. (E) TIDE score. (F–H) IPS score. (I) The proportion of 23-type immune cells of patients in the low- and high-risk group. (J) Immune function score. (K) Expression of immune checkpoints inhibitor (ICI) of in low- and high-risk group. The expression of ICI was transformed by log2 (expression + 1). (L) Correlation analysis of 5 prognostic HRLs and 23-type immune cells.

Correlation analysis of risk score and drug sensitivity

Targeted drug therapy has emerged as a promising approach in treating patients with cutaneous melanoma (CM). In light of the substantial differences in immunotherapy responses among CM patients in various risk subgroups, we conducted further investigations to assess the sensitivity of antineoplastic drugs in different risk groups. Our results, as depicted in Figure 8A–8H, demonstrate that the IC 50 values of Paclitaxel, AKT inhibitor VIII, Rapamycin, Pazopanib, Lapatinib, Crizotinib, and Sunitinib were significantly higher in the high-risk group, whereas the IC 50 of Sorafenib was significantly higher in the low-risk group. Furthermore, correlation analysis revealed a positive correlation between risk score and Paclitaxel, AKT inhibitor VIII, Rapamycin, Pazopanib, Lapatinib, Crizotinib, and Sunitinib, but a negative correlation with Sorafenib (Figure 8I–8P). Taken together, these results suggest that CM patients in different risk subgroups may exhibit distinct sensitivities to antineoplastic drugs, providing a new perspective for the development of precisely targeted drug and chemotherapy for CM.

Correlation analysis of risk score and drug sensitivity. IC 50 of (A) Paclitaxel, (B) Sorafenib, (C) AKT inhibitor VIII, (D) Rapamycin, (E) Pazopanib, (F) Lapatinib, (G) Crizotinib, and (H) Sunitinib. (I–P) Correlation analysis of risk score and frug sensitivity.

Figure 8. Correlation analysis of risk score and drug sensitivity. IC 50 of (A) Paclitaxel, (B) Sorafenib, (C) AKT inhibitor VIII, (D) Rapamycin, (E) Pazopanib, (F) Lapatinib, (G) Crizotinib, and (H) Sunitinib. (I–P) Correlation analysis of risk score and frug sensitivity.

The landscape of somatic gene mutations based on the HRLs prognostic signature

Tumor mutational burden (TMB) has emerged as a promising biomarker for predicting immunotherapy response in tumors. Our analysis of TMB indicated that patients in the low-risk group exhibited a higher TMB compared to those in the high-risk group (Figure 9A). Furthermore, Kaplan-Meier survival curve analysis revealed that patients with a high-risk score had a significantly lower overall survival (OS) rate than those with a low-risk score (Figure 9B). Additionally, the OS rate of patients with a low-risk score was significantly higher than those with a high-risk score in both H-TMB and L-TMB groups (Figure 9C). Analysis of the mutation frequencies of genes revealed that the low-risk group had higher mutation frequencies in most genes, including TTN, MUC16, BRAF, DNAH5, and PCLD (Figure 9D, 9E). These findings suggest that TMB may be a useful biomarker for predicting immunotherapy response in CM patients, and that patients in the low-risk group may benefit from immunotherapy treatment.

The tumor mutational burden landscape of patient in the low- and high-risk group. (A) TMB analysis in the low- and high-risk group. (B) The Kaplan-Meier survival curve analysis of patients with low- and high-TMB. (C) Kaplan-Meier survival curve of patient with L-TMB and H-TMB in the low- and high-risk group. (D, E) The mutant landscape of CM patients in the low- and high-risk group.

Figure 9. The tumor mutational burden landscape of patient in the low- and high-risk group. (A) TMB analysis in the low- and high-risk group. (B) The Kaplan-Meier survival curve analysis of patients with low- and high-TMB. (C) Kaplan-Meier survival curve of patient with L-TMB and H-TMB in the low- and high-risk group. (D, E) The mutant landscape of CM patients in the low- and high-risk group.

In vitro validation of the expression levels of five independent prognostic factors

To validate the results of the public database, we assessed the expression of five distinct prognostic factors in both normal and CM cells using the human fibroblasts cell line HFB4 and human melanoma cell line A375. Our findings indicate that EBLN3P, LINC00324, THCAT158, and USP30-AS1 were significantly overexpressed in HFB4 cells, while MIR205HG was significantly overexpressed in A375 cells (Figure 10).

Relative expression of five independent prognostic factors in HFB4 and A375 cell lines. Relative expression of (A) EBLN3P; (B) LINC00324; (C) MIR205HG; (D) THCAT158; (E) USP30-AS1. Data: mean±SD, *P

Figure 10. Relative expression of five independent prognostic factors in HFB4 and A375 cell lines. Relative expression of (A) EBLN3P; (B) LINC00324; (C) MIR205HG; (D) THCAT158; (E) USP30-AS1. Data: mean±SD, *P<0.05, **P<0.01, ***P<0.001, ****P<0.0001.

Discussion

CM is a malignant skin tumor with a poor prognosis, and its incidence has been increasing in recent years. CM has the ability to metastasize and evade immune and cytotoxic attacks, making conventional therapies, particularly chemotherapy, insufficient [16, 20]. While early-stage CM is usually treated with surgical resection, immunotherapy has shown positive outcomes for inoperable or metastatic cases [21, 22]. Notably, the mechanisms involved in immunotherapy are too broad and individual differences in different patients, and the prognosis of cancer patients through immunotherapy also varies from different population [18, 23]. Thus, exploring CM in-depth and discovering new biomarkers and diagnostic options are critical for improving patient outcomes.

The tumor microenvironment, particularly the hypoxic milieu, is known to play a significant role in the regulation of tumor cells, thereby influencing tumor progression and response to antitumor therapy [24, 25]. A mounting body of evidence indicates that the hypoxic microenvironment exerts a profound impact on tumor progression and response to therapy [26, 27]. The hypoxic microenvironment can facilitate tumor progression, metastasis, and heterogeneity, while also eliciting a range of immunosuppressive responses, such as modulating the levels of macrophages, natural killer cells, and T cells [28, 29]. Long non-coding RNAs (lncRNAs) are recognized to play a pivotal role in the communication between cancer cells and the hypoxic microenvironment.

Long non-coding RNAs (lncRNAs) have emerged as critical players in various aspects of tumorigenesis, particularly as novel biomarkers for tumor diagnosis and prognosis [30]. However, current research on the relationship between lncRNAs, hypoxia, and CM mainly focuses on individual or select lncRNAs, and a systematic study on the prognostic prediction of CM patients using high-risk lncRNA (HRL) signatures is lacking. Therefore, it is crucial to establish HRL signatures based on large-scale databases for accurate prognosis prediction in CM. In this study, we developed a novel risk model based on five HRLs. Univariate and multivariate Cox regression analyses revealed that the risk score derived from the HRL prognostic signature could serve as an independent prognostic factor, distinguishing it from other clinicopathological characteristics. Functional enrichment analysis demonstrated that HRLs may be involved in immune-related processes in CM. Results from immune infiltration analysis revealed that the HRL-based risk model could reflect the immune status of CM patients, providing insights for individualized immunotherapy. Additionally, drug sensitivity analysis revealed differences in response to antineoplastic drugs between the two CM subgroups identified by our risk model. Furthermore, mutation frequency analysis indicated that the low-risk group had higher mutation frequencies. Collectively, our findings provide a comprehensive view for the treatment of CM, contributing to the improvement of patient prognosis.

Recent evidence suggests that long non-coding RNAs (lncRNAs) play a significant role in tumor tumorigenesis, invasion, and metastasis, and their dysregulated expression is associated with the progression and recurrence of various cancers, including cutaneous melanoma (CM). In our study, we found that the expression levels of LINC00324, USP30-AS1, EBLN3P, and THCAT158 were significantly higher in the low-risk group, whereas the high-risk group exhibited higher expression of MIR205HG. LINC00324, a 2115 bp lncRNA located on chromosome 17p13.1, has been found to be abnormally expressed in multiple human cancer types, correlating with tumor initiation and progression [31, 32]. Ding et al. reported that LINC00324 is a protective factor in melanoma patients [33], which is consistent with our findings. USP30-AS1, transcribed from the antisense strand of the USP30 gene, has been implicated in autophagy and is a potential prognostic indicator in cancer, according to Sun et al. [34].

Chen et al. revealed that USP30-AS1 significantly inhibited cell proliferation, migration, and invasion in vitro, as well as tumor growth in vivo [35]. EBLN3P (Endogenous bornavirus-like nucleoprotein) is a recently identified lncRNA located on the forward strand of chromosome 9:37,079,935–37,086,874 [36]. EBLN3P has been found to be prognostic and incorporated into prognostic lncRNA signatures in multiple tumors [37, 38]. Clinical assays have demonstrated that high EBLN3P expression is positively correlated with tumor size, differentiation, and TNM stage, indicating a poor prognosis [39]. Although THCAT158 has not been frequently reported to be associated with tumors, it warrants further investigation. A study indicated that the lncRNA MIR205HG promoted the proliferation, migration, and invasion ability of Hepatoblastoma (HB) [40]. Taken together, these findings suggest that the candidate HRLs are involved in the progression of multiple human tumors. Therefore, it is necessary to construct a risk model based on HRLs to stratify CM patients by risk and predict their prognosis.

There is mounting evidence to suggest that hypoxia promotes tumor growth by suppressing the immune response [41, 42]. The potential involvement of long non-coding RNAs (lncRNAs) in immune regulation cannot be overlooked in hypoxic conditions. In our study, we observed a marked down-regulation of immune-related signaling pathways in high-risk CM patients. Enrichment analysis revealed that immune-related processes may mediate the role of hypoxia-related lncRNAs (HRLs) in CM patients. Furthermore, our immune function analysis demonstrated that Cluster B patients exhibited higher immune function scores compared with Cluster A patients, leading to differences in tumor growth, progression, infiltration, and angiogenesis across risk groups, ultimately resulting in poor prognosis. In summary, the differences in immune function suggest that the HRL signature we established may also reflect the infiltration of immune cells in CM and provide valuable information for clinical immunotherapy. Therefore, the HRL signature is not only a prognostic biomarker but also a predictor of tumor immune status, which can assist clinicians and physicians in better managing patients with CM.

Immunotherapy plays a crucial role in the treatment of various diseases, including cutaneous melanoma (CM). Emerging immunotherapeutic strategies for CM primarily focus on immune checkpoint inhibitors (ICIs), such as PD-1, PD-L1, CTLA4, and LAG3. In our study, we aimed to predict the response of patients with different risk levels to immunotherapy. The findings showed that patients in the high-risk group demonstrated a better response to immunotherapy. Hypoxic TME can directly affect the biological characteristics of invasive immune cells and their response to therapy. Evidence has been provided for possible mechanisms by which the hypoxic state promotes immune escape [43]. Hypoxia allows myeloid-derived suppressor cells (MDSCs) to function and have immunosuppressive activity, allowing cancer cells to evade immune surveillance and resist immune checkpoint blockades [44, 45]. In addition, as a major transcriptional regulator of tumor cell response to hypoxia, HIF-1α can interact with histone deacetylase 1 (HDAC1) and cause immune dysfunction [43, 46]. HIF-1α is also involved in the transformation of tumor-associated macrophages (tam) from antitumor phenotype (M1) to tumor-causing phenotype (M2) [47, 48]. The process of differentiation and functional implementation of T cell can also be inhibited by hypoxia through HIF-1α expression adjustment [49]. The above evidence shows the correlation between hypoxia and immune dysfunction and immune escape. We observed that patients in the high-risk group responded well to immunotherapy. Referring to the strong immunogenicity of CM, a possible explanation is that patients in the high-risk group have a strong dependence on the immune dysfunction and immune escape caused by hypoxia in the process of tumor development, so that the immune function recovery induced by immunotherapy can quickly kill tumor cells. The IPS results suggested that patients in the low-risk group benefited from anti-CTLA4, anti-PD-1, and anti-CTLA4/PD-1 immunotherapy, exhibiting higher expression of LAG3, CTLA4, PD−1, PDCD1LG2, and PD−L1. Lymphocyte-activating gene-3 (LAG-3), as a next-generation immune checkpoint, plays a critical role in regulating immune homeostasis by negatively regulating T cell activation and function [50]. CTLA4, on the other hand, is a potent inhibitor of T-cell proliferation, which prevents overactivation of the immune response [51].

Several studies have suggested that the programmed cell death-1 (PD-1)/programmed cell death ligand 1 (PD-L1) pathway plays a vital role in regulating immune responses. Targeting this pathway has been considered a breakthrough in the treatment of CM [52, 53]. Additionally, patients at different risk levels may have varying susceptibilities to anticancer drugs. High-risk group patients have been found to exhibit higher IC50 levels to pathway inhibitors such as AKT inhibitors, JNK inhibitors, and some drugs approved for anti-tumor treatment, such as Pazopanib, Lapatinib, Crizotinib, among others. These findings suggest that CM patients in different risk subgroups exhibit a promising response to antineoplastic drugs, highlighting the critical need for precisely targeted drugs and chemotherapy for CM.

Conclusions

In this study, we developed a novel risk model based on HRLs to classify CM patients into low- and high-risk groups, which demonstrated robust sensitivity and specificity as a prognostic predictor for CM. Our functional and immune cell infiltration analyses further validated the association of this risk model with immune response. Overall, our study offers unique insights into individualized immunotherapy treatment for CM.

Author Contributions

Q.X. was responsible for the study design and drafted the manuscript. G.W. and Y.S. were involved in data acquisition, analyses and interpretation. Y.S. finished the experimental part of the study. All authors approved the final manuscript.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Funding

This study was supported by the Clinical Research Program of Shandong University (Key Special Project on Acute and Critical Illness) (2021SDUCRCC005).

References

  • 1. 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. https://doi.org/10.3322/caac.21660 [PubMed]
  • 2. Rebecca VW, Somasundaram R, Herlyn M. Pre-clinical modeling of cutaneous melanoma. Nat Commun. 2020; 11:2858. https://doi.org/10.1038/s41467-020-15546-9 [PubMed]
  • 3. Rogala P, Czarnecka AM, Cybulska-Stopa B, Ostaszewski K, Piejko K, Ziętek M, Dziura R, Rutkowska E, Galus Ł, Kempa-Kamińska N, Calik J, Sałek-Zań A, Zemełka T, et al. Long-Term Outcomes of Targeted Therapy after First-Line Immunotherapy in BRAF-Mutated Advanced Cutaneous Melanoma Patients-Real-World Evidence. J Clin Med. 2022; 11:2239. https://doi.org/10.3390/jcm11082239 [PubMed]
  • 4. Thornton J, Chhabra G, Singh CK, Guzmán-Pérez G, Shirley CA, Ahmad N. Mechanisms of Immunotherapy Resistance in Cutaneous Melanoma: Recognizing a Shapeshifter. Front Oncol. 2022; 12:880876. https://doi.org/10.3389/fonc.2022.880876 [PubMed]
  • 5. Manoochehri Khoshinani H, Afshar S, Najafi R. Hypoxia: A Double-Edged Sword in Cancer Therapy. Cancer Invest. 2016; 34:536–45. https://doi.org/10.1080/07357907.2016.1245317 [PubMed]
  • 6. Semenza GL. Oxygen sensing, hypoxia-inducible factors, and disease pathophysiology. Annu Rev Pathol. 2014; 9:47–71. https://doi.org/10.1146/annurev-pathol-012513-104720 [PubMed]
  • 7. Semenza GL. Molecular mechanisms mediating metastasis of hypoxic breast cancer cells. Trends Mol Med. 2012; 18:534–43. https://doi.org/10.1016/j.molmed.2012.08.001 [PubMed]
  • 8. Vaupel P, Multhoff G. Hypoxia-/HIF-1α-Driven Factors of the Tumor Microenvironment Impeding Antitumor Immune Responses and Promoting Malignant Progression. Adv Exp Med Biol. 2018; 1072:171–5. https://doi.org/10.1007/978-3-319-91287-5_27 [PubMed]
  • 9. Hu K, Babapoor-Farrokhran S, Rodrigues M, Deshpande M, Puchner B, Kashiwabuchi F, Hassan SJ, Asnaghi L, Handa JT, Merbs S, Eberhart CG, Semenza GL, Montaner S, Sodhi A. Correction: Hypoxia-inducible factor 1 upregulation of both VEGF and ANGPTL4 is required to promote the angiogenic phenotype in uveal melanoma. Oncotarget. 2021; 12:519–20. https://doi.org/10.18632/oncotarget.27780 [PubMed]
  • 10. Jing X, Yang F, Shao C, Wei K, Xie M, Shen H, Shu Y. Role of hypoxia in cancer therapy by regulating the tumor microenvironment. Mol Cancer. 2019; 18:157. https://doi.org/10.1186/s12943-019-1089-9 [PubMed]
  • 11. Suthen S, Lim CJ, Nguyen PH, Dutertre CA, Lai HL, Wasser M, Chua C, Lim TK, Leow WQ, Loh TJ, Wan WK, Pang YH, Soon G, et al. Hypoxia-driven immunosuppression by Treg and type-2 conventional dendritic cells in HCC. Hepatology. 2022; 76:1329–44. https://doi.org/10.1002/hep.32419 [PubMed]
  • 12. Ouyang J, Zhong Y, Zhang Y, Yang L, Wu P, Hou X, Xiong F, Li X, Zhang S, Gong Z, He Y, Tang Y, Zhang W, et al. Long non-coding RNAs are involved in alternative splicing and promote cancer progression. Br J Cancer. 2022; 126:1113–24. https://doi.org/10.1038/s41416-021-01600-w [PubMed]
  • 13. Rafiee A, Riazi-Rad F, Havaskary M, Nuri F. Long noncoding RNAs: regulation, function and cancer. Biotechnol Genet Eng Rev. 2018; 34:153–80. https://doi.org/10.1080/02648725.2018.1471566 [PubMed]
  • 14. Bhan A, Soleimani M, Mandal SS. Long Noncoding RNA and Cancer: A New Paradigm. Cancer Res. 2017; 77:3965–81. https://doi.org/10.1158/0008-5472.CAN-16-2634 [PubMed]
  • 15. Jin KT, Lu ZB, Lv JQ, Zhang JG. The role of long non-coding RNAs in mediating chemoresistance by modulating autophagy in cancer. RNA Biol. 2020; 17:1727–40. https://doi.org/10.1080/15476286.2020.1737787 [PubMed]
  • 16. Davis LE, Shalin SC, Tackett AJ. Current state of melanoma diagnosis and treatment. Cancer Biol Ther. 2019; 20:1366–79. https://doi.org/10.1080/15384047.2019.1640032 [PubMed]
  • 17. Lázari LC, Wolf IR, Schnepper AP, Valente GT. LncRNAs of Saccharomyces cerevisiae bypass the cell cycle arrest imposed by ethanol stress. PLoS Comput Biol. 2022; 18:e1010081. https://doi.org/10.1371/journal.pcbi.1010081 [PubMed]
  • 18. Attrill GH, Ferguson PM, Palendira U, Long GV, Wilmott JS, Scolyer RA. The tumour immune landscape and its implications in cutaneous melanoma. Pigment Cell Melanoma Res. 2021; 34:529–49. https://doi.org/10.1111/pcmr.12926 [PubMed]
  • 19. Yu G, Wang LG, Han Y, He QY. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS. 2012; 16:284–7. https://doi.org/10.1089/omi.2011.0118 [PubMed]
  • 20. Gagliardi M, Cotella D, Santoro C, Corà D, Barlev NA, Piacentini M, Corazzari M. Aldo-keto reductases protect metastatic melanoma from ER stress-independent ferroptosis. Cell Death Dis. 2019; 10:902. https://doi.org/10.1038/s41419-019-2143-7 [PubMed]
  • 21. Costa Svedman F, Das I, Tuominen R, Darai Ramqvist E, Höiom V, Egyhazi Brage S. Proliferation and Immune Response Gene Signatures Associated with Clinical Outcome to Immunotherapy and Targeted Therapy in Metastatic Cutaneous Malignant Melanoma. Cancers (Basel). 2022; 14:3587. https://doi.org/10.3390/cancers14153587 [PubMed]
  • 22. Ferrara M, Adams W, Kotagiri A. Melanoma-associated retinopathy after starting immunotherapy for metastatic cutaneous melanoma. Eye (Lond). 2022; 36:2361–2. https://doi.org/10.1038/s41433-022-02057-8 [PubMed]
  • 23. Curran MA, Montalvo W, Yagita H, Allison JP. PD-1 and CTLA-4 combination blockade expands infiltrating T cells and reduces regulatory T and myeloid cells within B16 melanoma tumors. Proc Natl Acad Sci USA. 2010; 107:4275–80. https://doi.org/10.1073/pnas.0915174107 [PubMed]
  • 24. Ahmed F, Haass NK. Microenvironment-Driven Dynamic Heterogeneity and Phenotypic Plasticity as a Mechanism of Melanoma Therapy Resistance. Front Oncol. 2018; 8:173. https://doi.org/10.3389/fonc.2018.00173 [PubMed]
  • 25. Falcone I, Conciatori F, Bazzichetto C, Ferretti G, Cognetti F, Ciuffreda L, Milella M. Tumor Microenvironment: Implications in Melanoma Resistance to Targeted Therapy and Immunotherapy. Cancers (Basel). 2020; 12:2870. https://doi.org/10.3390/cancers12102870 [PubMed]
  • 26. Jha J, Singh MK, Singh L, Pushker N, Lomi N, Meel R, Chosdol K, Sen S, Bakhshi S, Kashyap S. Association of TYRP1 with hypoxia and its correlation with patient outcome in uveal melanoma. Clin Transl Oncol. 2021; 23:1874–84. https://doi.org/10.1007/s12094-021-02597-7 [PubMed]
  • 27. Malekan M, Ebrahimzadeh MA, Sheida F. The role of Hypoxia-Inducible Factor-1alpha and its signaling in melanoma. Biomed Pharmacother. 2021; 141:111873. https://doi.org/10.1016/j.biopha.2021.111873 [PubMed]
  • 28. Yuen VW, Wong CC. Hypoxia-inducible factors and innate immunity in liver cancer. J Clin Invest. 2020; 130:5052–62. https://doi.org/10.1172/JCI137553 [PubMed]
  • 29. Zheng S, Zou Y, Liang JY, Xiao W, Yang A, Meng T, Lu S, Luo Z, Xie X. Identification and validation of a combined hypoxia and immune index for triple-negative breast cancer. Mol Oncol. 2020; 14:2814–33. https://doi.org/10.1002/1878-0261.12747 [PubMed]
  • 30. Li J, Li Z, Leng K, Xu Y, Ji D, Huang L, Cui Y, Jiang X. ZEB1-AS1: A crucial cancer-related long non-coding RNA. Cell Prolif. 2018; 51:e12423. https://doi.org/10.1111/cpr.12423 [PubMed]
  • 31. Dong Y, Wan G, Yan P, Qian C, Li F, Peng G. Long noncoding RNA LINC00324 promotes retinoblastoma progression by acting as a competing endogenous RNA for microRNA-769-5p, thereby increasing STAT3 expression. Aging (Albany NY). 2020; 12:7729–46. https://doi.org/10.18632/aging.103075 [PubMed]
  • 32. Zou Z, Ma T, He X, Zhou J, Ma H, Xie M, Liu Y, Lu D, Di S, Zhang Z. Long intergenic non-coding RNA 00324 promotes gastric cancer cell proliferation via binding with HuR and stabilizing FAM83B expression. Cell Death Dis. 2018; 9:717. https://doi.org/10.1038/s41419-018-0758-8 [PubMed]
  • 33. Ding Y, Li T, Li M, Tayier T, Zhang M, Chen L, Feng S. A Novel Autophagy-Related lncRNA Gene Signature to Improve the Prognosis of Patients with Melanoma. Biomed Res Int. 2021; 2021:8848227. https://doi.org/10.1155/2021/8848227 [PubMed]
  • 34. Sun Z, Jing C, Xiao C, Li T. An autophagy-related long non-coding RNA prognostic signature accurately predicts survival outcomes in bladder urothelial carcinoma patients. Aging (Albany NY). 2020; 12:15624–37. https://doi.org/10.18632/aging.103718 [PubMed]
  • 35. Chen M, Chi Y, Chen H, Zhao L. Long non-coding RNA USP30-AS1 aggravates the malignant progression of cervical cancer by sequestering microRNA-299-3p and thereby overexpressing PTP4A1. Oncol Lett. 2021; 22:505. https://doi.org/10.3892/ol.2021.12766 [PubMed]
  • 36. Dai S, Li N, Zhou M, Yuan Y, Yue D, Li T, Zhang X. LncRNA EBLN3P promotes the progression of osteosarcoma through modifying the miR-224-5p/Rab10 signaling axis. Sci Rep. 2021; 11:1992. https://doi.org/10.1038/s41598-021-81641-6 [PubMed]
  • 37. Gong Z, Li Q, Li J, Xie J, Wang W. A novel signature based on autophagy-related lncRNA for prognostic prediction and candidate drugs for lung adenocarcinoma. Transl Cancer Res. 2022; 11:14–28. https://doi.org/10.21037/tcr-21-1554 [PubMed]
  • 38. Mathias C, Muzzi JC, Antunes BB, Gradia DF, Castro MA, Carvalho de Oliveira J. Unraveling Immune-Related lncRNAs in Breast Cancer Molecular Subtypes. Front Oncol. 2021; 11:692170. https://doi.org/10.3389/fonc.2021.692170 [PubMed]
  • 39. Xu XH, Song W, Li JH, Huang ZQ, Liu YF, Bao Q, Shen ZW. Long Non-coding RNA EBLN3P Regulates UHMK1 Expression by Sponging miR-323a-3p and Promotes Colorectal Cancer Progression. Front Med (Lausanne). 2021; 8:651600. https://doi.org/10.3389/fmed.2021.651600 [PubMed]
  • 40. Yin L, Zhang Y, Zheng L. Analysis of differentially expressed long non-coding RNAs revealed a pro-tumor role of MIR205HG in cervical cancer. Mol Med Rep. 2022; 25:42. https://doi.org/10.3892/mmr.2021.12558 [PubMed]
  • 41. Rama I, Bruene B, Torras J, Koehl R, Cruzado JM, Bestard O, Franquesa M, Lloberas N, Weigert A, Herrero-Fresneda I, Gulias O, Grinyó JM. Hypoxia stimulus: An adaptive immune response during dendritic cell maturation. Kidney Int. 2008; 73:816–25. https://doi.org/10.1038/sj.ki.5002792 [PubMed]
  • 42. Castillo-Rodríguez RA, Trejo-Solís C, Cabrera-Cano A, Gómez-Manzo S, Dávila-Borja VM. Hypoxia as a Modulator of Inflammation and Immune Response in Cancer. Cancers (Basel). 2022; 14:2291. https://doi.org/10.3390/cancers14092291 [PubMed]
  • 43. Zhuang Y, Liu K, He Q, Gu X, Jiang C, Wu J. Hypoxia signaling in cancer: Implications for therapeutic interventions. MedComm (2020). 2023; 4:e203. https://doi.org/10.1002/mco2.203 [PubMed]
  • 44. Noman MZ, Desantis G, Janji B, Hasmim M, Karray S, Dessen P, Bronte V, Chouaib S. PD-L1 is a novel direct target of HIF-1α, and its blockade under hypoxia enhanced MDSC-mediated T cell activation. J Exp Med. 2014; 211:781–90. https://doi.org/10.1084/jem.20131916 [PubMed]
  • 45. He Q, Liu M, Huang W, Chen X, Zhang B, Zhang T, Wang Y, Liu D, Xie M, Ji X, Sun M, Tian D, Xia L. IL-1β-Induced Elevation of Solute Carrier Family 7 Member 11 Promotes Hepatocellular Carcinoma Metastasis Through Up-regulating Programmed Death Ligand 1 and Colony-Stimulating Factor 1. Hepatology. 2021; 74:3174–93. https://doi.org/10.1002/hep.32062 [PubMed]
  • 46. Blouw B, Song H, Tihan T, Bosze J, Ferrara N, Gerber HP, Johnson RS, Bergers G. The hypoxic response of tumors is dependent on their microenvironment. Cancer Cell. 2003; 4:133–46. https://doi.org/10.1016/s1535-6108(03)00194-6 [PubMed]
  • 47. Zhang J, Zhang Q, Lou Y, Fu Q, Chen Q, Wei T, Yang J, Tang J, Wang J, Chen Y, Zhang X, Zhang J, Bai X, Liang T. Hypoxia-inducible factor-1α/interleukin-1β signaling enhances hepatoma epithelial-mesenchymal transition through macrophages in a hypoxic-inflammatory microenvironment. Hepatology. 2018; 67:1872–89. https://doi.org/10.1002/hep.29681 [PubMed]
  • 48. Almendros I, Wang Y, Becker L, Lennon FE, Zheng J, Coats BR, Schoenfelt KS, Carreras A, Hakim F, Zhang SX, Farré R, Gozal D. Intermittent hypoxia-induced changes in tumor-associated macrophages and tumor malignancy in a mouse model of sleep apnea. Am J Respir Crit Care Med. 2014; 189:593–601. https://doi.org/10.1164/rccm.201310-1830OC [PubMed]
  • 49. McNamee EN, Korns Johnson D, Homann D, Clambey ET. Hypoxia and hypoxia-inducible factors as regulators of T cell development, differentiation, and function. Immunol Res. 2013; 55:58–70. https://doi.org/10.1007/s12026-012-8349-8 [PubMed]
  • 50. Grebinoski S, Zhang Q, Cillo AR, Manne S, Xiao H, Brunazzi EA, Tabib T, Cardello C, Lian CG, Murphy GF, Lafyatis R, Wherry EJ, Das J, et al. Autoreactive CD8+ T cells are restrained by an exhaustion-like program that is maintained by LAG3. Nat Immunol. 2022; 23:868–77. https://doi.org/10.1038/s41590-022-01210-5 [PubMed]
  • 51. Lanz AL, Riester M, Peters P, Schwerd T, Lurz E, Hajji MS, Rohlfs M, Ley-Zaporozhan J, Walz C, Kotlarz D, Klein C, Albert MH, Hauck F. Abatacept for treatment-refractory pediatric CTLA4-haploinsufficiency. Clin Immunol. 2021; 229:108779. https://doi.org/10.1016/j.clim.2021.108779 [PubMed]
  • 52. Gao Y, Li S, Xu D, Chen S, Cai Y, Jiang W, Zhang X, Sun J, Wang K, Chang B, Wang F, Hong M. Prognostic value of programmed death-1, programmed death-ligand 1, programmed death-ligand 2 expression, and CD8(+) T cell density in primary tumors and metastatic lymph nodes from patients with stage T1-4N+M0 gastric adenocarcinoma. Chin J Cancer. 2017; 36:61. https://doi.org/10.1186/s40880-017-0226-3 [PubMed]
  • 53. Tong G, Cheng B, Li J, Wu X, Nong Q, He L, Li X, Li L, Wang S. MACC1 regulates PDL1 expression and tumor immunity through the c-Met/AKT/mTOR pathway in gastric cancer cells. Cancer Med. 2019; 8:7044–54. https://doi.org/10.1002/cam4.2542 [PubMed]
Download PDF
Corresponding Author

Qingjia Xu
xuqingjia@qiluhospital.com

Keywords
cutaneous melanoma hypoxia-related long noncoding RNAs consensus clustering analysis immune infiltration landscape prognosis
Paper Sections