ORIGINAL RESEARCH article

Front. Oncol., 30 July 2021

Sec. Skin Cancer

Volume 11 - 2021 | https://doi.org/10.3389/fonc.2021.704543

N6-Methyladenosine Regulators and Related LncRNAs Are Potential to be Prognostic Markers for Uveal Melanoma and Indicators of Tumor Microenvironment Remodeling

  • 1. School of Biomedical Engineering, Capital Medical University, Beijing, China

  • 2. Beijing Key Laboratory of Fundamental Research on Biomechanics in Clinical Application, Capital Medical University, Beijing, China

Abstract

Uveal melanoma (UM) is one of the most common malignant intraocular tumors in adults. Few studies have investigated the effect of N6-methyladenosine (m6A) RNA methylation regulators and related long noncoding RNAs (lncRNAs) on the tumor microenvironment (TME) and survival time of patients with UM. Based on the transcriptome and clinical data from The Cancer Genome Atlas, we systematically identified m6A regulators. Then, we constructed an m6A regulators-based signature to predict the prognostic risk using univariate and LASSO Cox analyses. The signature was then validated by performing Kaplan-Meier, and receiver operating characteristic analyses. Through the correlation analysis, m6A regulators-related lncRNAs were identified, and they were divided into different clustering subtypes according to their expression. We further assessed differences in TME scores, the survival time of patients, and immune cell infiltration levels between different clustering subtypes. Finally, we screened out the common immune genes shared by m6A-related lncRNAs and determined their expression in different risk groups and clustering subtypes. For further validation, we used single-cell sequencing data from the GSE139829 dataset to explore the expression distribution of immune genes in the TME of UM. We constructed a prognostic risk signature representing an independent prognostic factor for UM using 3 m6A regulators. Patients in the low-risk group exhibited a more favorable prognosis and lower immune cell infiltration levels than patients in the high-risk group. Two subtypes (cluster 1/2) were identified based on m6A regulators-related lncRNAs. The TME scores, prognosis, and immune cell infiltration have a marked difference between cluster 1 and cluster 2. Additionally, 13 common immune genes shared by 5 lncRNAs were screened out. We found that these immune genes were differentially expressed in different risk groups and clustering subtypes and were widely distributed in 3 cell types of TME. In conclusion, our study demonstrated the important role of m6A regulators and related lncRNAs in TME remodeling. The signature developed using m6A regulators might serve as a promising parameter for the clinical prediction of UM.

Introduction

Uveal melanoma (UM), which is the second most common type of melanoma, originates from melanocytes in the intraocular uvea (). Currently, surgery and radiotherapy are the most effective methods to treat local tumors (). However, the overall mortality in patients with UM is more than 50%, because it is highly susceptible to early metastasis (, ). Therefore, new treatments such as immunotherapy or targeted therapy (, ) are being developed, which requires the identification of several potential prognostic biomarkers and therapeutic targets for UM.

The most prevalent RNA modification is N6-methyladenosine (m6A), which involves methylation of the sixth N atom of adenine (). m6A methylation is a dynamic process regulated by methyltransferases (writers) and demethylases (erasers), whereas binding proteins (readers) bind to m6A methylation sites (, ). Methyltransferases such as METTL3/14/16, WTAP, RBM15, and VIRMA promote m6A methylation (). On the other hand, demethylases, which include FTO and ALKBH5, inhibit m6A methylation (, ). Binding proteins, such as YTHDC1/2, YTHDF1-3, IGF2BP1-3, and HNRNPC, bind to the m6A modified site to form a complex that mediates its biological function (). These m6A methylases are primarily involved in mammalian development, immune response, tumorigenesis, and metastasis, and stem cell differentiation (). However, the prognostic role of m6A methylases in UM development has not been sufficiently investigated. Besides, the involvement of m6A methylases in the tumor microenvironment (TME) remains to be thoroughly explored.

Long noncoding RNAs (lncRNAs) are frequently defined as RNAs that have a transcript length exceeding 200 nucleotides and do not encode proteins (). They can regulate gene expression at epigenetic, transcriptional, and post-transcriptional levels (). A recent study has reported that lncRNAs promote tumor development by altering the immune microenvironment (). Increasing evidence has demonstrated the TME, which mainly consists of stromal and immune cells, plays an important role in tumor progression (). Stromal cells may contribute to tumor angiogenesis and extracellular matrix reorganization, whereas immune cells may contribute to TME via dysregulation of immune-mediated responses (). Therefore, immune cell infiltration in the TME may serve as a potential target for immunotherapy. However, the involvement of lncRNAs in immune cell infiltration in UM remains unclear.

This study aimed to systematically explore m6A regulators and related lncRNAs involved in the TME in UM and developed an m6A regulators-based signature for improving the accuracy of prognosis in patients with UM. We also established clustering subtypes based on m6A regulators-related lncRNAs to determine the relationships between the clustering subtypes, TME scores, prognosis, and immune cell infiltration, and further explained the mechanism of action of m6A regulators. Finally, we explored the expression of 13 immune genes shared by 5 lncRNAs in different risk groups and clustering subtypes.

Materials and Methods

Data Collection and Preparation

The RNA-sequencing transcriptome data of 80 patients with UM and corresponding clinical data were downloaded from The Cancer Genome Atlas (TCGA) data portal (https://portal.gdc.cancer.gov/). GTF files were downloaded from Ensembl (https://asia.ensembl.org) to distinguish between lncRNAs and mRNAs for subsequent analyses. The list of immune genes was downloaded from the ImmPort database (https://www.immport.org).

Generation of TME Scores and Tumor-Infiltrating Immune Cells

The ESTIMATE algorithm in the R “estimate” package was used to calculate the TME scores of 80 UM patients. The Kaplan-Meier survival analysis was conducted to compare the difference in survival time using R “survMiner” and “survival” packages. The fraction of 22 immune cell types in each sample was estimated using CIBERSORT. The association between TME scores and tumor-infiltrating immune cells (TICs) was established using correlation analysis.

Construction and Validation of m6A Regulators-Based Signature

The m6A regulators were identified from the published literature (). The m6A regulators contain 8 writers (METTL3, METTL14, METTL16, WTAP, VIRMA, ZC3H13, RBM15, and RBM15B), 13 readers (YTHDC1, YTHDC2, YTHDF1, YTHDF2, YTHDF3, HNRNPC, FMR1, LRPPRC, HNRNPA2B1, IGF2BP1, IGF2BP2, IGF2BP3, and RBMX), and 2 erasers (FTO and ALKBH5). The expression data of m6A regulators were extracted from the mRNA expression data of TCGA. The m6A regulators, which were previously identified using the univariate Cox regression analysis, were further subjected to the LASSO Cox regression analysis using the “glmnet” package. The minimum 10-fold cross-validation was used to select the best penalty parameter λ. Then, the risk score of each patient was calculated using a linear combination of m6A regulators expression weighted by the multivariate Cox regression analysis. According to the median risk score, the samples were divided into high-risk and low-risk groups. Subsequently, the Kaplan-Meier survival analysis was performed to compare the survival difference. Receiver operating characteristic (ROC) analysis was performed to evaluate the prognostic value of the signature using the “timeROC” package.

Generation and Validation of Clustering Subtypes From m6A Regulators-Related LncRNAs

We screened m6A regulators-related lncRNAs by Pearson’s correlation analysis. The process used the criteria of |correlation coefficient| > 0.4 and p < 0.001. The expression data of m6A regulators-related lncRNAs were extracted from the lncRNA expression data of TCGA. To clarify the biological characteristics of m6A regulators-related lncRNAs, the R “ConsensusClusterPlus” package was used to divide the samples into different clustering subtypes according to the expression of lncRNAs. Principal component analysis (PCA), Kaplan-Meier survival analysis, TME scores, and TICs profiles were performed for different clustering subtypes.

Identification and Validation of Immune Genes Shared by 5 LncRNAs

The Pearson’s correlation analysis was performed to screen common immune genes shared by 5 lncRNAs. The Kaplan-Meier survival analysis was used to compare the survival difference between high- and low-expression of immune genes. The differential expression analysis of immune genes was performed using the R “limma” package in different risk groups and clustering subtypes. A p-value < 0.05 indicated statistical significance. At present, single-cell sequencing technology has been widely used to explore the heterogeneity of TME. To characterize immune genes expression distribution in TME of UM, we search for single-cell sequencing data of UM from Tumor Immune Single-cell Hub (TISCH) (). The GSE139829 dataset including 59,915 cell sequencing data from 11 samples was collected to perform gene expression distribution ().

Results

The Correlation Between TME Scores With the Survival of UM Patients and Immune Cell Infiltration

To establish the correlation of ImmuneScore and StromalScore with the survival time, we performed the Kaplan-Meier survival analysis. A high score of immune and stromal cells signified large numbers of these cells in the TME. As shown in Figure 1A, the overall survival (OS) in the low ImmuneScore group was longer than that in the high ImmuneScore group. Similarly, StromalScore and ESTIMATEScore showed a negative correlation with the OS (Figures 1B, C). To confirm exact changes in the genetic profiles in the TME about immune and stromal cell components, variance analysis of high and low scores was performed. As shown in Figures 1D, E, 700 common differentially expressed genes (DEGs) were upregulated, and 74 common DEGs were downregulated. To further explore the interaction between TME scores and the 22 immune cell types, we first estimated the 22 types of TICs with abundance distribution in all the tumor samples and then calculated the correlation index between the TME scores and TICs. The results showed that ImmuneScore is correlated with CD8+ T cells, resting memory CD4+ T cells, activated memory CD4+ T cells, helper T cells (follicular), Tregs, resting NK cells, monocytes, M0 macrophages, M1 macrophages, resting mast cells, and eosinophils (Figure 1F).

Figure 1

Construction and Validation of m6A Regulators-Based Signature

To clarify the biological function of m6A regulators in the prognosis of patients with UM, we comprehensively investigated the prognostic value of m6A regulators based on the expression and clinical data (Table S1). Of these m6A regulators, seven exhibited a prognostic value based on OS (Figure 2A), while nine displayed a prognostic value based on progression-free survival (PFS) (Figure 2B). The Venn plot results indicated that 6 m6A regulators (RBM15B, IGF2BP2, YTHDF1, METTL16, VIRMA, and YTHDF3) were identified based on OS and PFS (Figure 2C). To avoid overfitting, we performed the LASSO Cox analysis and selected 3 of the 6 m6A regulators were to establish a risk signature (Figures 2D, E). Thus, we established a predictive model: risk score = (RBM15B * −0.14284) + (YTHDF3 * 0.02121) + (IGF2BP2 * −0.11533). The distribution of risk score (Upper), patients’ survival time (Middle), and heat map analysis (Bottom) of the 3 prognostic m6A regulators were shown based on the OS (Figures 2F) and PFS (Figure 2G). Results of the heat maps (Figures 2F, G) was survival curves (Figures 3A–C) suggested that YTHDF3 was likely to be a high-risk factor because it was upregulated in the high-risk group. However, the highly expressed RBM15B and IGF2BP2 in the low-risk group might be protective factors. As shown in Figures 3D, E, the OS and PFS of patients in the low-risk group were longer than in the high-risk group. To evaluate the prognostic accuracy of the 3 m6A regulators-based signature, we performed the ROC analysis based on OS and PFS. Areas under the ROC curves of 1, 2, and 3 years were 0.774, 0.811, and 0.843, respectively (Figure 4A). The PFS prediction of 3 m6A regulators-based signature was also accurate (Figure 4B). These results suggested that the 3 m6A regulators-based signature might serve as a promising parameter for prognostic prediction of UM.

Figure 2

Figure 3

Figure 4

The Association Between Risk Score With Immune Cell Infiltration and TME Scores

The correlation between the risk score and the immune cell infiltration levels was calculated to establish the association between the 3 m6A regulators-based signature with the TME. ImmuneScore, StromalScore, and ESTIMATEScore in the high-risk group were significantly higher than in the low-risk group (Figures 5A–C). The results showed that the risk score was significantly negatively correlated with the infiltration levels of naive B cells (Figure 5D), eosinophils (Figure 5F), monocytes (Figure 5G), and plasma cells (Figure 5H). Only the infiltration level of resting dendritic cells was positively correlated with the risk score (Figure 5E). These results indicated that the 3 m6A regulators were involved in the immune cell infiltration of UM.

Figure 5

The Correlation Between Consensus Clustering of m6A Regulators-Related LncRNAs and TME Scores, Survival Time, and Immune Cell Infiltration

Using the correlation analysis, we identified 514 lncRNAs based on RNA-Seq data and constructed a network between m6A regulators and lncRNAs (Figure 6A). A total of 66 of the 514 lncRNAs had a prognostic value based on the OS (Figure S1A), while 70 of the 514 lncRNAs had a prognostic value based on the PFS (Figure S1B). The Venn plot results showed that 38 lncRNAs were common according to OS and PFS (Figure 6B). Based on the similarity identified by consensus clustering using the ‘ConsensusClusterPlus’ package, we found that k = 2 was the optimal clustering stability value (Figures 6C, D). The 80 UM samples were well-differentiated into two subtypes according to the expression of 38 lncRNAs (Figure 6E).

Figure 6

The TME scores of cluster 1 were lower than those of cluster 2 (Figures 7A–C), while the OS (Figure 7D) and PFS (Figure 7E) of cluster 1 were notably longer than those of cluster 2. Subsequently, the 22 immune cell levels for the two subtypes were calculated. The results showed that cluster 1 had higher immune infiltration levels of plasma cells, and monocytes while there were higher immune cell infiltration levels of activated memory CD4+ T cells, follicular helper T cells, M1 macrophages, and resting dendritic cells (Figure 7F).

Figure 7

Identification and Validation of Immune Genes Shared by 5 LncRNAs

First, we identified the common lncRNAs from 5 lists that were DEGs in ImmuneScore (Immune-DEGs) and StromalScore (Stromal-DEGs), m6A regulators-related lncRNAs (m6A-lncRNAs), and immune gene-specific lncRNAs based on OS (OS-immune) and PFS (PFS-immune). A total of 5 lncRNAs, namely AC008555.4, AC018529.1, AC104129.1, CYTOR, and MIR4435-2HG, were found to be common across the 5 lists (Figure 8A). High expression of AC018529.1 (Figure 8B), MIR4435-2HG (Figure 8C), AC104129.1 (Figure 8D), and CYTOR (Figure 8F) was related to the poor prognosis of patients, while high expression of AC008555.4 was associated with good prognosis (Figure 8E). Next, 13 immune genes having 5 common lncRNAs were screened using the correlation analysis (Figure 9A). The 13 immune genes were ADGRE5, C2, CD79B, CTSC, GEM, JAG2, LYN, MAFB, MBP, MR1, PREX1, RUNX1, and TCF12. Except for C2, the other 12 immune genes were associated with a poor prognosis (Figures 9B–M). To further verify whether the immune genes were differentially expressed in different groups based on m6A regulators and lncRNAs, we extracted the expression data from RNA-Seq data and performed the differential expression analysis. The expression levels of 13 immune genes in the high-risk group were upregulated than those in the low-risk group (Figure 10). Similarly, the 13 immune genes levels in cluster 2 were upregulated compared with cluster 1 (Figure 11). For further validation, we used TISCH to depict the expression distribution of 13 immune genes in the TME of UM. The overall distribution of the 3 cell types in the GSE139829 dataset was shown in Figure 12A. Through analysis, we found that the expression distribution of ADGRE5 (Figure 12B), CTSC (Figure 12E), LYN (Figure 12H), MAFB (Figure 12I), PREX1 (Figure 12L), and RUNX1 (Figure 12M) were abundant in immune cells. CD79B (Figure 12D), GEM (Figure 12F), and JAG2 (Figure 12G) expression distribution were concentrated in malignant cells. The gene expression of C2 (Figure 12C), MBP (Figure 12J), MR1 (Figure 12K), and TCF12 (Figure 12N) were evenly distributed in immune cells and malignant cells. These results suggested that these immune genes may be the downstream regulators of m6A regulators and related lncRNAs participated in TME remodeling.

Figure 8

Figure 9

Figure 10

Figure 11

Figure 12

Discussion

The TME plays a key role in different stages of tumorigenesis. Eyes are an immune-privileged site but inflammation can develop in an ocular tumor TME (). UM is homogeneous without much stromal tissue, and therefore, it may be affected by immune cells (). Compared with other malignancies, the presence of infiltrating macrophages and T cells in UM is associated with a poorer rather than a better prognosis (), which was consistent with our findings. Moreover, previous studies have suggested that tumor-infiltrating macrophages and T cells are independent predictors for the prognosis of patients with UM (, ). In this study, results of the transcriptome analysis of UM data indicated that UM patients with high ImmuneScore had a poor prognosis. Besides, ImmuneScore was found to be significantly associated with many TICs such as T cells and macrophages. These results suggested that the TME played an important role in UM. Clarifying the mechanisms of the TME will provide novel insight into the development of highly effective immunotherapeutic strategies.

Post-transcriptional regulation is important for regulating the gene expression processes, which determine cellular function. Decades of research have identified more than 100 types of ribonucleosides that are post-transcriptionally modified (). m6A methylation is one of the most prevalent post-transcriptional modifications found in eukaryotic mRNAs and lncRNAs (, ). More studies have reported that m6A regulators extensively participate in diverse biological processes and prognoses in different cancers (, , , ). A recent study has suggested that METTL3-mediated m6A methylation modulates UM cell proliferation, migration, and invasion by targeting c-Met (). As far as we know, the role of m6A methylation in UM has less been studied, and the effect of m6A methylation on the TME of UM has not been fully understood.

In this study, we found that m6A regulators were related to the prognosis and TME of patients with UM. We established a prognostic risk signature using 3 m6A regulators based on OS. The signature helped differentiate UM patients into high- and low-risk groups and could serve as an independent risk factor for UM prognosis. The high-risk group was positively correlated with immune cell infiltration levels. Among the 3 m6A regulators, IGF2BP2 acts as m6A readers to enhance mRNA stability and translation and plays an important role in tumors (). YTHDF3 functions as oncogenes in breast cancer (). A recent study has investigated that ocular melanoma samples show decreased m6A levels, indicating a poor prognosis (). In our study, patients with high RBM15B expression had a good prognosis. RBM15B acts as a methyltransferase and thus promotes the level of m6A RNA methylation. Therefore, it is reasonable to speculate that high levels of m6A methylation are beneficial to patients’ survival, which is consistent with the current study. These findings indicated that m6A regulators played an important role in the development and progression of cancer. However, the underlying mechanisms of m6A in tumor development still need to be further clarified.

Currently, some studies have clarified the role of lncRNAs in UM. Among these lncRNAs, lncRNA PVT1 and R2RX7-V3 function as novel oncogenes and promote tumorigenesis (, ), whereas lncRNAs CANT1 and PAUPAR suppress tumorigenesis in malignant UM (, ). However, no study has analyzed the effect of m6A regulators-related lncRNAs on the TME and prognosis in UM. Here we identified m6A regulators-related lncRNAs by performing the correlation analysis and further screened lncRNAs based on OS and PFS. The cluster 1/2 subtypes identified through consensus clustering based on the expression of m6A regulators-related lncRNAs were also related to ImmuneScore, the prognosis of patients, and immune cell infiltration levels. Finally, 5 m6A regulators-related lncRNAs were found to be associated with the OS of UM patients. Among the 5 lncRNAs, lncRNA MIR4435-2HG targets desmoplakin and promotes growth and metastasis of gastric cancer by activating Wnt/β-catenin signaling (). LncRNA-CYTOR and Wnt/β-Catenin signaling form a positive feed-forward loop to promote the metastasis of colon cancer (). Through the correlation analysis, we screened 13 downstream immune genes targeted by 5 lncRNAs. They were all found to be upregulated in the high-risk group and cluster 2. At present, RUNX1, MR1, and PREX1 have been reported to be associated with T cells (). CD79B is an important driver of immune-privileged site-associated diffuse large B cell lymphoma (). JAG2 has been found to be overexpressed in malignant plasma cells from multiple myeloma patients and cell lines (). These results suggested that lncRNAs may affect the immune cell infiltration through 13 common immune genes. Finally, we constructed a Sankey diagram that depicted the relationship between m6A methylation regulators, lncRNAs, and immune genes (Figure 13). These findings require further validation and may provide invaluable insights into the future treatment of patients with UM.

Figure 13

m6A methylation is a prevalent form of RNA modification that may provide a novel approach for tumor treatment. However, many key aspects, such as the regulatory mechanisms of m6A regulators and the unidentified relationship between m6A regulators and TME, remain to be explored. Therefore, in this study, we systematically explored the relationship between m6A regulators with the prognosis and TME in UM, further identified the potential lncRNAs and immune genes. However, further validation based on more clinical samples is required, and thus clinical samples will be collected to determine the level of m6A methylation and the association between the expression of m6A regulators and patients’ survival in the future. Furthermore, the downstream regulatory mechanisms of m6A regulators will be investigated to screen possible targets by methylated RNA immunoprecipitation sequencing (MeRIP-seq) and RNA-binding protein immunoprecipitation-quantitative polymerase chain reaction (RIP-qPCR). Tumorigenesis in animals and phenotypes of cell lines are necessary to explore the function of m6A regulators in this process.

In conclusion, our study provided an m6A regulators-based signature for prognostic prediction of UM and confirmed that m6A regulators and related lncRNAs played an important role in TME remodeling. These findings might provide promising targets for improving the survival of UM patients.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Statements

Data availability statement

The datasets analysed for this study can be found in The Cancer Genome Atlas (https://portal.gdc.cancer.gov/).

Author contributions

ZL (1st author) and ZL (5th author) conceived the project. ZL (1st author) performed the computational analyses and wrote the manuscript. SL designed and produced the figures. SH and TW contributed to the literature search for the manuscript. SL and ZL (5th author) revised the manuscript. All authors contributed to the article and approved the submitted version.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fonc.2021.704543/full#supplementary-material

Supplementary Figure 1

Identification of m6A regulators-related lncRNAs related to OS and PFS. (A, B) Forest plots for the univariate Cox analysis of prognosis based on OS (A) and PFS (B). Colored dots represent hazard ratio, and the horizontal lines across the hazard ratio represent 95% confidence interval. OS, overall survival; PFS, progression free survival.

References

Summary

Keywords

m6A RNA methylation regulators, long noncoding RNAs, tumor microenvironment, immune cell infiltration, uveal melanoma

Citation

Liu Z, Li S, Huang S, Wang T and Liu Z (2021) N6-Methyladenosine Regulators and Related LncRNAs Are Potential to be Prognostic Markers for Uveal Melanoma and Indicators of Tumor Microenvironment Remodeling. Front. Oncol. 11:704543. doi: 10.3389/fonc.2021.704543

Received

03 May 2021

Accepted

14 July 2021

Published

30 July 2021

Volume

11 - 2021

Edited by

Vijayasaradhi Setaluri, University of Wisconsin-Madison, United States

Reviewed by

Chandra K. Singh, University of Wisconsin-Madison, United States; Sandhya Annamaneni, Osmania University, India

Updates

Copyright

*Correspondence: Zhicheng Liu,

This article was submitted to Skin Cancer, a section of the journal Frontiers in Oncology

Disclaimer

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.

Outline

Figures

Cite article

Copy to clipboard


Export citation file


Share article

Article metrics