Comprehensive analysis of N6-methyladenosine regulators with the tumor immune landscape and correlation between the insulin-like growth factor 2 mRNA-binding protein 3 and programmed death ligand 1 in bladder cancer

doi:10.21203/rs.3.rs-853237/v1

Download PDF

Research Article

Comprehensive analysis of N⁶-methyladenosine regulators with the tumor immune landscape and correlation between the insulin-like growth factor 2 mRNA-binding protein 3 and programmed death ligand 1 in bladder cancer

https://doi.org/10.21203/rs.3.rs-853237/v1

This work is licensed under a CC BY 4.0 License

Journal Publication

published 11 Feb, 2022

Read the published version in Cancer Cell International →

You are reading this latest preprint version

Background

N⁶-methyladenosine (m⁶A) is one of the most abundant post-transcriptional modifications of RNA. However, there is limited information about the potential roles of m⁶A regulators in tumor immunity. Therefore, in this study, we aimed to determine the functions of m⁶A regulators in bladder cancer as well as their association with the tumor immune landscape.

Methods

We reported that the variation and expression levels of m6A regulators in TCGA database and GTEx database of bladder cancer. Clusters, risk score patterns, and nomograms were constructed to evaluate the function and prognostic value of m⁶A regulators. Further, we constructed nomogram to evaluate the prognosis of individual patient. The correlation between insulin-like growth factor 2 mRNA-binding protein 3 (IGF2BP3) and programmed cell death ligand 1 (PD-L1) was evaluated both in vitro and in vivo.

Results

We found that the tumor grade and DNA damage pathways were strongly correlated with distinct clusters. Furthermore, two risk score groups with six m⁶A regulators were identified using the least absolute shrinkage and selection operator (LASSO) and multivariate Cox regression analysis, which could be regarded as independent prognostic markers in patients with bladder cancer. The risk score pattern was linked to the tumor immune landscape, indicating a correlation between immune checkpoints and m⁶A regulators. Moreover, an m⁶A regulator, IGF2BP3, was found to be highly expressed in the tumor samples, regulating both the total and membrane-bound PD-L1 expression levels.

Conclusions

the results of this study revealed that the m⁶A clusters and patterns play crucial roles in the regulation of tumor immunity, which may be used to develop comprehensive treatment strategies for the management of bladder cancer.

Cancer Biology

Molecular Genetics

Oncology

m6A

immune landscape

bladder cancer

IGF2BP3

PD-L1

Of the 171 post-transcriptional modifications of RNA, including 5-methylcytosine (m⁵C), 7-methylguanosine (m⁷G), m¹G, m²G, and m⁶G, N¹-methyladenosine (m¹A) and N⁶-methyladenosine (m⁶A) have been identified in living organisms ¹. m⁶A is one of the most prominent and abundant internal modifications in eukaryotic mRNA and long noncoding RNA (lncRNA), accounting for 0.2–0.4% of the total adenosine residues and half of the total ribonucleotides in mammalian RNA ^2–4. Similar to dynamic and reversible epigenetic modifications of genomic DNA and proteins, m⁶A RNA modification is a reversible process in mammalian cells, which may be regulated by three vital factors, namely the methyltransferases, demethylases, and binding proteins, which are also known as writers, erasers, and readers, respectively. They can add, remove, or read an m⁶A site ^{5, 6}. Writing is the process of adding methylated modifications to RNA, including methyltransferase-like (METTL)-3, METTL14, METTL16, KIAA1429, WT1-associated protein (WTAP), RNA-binding motif protein (RBM)-15, RBM15B, Cbl proto-oncogene like 1 (CBLL1), and zinc finger CCCH-type containing 13 (ZC3H13) ^{5, 7, 8}, while the reversible process is mediated by erasers, including the fat mass and obesity-associated protein (FTO) and alkB homolog 5 (ALKBH5). Moreover, m⁶A indirectly affects RNA processing by recruiting specific reader proteins, including nuclear m⁶A readers, YTH domain containing (YTHDC)-1, heterogeneous nuclear ribonucleoprotein A2/B1 (HNRNPA2B1), heterogeneous nuclear ribonucleoprotein C (HNRNPC) etc., and cytoplasmic m⁶A readers, YTH N⁶-methyladenosine RNA binding protein (YTHDF)-1/2/3 and YTHDC2. m⁶A modifications are related to various biological functions, such as RNA splicing, translocation, stability, and translation ^{9, 10}, as well as multiple dysregulated biological processes, including aberrant proliferation, promotion of tumor metastasis, and inhibition of apoptosis ^11–14.

Tumor progression is attributed to multiple genetic and epigenetic variations in tumor cells ⁸. However, increasing evidence has shown that evading tumor surveillance is a hallmark of tumor development ¹⁵. The immunogenic interaction between the host tissues and the tumors, and the ability of the tumor to evade immune recognition could determine the prognosis of patients ¹⁶. Immunotherapy targeting immunological checkpoints, such as programmed cell death 1 (PD-1), programmed cell death 1 ligand 1 (PD-L1), and cytotoxic T lymphocyte-associated antigen-4 (CTLA-4), has been used as a potential therapeutic strategy for cancers. Unfortunately, the majority of patients receiving immunotherapeutic agents, anti-PD-1, anti-PD-L1, or anti-CTLA-4 antibodies gained no or minimal clinical benefits, different from the results of in vitro experiments ^17–19, suggesting that potentially undiscovered mechanisms play a role in regulating the immune system. As cancers are not solely neoplastic cells, they contain the tumor microenvironment (TME), which can be divided into immune and non-immune infiltrates, such as cytotoxic T cells, natural killer cells, dendritic cells (DCs), tumor-associated macrophages (TAMs), endothelial cells, and stromal cells ²⁰. The TME is a highly complex ecosystem, and various biological behaviors change through direct and indirect interactions with the TME components. As the understanding of complexity of TME has deepened, more evidence has shown that tumor-infiltrating cells play either tumor-suppressive or tumor-promoting roles, thus influencing cancer initiation and progression. For instance, M1 macrophages mainly produce pro-inflammatory cytokines that potentiate the anti-tumor immune response, while M2 macrophages promote ECM deposition and immunosuppression ²¹. Thus, it is important to clarify the immune infiltration at the tumor site and the biomarkers associated with TME, which might help to individually evaluate patients who could benefit from immunotherapies and will broaden our understanding of tumor immunity.

Bladder cancer (BCa) is more common in men than in women, and ranks 4th in incidence and 8th in mortality among the male population, according to the latest published cancer statistics ²². Approximately 75% of patients with newly diagnosed BCa are non-muscle invasive bladder cancer (NMIBC), and 25% have muscle invasive bladder cancer (MIBC) ²³. Bacillus Calmette-Guérin (BCG), which can activate human innate and adaptive immune responses, intravesical instillation is the current gold standard clinical treatment for NMIBC. Meanwhile, anti-PD-1/PD-L1 immunotherapy is the hotspot for advanced MIBC, and the response rate of immunotherapies is determined by various conditions, including tumor immunity and cancer cell immunogenicity ^{24, 25}. Recent studies have revealed an interaction between m⁶A modification and the TME. A study by Jiang et al. revealed that when co-cultured with M2 macrophages, the expression levels of ALKBH5 and the toll-like receptor (TLR)-4 increased in ovarian cancer cells, and TLR4 upregulated ALKBH5 expression and increased Nanog expression via mRNA demethylation ²⁶. Moreover, TNF-α inhibits the differentiation of mesenchymal stem cells by repressing FTO expression and FTO-mediated demethylation of Nanog mRNA levels and decreasing Nanog mRNA expression levels ²⁷. Zhang et al. integrated gastric cancer samples to establish m⁶A modification patterns and scoring systems, and found that TME characteristics were highly consistent with the patterns, suggesting that m⁶A modification played an insignificant role in tumor immunity in gastric cancer ⁸.

Several studies have demonstrated the function of m⁶A regulators in bladder cancer. For instance, METTL3 might act as an oncogene by interacting with the DiGeorge syndrome critical region 8 (DGCR8) and accelerating the pri-miR221/222 maturation to promote tumor proliferation ²⁸. Meanwhile, METTL3 plays a role in BCa progression by promoting the cancer cell growth and invasion by regulating a network that involves the AF4/FMR2 family member 4 (AFF4), nuclear factor-kappa B (NF-κB), and Myc ²⁹. Another mechanism of the m⁶A regulator METTL3/YTHDF2 mA axis directly degrades the mRNAs of SET domain containing 7 (SETD7) and Kruppel-like factor 4 (KLF4), contributing to the progression of BCa ³⁰. Analysis of the expression levels of METTL3 and CDCP1 in patients with BCa revealed that METTL3 and CDCP1 were strongly upregulated in the tumor samples, and the METTL3-CDCP1 axis could increase the tumor proliferation, migration, and invasion ³¹. However, the above studies only studied one or two m⁶A regulators and did not investigate the connection between m⁶A regulators and tumor immunity. Therefore, in this study, we comprehensively investigated the tumor immune landscape associated with m⁶A regulators, established a set of scoring patterns, and evaluated the prognostic value of this pattern for individual patients.

Dataset source acquisition

All data were obtained from The Cancer Genome Atlas (TCGA)-Urothelial Bladder Carcinoma (BLCA) database (https://portal.gdc.cancer.gov/) and the Genotype Tissue Expression project (GTEx) database (https://gtexportal.org/), including RNA sequencing data (fragments per kilobase of transcript per million fragments sequenced (FPKM) value) of gene expression, copy number variation (CNV), somatic mutation, and clinical information. The GTEx database includes over 10,000 bulk RNA-seq samples representing 53 different tissues (corresponding to 30 organs) obtained from 635 pre-healthy individuals, to link the influence of genetic variants on gene expression levels via quantitative trait loci analysis (eQTL) ³². The expression levels of normal samples from TCGA and GTEx databases were integrated and used for comparison with tumor samples. The Cancer Cell Line Encyclopedia (CCLE) database was used to evaluate the expression levels of the m⁶A regulator in several cell types, and the correlation with PD-L1 at the bladder cancer cell level. The protein levels of m⁶A regulators were determined using the Human Protein Atlas database (https://www.proteinatlas.org/). The Oncomine database was used to determine the insulin-like growth factor 2 mRNA-binding protein 3 (IGF2BP3) expression levels in several cancer types (https://www.oncomine.org/). The gene expression profiling interactive analysis (GEPIA) database was used to evaluate the prognostic value of IGF2BP3 (http://gepia.cancer-pku.cn/). The workflow of this study is shown in Supplementary Figure 1.

Selection of m⁶A regulators

A total of 24 m⁶A RNA methylation regulators were extracted from the database according to the relevant m⁶A studies ^{7, 8}, including nine writers (METTL3, METTL14, METTL16, RBM15, RBM15B, WTAP, KIAA1429, CBLL1, and ZC3H13), two erasers (ALKBH5 and FTO), and 13 readers (YTHDC1, YTHDC2, YTHDF1, YTHDF2, YTHDF3, IGF2BP1, IGF2BP2, IGF2BP3, HNRNPA2B1, HNRNPC, FMR1, leucine rich pentatricopeptide repeat containing (LRPPRC), and ELAV-like RNA binding protein 1 (ELAVL1)). The extracted information was used for further analyses.

Consensus clustering analysis

Unsupervised clustering of BCa samples was conducted to identify different m⁶A regulatory patterns based on their mRNA expression levels and classify the patients for further analysis. The “ConsensuClusterPlus” R package was used to perform the consensus clustering algorithm, which could determine the number of clusters and assess their stability ³³. Principal component analysis (PCA) was conducted to verify the grouping results.

Identification of differentially expressed genes (DEGs) and construction of the protein-protein interaction (PPI) network

The “limma” R package was applied to identify the differentially m⁶A expressed genes between normal and tumor samples ³⁴. The significance criteria for determining DEGs was set as the P value < 0.05, and |log2FC|>1. The Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) database and Cytoscape software were used to retrieve and construct a PPI network of the m⁶A regulator network.

Gene set variation analysis (GSVA)

To investigate the changes in pathway activity between different groups, GSVA was performed using the “GSVA” R packages. GSVA is a non-parametric, unsupervised method for estimating the variation in gene set enrichment and biological process activity of samples from an expression dataset ³⁵. In the present study, the gene sets of “c2.cp.kegg.v7.1.symbols” were downloaded from the Molecular Signatures Database (MsigDB) for running the GSVA analysis. We defined an adjusted P value < 0.05, and |log2FC| ≥0.08, as statistically significant.

Estimation of immune signatures, TME cell infiltration level and tumor purity in BCa

The single-sample gene-set enrichment analysis (ssGSEA) algorithm was used to quantify the relative abundance of immune signatures in the TME. The enrichment levels of 29 immune signatures were quantified for each sample. The gene sets representing each immune signature are shown in Supplementary Table 1, including B cells, Th2 cells, NK cells MHC class I cells, CD8+ T cells, and so on ³⁶. The “sparcl” R package was employed to divide the samples into three groups, including immunity high, medium and low for further analyses. ESTIMATE was used to evaluate TME cell infiltration level (including immune and stromal scores) and tumor purity for each sample ³⁷.

Cox regression analyses

Univariate and multivariate Cox regression analyses were used to assess the prognostic value of m⁶A regulators, and hazard ratios (HRs) >1 or <1 were regarded as risk and protective genes, respectively. The least absolute shrinkage and selection operator (LASSO) Cox regression algorithm was used to construct the optimal prognostic model out of the m⁶A regulators using the “glmnet” package in R. The LASSO analysis performed predictor selection, minimized over-fitting, selected genes to reduce bias, and developed the best survival-associated risk pattern ³⁸. After testing for collinearity, the sum of the Cox coefficient and gene values was calculated using the risk score based on the following formula:

Each patient was assigned a risk score based on integrative m⁶A regulator expression patterns. Tree-fold cross-validation and 1000 iterations were conducted to reduce the potential instability of the results.

The correlation between the expression levels of m⁶A regulators and mutation with immune cell infiltration

To explore the relationship between m⁶A regulators involved in the risk score pattern and the infiltrating immune cells, we utilized the Tumor IMmune Estimation Resource (TIMER) web tool (https://cistrome.shinyapps.io/timer/) to calculate the correlation coefficients between m⁶A regulator expression and mutations with infiltrated immune cells including B cells, CD4+ T cells, CD8+ T cells, macrophages, neutrophils, and dendritic cells ^{39, 40}.

Construction and validation of a predictive nomogram

Based on the results of the multivariable Cox regression model, a nomogram based on independent prognostic factors was constructed to predict 3- and 5-year OS. The nomogram provides a graphic representation linking individual patient factors to predict the survival probability of BCa patients ⁴¹. In addition, a bootstrapped resample with 1000 iterations was applied to verify the accuracy of the nomogram. Furthermore, the performance of the prognostic models was evaluated by receiver operating characteristic (ROC) analyses, and the concordance index (c-index) was measured to quantify the nomogram discrimination. A scale of 1.0, represents perfect predictions, and 0.5, the equivalent of a coin toss. The calibration of the model was assessed by comparing the observed survival with the predicted survival by plotting a calibration curve. The 45° line indicates a perfect calibration. Any deviation above or below the 45° line indicates underprediction or overprediction, respectively. Due to the limited conditions, no extra BCa cohort could be used as an externally validated database to evaluate the efficacy of model validation and prediction. Thus, only internal validation was conducted to evaluate the nomogram model.

Tissue specimens

This study was approved by the Medical Ethics Committee of the Shandong University School of Clinical Medicine. Twenty human BCa tissues were collected at the Qilu Hospital of Shandong University. Informed consent was obtained from all patients.

Cell cultures and manipulation

T24, 5637, and UMUC3 cell lines were purchased from the American Type Culture Collection (ATCC). T24 and 5637 cells were cultured in the Roswell Park Memorial Institute (RPMI)-1640 medium (Gibco; 11875093). UMUC3 cells were cultured in Dulbecco’s modified Eagle’s medium (DMEM; Gibco; 11995065). All media were supplemented with 10% fetal bovine serum (Gibco; 10099-141C). The cells were incubated at 37 °C in a humidified atmosphere with 5% CO₂.

Stable IGF2BP3 knockdown, overexpression cell lines, and their controls were generated as described previously ⁴². Lentiviruses were purchased from GeneChem Inc. (Shanghai, China).

Western blotting and antibodies

Western blotting was performed as previously described ⁴². The primary antibodies included anti-IGF2BP3 (Abcam; ab177477), anti-PD-L1 (Proteintech; 28076-1-AP), and anti-GAPDH (Proteintech; 60004-1-Ig).

RNA extraction and reverse transcription-polymerase chain reaction (RT-PCR)

Extraction of total RNA and RT-PCR were performed as previously described ⁴². Primers used were purchased from Sangon Biotech (Shanghai, China), and primer sequences are shown in Supplementary Table 2.

Flow cytometric assay

Flow cytometric assays were performed as previously described ⁴³. Briefly, single-cell suspensions were freshly prepared from the inidicated cells. Cells were washed once with PBS and stained with anti-PD-L1 antibody (Biolegend, 329705) for 30 min. The samples were analyzed on a FACSCanto II (BD Bioscience, USA) using FlowJo 7.6.5.

Statistical analysis

All analyses were conducted using the R v.4.0.0 and SPSS v.20.0 software. The “RCircos” R package was used to plot the CNV landscape of m⁶A regulators. The somatic mutation landscape was assessed using the “maftools” R package to plot the mutation summary, waterfall, and gene cloud figures. The correlations among different m⁶A regulators were computed by “corrplot” R package. “Survival” R package was adopted to analyze Kaplan-Meier curve analysis. The “forestplot” R package was conducted to visualize the univariate and multivariate prognostic analysis for risk score. The specificity and sensitivity of risk score were assessed using the ROC curve and the area under the curve (AUC), which were quantified by the “pROC” R package. Data from the two groups were evaluated using a two-tailed unpaired Student’s t-test. Categorical data were analyzed using the chi-square test. The correlation between continuous variables was assessed using the Spearman’s correlation analysis. Survival analysis was performed using log-rank test. Statistical significance was set at P < 0.05.

Landscape of somatic mutation and CNV mutation of m⁶A regulators in BCa

A total of 24 m⁶A regulators in BCa were used in the present study. We first clarified the incidence of CNV and somatic mutations in m⁶A regulators in BCa. Among the 412 samples, 116 (28.16%) experienced somatic mutations in m⁶A regulators (Figure 1A). METTL3 showed the highest mutation frequency among the m⁶A regulators (approximately 4%), while two writers (RBM15B and METTL16) and two readers (HNRNPC and FMR) did not exhibit any somatic mutations in BCa. Correlation analyses revealed that most m⁶A somatic mutations did not exhibit a co-occurrence relationship, except for FMR1 and YTHDF2, YTHDF1 and KIAA1429, WTAP and METTL3, ZC3H13, and RBM15 (Figure 1B). Next, we summarized the CNV mutation frequency among the m⁶A regulators, and KIAA1429, YTHDF1, YTHDF3, and IGF2BP2 had a widespread frequency of CNV amplification, while METTL16 and ALKBH5 showed high CNV deletion frequency (Figure 1C). We also explored the CNV mutation in normal tissues, and only 7 m⁶A regulators had a CNV mutation burden, with an extremely low frequency (Supplementary Figure 2). The location of the CNV mutation in m⁶A regulators on different chromosomes is shown in Figure 1D.

Profiles of mRNA and protein expression level of m⁶A methylation regulators in BCa

After exploring the mutation of m⁶A regulators, we investigated the mRNA expression levels of m⁶A regulators in normal bladder samples and tumor samples. The GTEx and TCGA datasets were merged for further analysis, with 28 normal samples and 411 tumor samples. As shown in Figures 2A, B, mRNA expression levels of CBLL1, ELAVL1, HNRNPA2B1, IGF2BP1, IGF2BP2, IGF2BP3, LRPPRC, RBM15, RBM15B, YTHDF1, YTHDF2, and YTHDF3 were significantly higher in tumor samples than in healthy samples, while the expression levels of FTO, METTL14, METTL16, WTAP, YTHDC1, YTHDC2, and ZC3H13 were decreased in tumor samples. Due to the functional similarity or complementation, the comprehensive landscape of m⁶A regulator connections was depicted by Spearman correlation analysis, STRING website, and Cytoscape software, METTL3 and YTHDF3 showed the strongest positive correlation, while METTL3 and IGF2BP2 showed the strongest negative relevance (Figures 2C–D). Not only did the m⁶A regulators with similar functions show a significant correlation, but a remarkable interaction was shown among writers, erasers, and readers. Moreover, correlations between writers and erasers were investigated to determine whether tumors with high eraser expression levels exhibited low writer expression levels. The results revealed that tumors with high expression of CBLL1 and METTL14 showed a high expression of FTO, while the high expression of CBLL1 and METTL14 showed low expression of ALKBH5. Tumors with high expression of ZC3H13 and WTAP showed high expression of FTO. However, ZC3H13 and WTAP did not interfere with ALKBH5 expression. The remaining writer genes did not affect the eraser genes ALKBH5 and FTO (Supplementary Figure 3A). Immunohistochemistry staining images of m⁶A regulator proteins were retrieved from the Human Protein Atlas, revealing cellular sublocalization and intensity (Figure 2E and Supplementary 3B). The above results revealed that cross-talk among m⁶A regulators might construct important modification patterns.

Identification of m⁶A regulators in two subgroups using consensus clustering

Further, the “ConsensusClusterPlus” R package was used to classify the subgroup number according to the m⁶A regulator expression levels. “k” was used to represent the number of subgroups. The empirical cumulative distribution function was plotted to analyze the optimal k value at which the cluster model achieved maximum stability (Figure 3A and Supplementary Figure 4A). The results showed that k=2 gained the most powerful clustering efficacy, and the samples were divided into two subclusters using unsupervised clustering (Figure 3B and Supplementary Figures 4B–F). PCA analysis was used to judge the classification, and the two clusters could gather together (Figure 3C). Prognostic analysis for the two clusters did not show a statistically significant difference, but a trend in overall survival (OS) (Figure 3D). We explored PCA analysis and prognostic analysis for k=3 and 4, and no significant benefits in OS were found (data not shown). However, the correlation analysis showed that clustering was associated with the tumor grade (Figure 3E). The m⁶A expression profiles showed that all m⁶A regulators were upregulated in cluster 2, with the exception of IGF2BP1 (Supplementary Figure 4G). GSVA enrichment analysis was performed to explore the biological behaviors of the two clusters. As shown in Figures 3F–G, cluster 1 presented enrichment pathways associated with metabolism, such as linoleic acid, arachidonic acid, retinol, drug, and xenobiotic metabolism. Cluster 2 was remarkably enriched in DNA damage, including mismatch repair, DNA replication, cell cycle, nucleotide excision repair, and spliceosome. Taken together, the results reveal that clusters based on m⁶A regulator expression could provide a few fresh outlooks for further study.

Characteristics of risk score pattern based on m⁶A regulators

To further explore the prognostic value of m⁶A regulators in BCa, we first performed a univariate Cox regression analysis on the expression levels of m⁶A regulators. The results showed that high expression of IGF2BP3 (hazard ratio [HR]=1.163, 95% confidence interval [CI]=1.059-1.278) and LRPPRC (HR=1.029, 95% CI=1.005-1.053) had worse survival outcomes in patients with BCa, whereas YTHDC1 (HR=0.929, 95% CI=0.885-0.975) and WTAP (HR=0.964, 95% CI=0.931-0.999) were regarded as protective markers for BCa (Figure 4A).

LASSO Cox regression analysis was performed to determine the optimal genes for selecting predictors and building the most regularized and parsimonious risk score pattern, and we only chose the prognostic value of m⁶A regulators p<0.1 for further analysis. The genes and their coefficients obtained from the LASSO analysis were used to calculate the risk scores for individual patients (Figure 4B, Supplementary Figure 5A, and Supplementary Table 3). The final LASSO model with the best optimal lambda included six m⁶A regulators (IGF2BP3, LRPPRC, YTHDC1, YTHDF2, and WTAP). To investigate the prognostic value of the risk score pattern, BCa patients were divided into low-and high-risk groups, and the Kaplan-Meier curve revealed that the patients in the high-risk group had a worse survival than patients in the low-risk group (Figure 4C). As shown in Figures 4D–E, green (low risk or alive) and red (high risk or dead) dots demonstrated significant differences between the low-and high-risk groups. ROC analysis showed a solid risk score pattern with AUC=0.64, indicating that the risk score could predict the OS of patients with BCa (Figure 4F). Next, we explored the correlation between the risk score pattern and clinical characteristics, and the risk score pattern was related to tumor grade, tumor stage, T status, M status, and N status (Figure 4G and Supplementary Table 4). Moreover, the expression of risk genes (IGF2BP3, LRPPRC, and FTO) was higher in patients at high risk, while YTHDC1, YTHDF2, and WTAP tended to be expressed in the low-risk group. Univariate and multivariate Cox analyses were performed to determine the independent prognostic value of the risk score pattern, and patient age, tumor T status, and risk score were independent prognostic predictors in patients with BCa (Figures 4H–I).

To better understand the function of the risk score pattern, we analyzed the GO analysis of DEGs based on expressions in low-and high-risk score groups. GO analysis indicated that upregulated genes in the high-risk group were enriched in malignancy-related biological processes, including extracellular matrix organization, extracellular structure organization, and antimicrobial humoral response, and downregulated genes were enriched in hormone metabolic processes and terpenoid metabolic processes (Supplementary Figures 5B–C). GSVA enrichment analysis was conducted to explore the different pathways between the two groups. The results revealed that the high-risk group was significantly related to the malignant pathways, including gap junction, focal adhesion, and ECM receptor interaction (Supplementary Figures 5D–E). Then, we investigate the distribution differences of somatic mutation between low and high risk score using “maftools” R package. As shown in Supplementary Figures 5F–K, the high-risk score group exhibited more somatic mutation burden than the low-risk score group, and missense mutation was the most common variant classification; the most common variant type was SNP, and C>T transversion was the most common type of SNV class. Moreover, the top three mutated genes were TP53, TTN, and KDM6A in the low-risk group and TP53, TTN, and KMT2D in the high-risk group, respectively. Taken together, the risk score pattern based on m⁶A regulators could be regarded as an independent prognostic factor in patients with BCa, and the high-risk group gained more malignant behaviors and more mutation burden.

Characteristics of immune landscape with risk score pattern

To explore the potential relationship between immunity and risk score pattern, we first divided samples into three clusters, immunity low, median and high using the ssGSEA score to quantify the immune cell types, functions and pathways, and the differences in 29 immune-associated gene sets were shown in three distinct immunity clusters (Supplementary Figures 6A–B). Next, we investigated the correlation between the immune landscape and the risk score pattern. As shown in Figure 5A, the enrichment of the immune landscape in the high-risk group was higher than that in the low-risk group. Moreover, the percentage of low immunity samples in the high-risk group was significantly lower than that in the low-risk group and more median immunity samples in the high-risk group than in the low-risk group (Figure 5B). In addition, comparing the stromal score, immune score, tumor purity, and ESTIMATE score between the two distinct risk score groups, we found that the high-risk group had significantly higher stromal scores, immune scores, and ESTIMATE scores, and lower tumor purity (Figures 5C–F). Taken together, these results suggest that the risk score pattern has a strong relationship with the immune landscape, and the potential mechanisms of m⁶A regulators in tumorigenesis and progression may be associated with tumor immunity.

Next, we explored the immune characteristics of independent m⁶A regulators in the risk group using the TIMER database to investigate the correlation between m⁶A regulator expression and immune cells, including B cells, CD8+ T cells, CD4+ T cells, macrophages, neutrophils, dendritic cells, and tumor purity. As shown in Figure 5G and Supplementary Figure 6C The expression of IGF2BP3 was positively correlated with macrophages, neutrophils, and dendritic cells and negatively correlated with tumor purity. The expression of LRPPRC was positively correlated with B cells, CD8+ T cells, neutrophils, and dendritic cells, and negatively correlated with CD4+ T cells. As for FTO, B cells, CD8+ T cells, macrophages, neutrophils, and dendritic cells were identified as significant co-expression cells. The expression of WTAP was negatively correlated with tumor purity, but positively correlated with CD8+ T cells, CD4+ T cells, neutrophils, and dendritic cells. The expression of YTHDC1 was only correlated with tumor purity, B cells, and macrophages. As for YTHDF2, tumor purity, B cells, CD8+ T cells, and neutrophils showed a strong correlation. The SCNA module, which was defined by GISTIC 2.0, was conducted to provide a comparison of immune infiltration levels in BCa with different somatic copy number alterations for m⁶A regulators. As shown in Figure 5H and Supplementary Figure 6D, IGF2BP3 amplification was associated with dendritic cells, and YTHDC1 deletion was related to B cells, amplification related to CD8+ T cells, neutrophils, and dendritic cells. Moreover, FTO amplification had a connection with B cells and macrophages, and deletion had a connection with CD8+ T cells. Interestingly, LRPPRC deletion and amplification were both associated with CD4+ T cells, neutrophils, and dendritic cells. The YTHDF2 mutation was associated with immune cells, except for CD8+ T cells and macrophages. The WTAP mutation is only related to CD4+ T cells and neutrophils. Furthermore, we investigated the co-expression of m⁶A regulators in the risk score model and several immune checkpoints (Figure 5I). The results indicate that m⁶A regulators are correlated with most immune checkpoints, including PD-L1 (also known as CD274). In summary, these results strongly indicate that the risk score pattern based on m⁶A regulators is significantly correlated with the tumor immune landscape.

Construction and validation of nomogram

A nomogram was established based on the independent factors using a multivariate Cox regression model to predict OS in patients with BCa (Figure 6A). The AUCs of the nomogram for predicting the 3- and 5-year OS were 0.685 and 0.695, respectively (Figures 6B–C). The c-index of the nomograms for OS in the training set was 0.676. As shown in Figures 6D–E, calibration plots were generated to validate the similarities between the actual survival rate and the survival prediction by the nomogram, and the results demonstrated that the 3- and 5-year survival rates predicted by the nomogram closely corresponded with the actual survival rates in the training set.

Moreover, 30 percent of patients with BCa were selected in the internal validation set. The AUCs in the validation set for predicting the 3- and 5-year OS were 0.747 and 0.723, respectively (Supplementary Figures 7A–B). The c-index of the nomogram in the validation set was 0.688. The results of the calibration plot suggested that the predicted 3- and 5-year survival rates were consistent with the actual survival rate within an acceptable margin of error in patients with BCa (Supplementary Figures 7C–D).

Characteristics of IGF2BP3 expression in cancers

Because of the important role of IGF2BP3 in risk score patterns, we used TCGA, GTEx CCLE, and Oncomine datasets to further understand IGF2BP3 in normal and tumor tissues. As shown in Figure 7A, the expression of IGF2BP3 was higher in BCa, cholangiocarcinoma (CHOL), colon adenocarcinoma (COAD), esophageal carcinoma (ESCA), head and neck squamous cell carcinoma (HNSC), kidney chromophobe (KICH), kidney renal clear cell carcinoma (KIRC), kidney renal papillary cell carcinoma (KIRP), liver hepatocellular carcinoma (LIHC), lung adenocarcinoma (LUAD), lung squamous cell carcinoma(LUSC), stomach adenocarcinoma (STAD), uterine corpus endometrial carcinoma (UCEC), and low expression in thyroid carcinoma (THCA), compared to their corresponding normal tissues (Figure 7A). Moreover, the CCLE dataset was used to evaluate the expression levels of IGF2BP3 in various tumor cell lines. The results showed that the top three expression levels in tumor cell lines were liver cancer, lymphoma, and medulloblastoma. IGF2BP3 seemed to be positively associated with PD-L1 expression in BCa cell lines (Figures 7B–C).

The Oncomine database was used to determine the expression level of IGF2BP3. And as shown in Figure 7D, IGF2BP3 in most cancer types showed high expression levels, except for kidney cancer and myeloma, which were opposite to the TCGA database. Furthermore, the correlation between IGF2BP3 and PD-L1 in patients with BCa showed a trend similar to that of bladder cancer cell lines from CCLE (Figure 7E). Next, we investigated the prognostic value of IGF2BP3 in BCa using the GEPIA website, and the results revealed that patients with high expression of IGF2BP3 had worse prognosis in BCa (Figure 7F). The GETx database indicated that the expression level of IGF2BP3 in male and female bone marrow was significantly high. To compare the expression differences between males and females, and there was no difference in the expression of IGF2BP3 in most female and male tissues, except for blood vessels, brain, breast, and lung (Figures 7G–J). Taken together, these results reveal that the expression of IGF2BP3 is high in various tumors and is associated with PD-L1, which may be a potential target for anti-PD-L1 immunotherapy.

The correlation between IGF2BP3 and PD-L1 in BCa cells and tumor specimen

Given that IGF2BP3 had a strong correlation with PD-L1 analyzed using a public database, we examined the expression levels of IGF2BP3 and PD-L1 in vitro. Stable IGF2BP3 overexpression and knockdown of T24, 5637, and UMUC3 cells were established, and the results revealed that overexpression of IGF2BP3 significantly increased, while knockdown of IGF2BP3 decreased both the protein and mRNA levels of PD-L1 in BCa cells (Figures 8A–F). Further, flow cytometric assay showed that overexpression of IGF2BP3 significantly enhanced membrane-bound PD-L1 expression, and knockdown of IGF2BP3 decreased membrane-bound PD-L1 expression in T24 cells (Figures 8G–H). The correlation between IGF2BP3 and PD-L1 expression was analyzed using BCa specimens. As shown in Figure 8I, positively correlated expression between IGF2BP3 and PD-L1 was found in 14/20 (70%) tumor specimens. Taken together, these data demonstrate that IGF2BP3 regulates both total and membrane-bound PD-L1 expression levels in BCa.

Although several studies have explored the m⁶A regulators in tumorigenesis and tumor development, the comprehensive analysis of m⁶A regulators with tumor immune landscape in bladder cancer has been poorly investigated. Here, we reveal that distinct clusters and risk groups are associated with tumor immunity and have prognostic value for patients with BCa. Furthermore, PD-L1 was identified as a potential target of IGF2BP3, and IGF2BP3 can regulate both total and membrane-bound PD-L1 expression levels.

Modified RNA bases have been discovered for over six decades. After the first RNA demethylase (FTO) was identified, the research field was revived, and the formation of m⁶A is a reversible process ^{6, 44}. The m⁶A regulators can be divided into three functional groups: writers, erasers, and readers. Increasing evidence has shown that m⁶A regulators play a crucial role in various pathophysiological processes, including circadian rhythms, spermatogenesis, DNA damage response, tumorigenesis, and tumor progression ^45–48.

The TME consists of immune cells and non-immune cells, which could play a crucial role in tumor growth and progression, and tumor-infiltrated immune cells are highly associated with tumorigenesis, angiogenesis, and metastasis ⁴⁹. Meanwhile, the imbalance between tumor cell growth and elimination might activate immunosurveillance. Therefore, understanding the crosstalk between TME and tumor cells might be useful for assessing the prognosis and improving the response rate of immunotherapy for individual patients with various cancers.

Bladder is an organ highly correlated with immune, and immunotherapy, including anti-PD-L1 and anti-CTLA4 becomes a hotspot in advanced BCa treatment. The Food and Drug Agency (FDA) and the European Medicines Agency (EMA) granted accelerated approval to atezolizumab and pembrolizumab as first-line metastatic cisplatin-unfit BCa. However, anti-PD-L1 treatment showed limited efficacy in the first-line phase III clinical trials ₁₈, with a relatively low response rate of approximately 20% ⁵⁰. Moreover, several studies have recently demonstrated that PD-1 and PD-L1 expression are not reliable biomarkers for predicting the benefits of immunotherapy. ^{51, 52} A retrospective study demonstrated that patients progressing to frontline PD-1/PD-L1 immunotherapy were even at risk of early death, excluding them from experiencing potential benefit from subsequent systemic treatment ¹⁹.

Most studies have focused on m⁶A regulators or immunotherapy; however, the correlation between m⁶A regulators and tumor immunity has not been fully recognized, and only a few studies have demonstrated the potential relationship between m⁶A regulators and TME anti-tumor immune responses in various cell types, such as gastric cancer, melanoma, and dendritic cells ^{8, 53, 54}. Here, we first identified two distinct m⁶A clusters and constructed a risk score pattern based on m⁶A regulators to reveal the potential pathways and functional processes, predict the prognosis of patients with BCa, and investigate the correlation between m⁶A regulators and tumor immunity. Moreover, we analyzed one of the m⁶A regulators in pattern, IGF2BP3, and identified its expression level, prognostic value, and association with PD-L1. Clarifying the role of risk score pattern with TME will contribute to broadening the understanding of TME antitumor immune response and suggest appropriate effective immunotherapy strategies for individual BCa patients.

Although the prognostic value of the m⁶A cluster was limited, the m⁶A cluster was associated with tumor grade, and two clusters showed significantly distinct pathway enrichment. Cluster 1 was characterized by metabolism, and cluster 2 was characterized by DNA damage. Moreover, the correlation analysis revealed that the most significant positive and negative relevance were METTL3 with YTHDF3 and IGF2BP2, respectively, showing different functions in BCa. The function of m⁶A regulators primarily depends on reader proteins. YTHDF2 could induce mRNA degradation and YTHDF1 and YTHDF3 could initiate mRNA translation, while the IGF2BP family could enhance the stability of target mRNA ^55–58. A comprehensive analysis of m⁶A regulators revealed that the mRNA expression of METTL4 and YTHDF3 was higher in high-grade tumors than in low-grade tumors, and YTHDC1 was upregulated in the I/II stage, compared to the III/IV stage ⁵⁹.

Furthermore, the risk score pattern based on 6 m⁶A regulators revealed its prognostic value for OS in patients with BCa, and the risk score pattern was highly associated with pathological features, such as T status, M status, N status, clinical stage, and tumor grade. Moreover, in the present study, we found an association between risk score pattern and TME, and high immune score, high stromal score, high ESTIMATE score, and low tumor purity were found in the high-risk score group. We next found that the expression and mutation of individual m⁶A regulators in the risk score pattern was associated with immune cells and immune checkpoints, which could underlie part of the mechanism of the risk score pattern. A nomogram was constructed to evaluate the prognostic value of individual patients for predicting 3- and 5-year survival times. If the physicians were able to estimate whether individual patients had shorter or longer than the median OS according to their expression of m⁶A regulators in tumor tissues, it would be useful for patients with different treatment strategies ⁶⁰.

IGF2BP3 was highly expressed in the high-risk group, and a recent study demonstrated that IGF2BP3 could be regarded as an independent prognostic factor in NMIBC, which could present a subgroup of patients with high probability of relapse, progression, and metastasis ⁶¹. A comprehensive study has reported that the expression of IGF2BP3 was detected in 76 different normal tissue types and 3889 cancer samples from 95 different tumor categories, IGF2BP3 overexpression has been found in various cancer types, and IGF2BP3 is typically associated with aggressive tumor features ⁶². IGF2BP3 has been shown to directly interact with ULBP2 mRNA, thereby reducing ULBP2 surface expression. IGF2BP3 indirectly interacts with MICB. The IGF2BP3-mediated pathway leads to impaired NK cell recognition of transformed cells to facilitate tumor immune escape ⁶³. In the present study, we found a positive correlation between IGF2BP3 and PD-L1, and IGF2BP3 could regulate total and membrane-bound PD-L1 expression levels in BCa cells, which implies the potential role of IGF2BP3 in anti-PD-L1 immunotherapy.

This study has several limitations. First, because of the limited clinical database on BCa, only TCGA patients with clinical characteristics were included. Second, immunohistochemical staining of m⁶A regulators was obtained from the public database, and the protein levels of m⁶A regulators will be explored in further studies. Third, our nomogram only underwent internal validation; it could be more powerful to obtain an external validation with a large multicenter cohort. Despite considering the limitations of the present study, our findings provide novel insights for m⁶A regulator clusters, risk score based on m⁶A regulators, and identified the association between tumor immunity and m⁶A regulators.

In summary, the present study investigated the cluster and prognosis of m⁶A regulators in BCa and found that the expression of m⁶A regulators is highly correlated with clinicopathological characteristics. We also constructed a risk score pattern and nomogram to evaluate the OS of patients with BCa. Moreover, we illustrated the relationship between m⁶A regulators and the TME. Therefore, our study provides important ideas for improving the clinical outcomes of patients with BCa, which may be used to develop different immunotherapies based on the expression levels of m⁶A regulators.

Acknowledgements

Not applicable.

Author Contributions

Conception and design: BS, JC, XL conceived the study; JC, XJ, YX performed data analysis; WW, GZ designed the figures and tables; SC, JC, XL wrote the manuscript; BS and SC reviewed and supervised the manuscript. Final approval of manuscript was gained form all authors, and all authors agree to be accountable for the content of the work.

Funding

This work was supported by the National Natural Science Foundation of China (Grant No. 81900637 to BS, Grant No. 81800672 to SC), the Tai Shan Scholar Foundation (ts201511092 to BS) and Primary Research & Development Plan of Shandong Province (2019GSF108123 to SC). And the authors have declared that no conflicts of interest exists.

Availability of data and materials

The data and materials in the current study are available from the corresponding author on reasonable request.

Ethics approval and consent to participate

The experimental protocol was established, according to the ethical guidelines of the Helsinki Declaration and was approved by the Medical Ethics Committee of the Shandong University School of Clinical Medicine. Informed consent was obtained from all patients.

Consent for publication

Not applicable.

Competing interests

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.

Boccaletto P, Machnicka MA, Purta E, Piatkowski P, Baginski B, Wirecki TK, de Crécy-Lagard V, Ross R, Limbach PA, Kotter A, et al. MODOMICS: a database of RNA modification pathways. 2017 update. Nucleic Acids Res. 2018;46:D303-d307. doi:10.1093/nar/gkx1030.
Desrosiers R, Friderici K, Rottman F. Identification of methylated nucleosides in messenger RNA from Novikoff hepatoma cells. Proc Natl Acad Sci U S A. 1974;71:3971–5. doi:10.1073/pnas.71.10.3971.
Alarcón CR, Lee H, Goodarzi H, Halberg N, Tavazoie SF. N6-methyladenosine marks primary microRNAs for processing. Nature. 2015;519:482–5. doi:10.1038/nature14281.
Zhao BS, Roundtree IA, He C. Post-transcriptional gene regulation by mRNA modifications. Nat Rev Mol Cell Biol. 2017;18:31–42. doi:10.1038/nrm.2016.132.
Yang Y, Hsu PJ, Chen YS, Yang YG. Dynamic transcriptomic m(6)A decoration: writers, erasers, readers and functions in RNA metabolism. Cell Res. 2018;28:616–24. doi:10.1038/s41422-018-0040-8.
Jia G, Fu Y, Zhao X, Dai Q, Zheng G, Yang Y, Yi C, Lindahl T, Pan T, Yang YG, et al. N6-methyladenosine in nuclear RNA is a major substrate of the obesity-associated FTO. Nat Chem Biol. 2011;7:885–7. doi:10.1038/nchembio.687.
Lan Q, Liu PY, Haase J, Bell JL, Hüttelmaier S, Liu T. The Critical Role of RNA m(6)A Methylation in Cancer. Cancer Res. 2019;79:1285–92. doi:10.1158/0008-5472.Can-18-2965.
Zhang B, Wu Q, Li B, Wang D, Wang L, Zhou YL. m(6)A regulator-mediated methylation modification patterns and tumor microenvironment infiltration characterization in gastric cancer. Mol Cancer. 2020;19:53. doi:10.1186/s12943-020-01170-0.
Liu Y, Guo X, Zhao M, Ao H, Leng X, Liu M, Wu C, Ma J, Zhu J. Contributions and prognostic values of m(6) A RNA methylation regulators in non-small-cell lung cancer. J Cell Physiol. 2020;235:6043–57. doi:10.1002/jcp.29531.
Wang S, Sun C, Li J, Zhang E, Ma Z, Xu W, Li H, Qiu M, Xu Y, Xia W, et al. Roles of RNA methylation by means of N(6)-methyladenosine (m(6)A) in human cancers. Cancer Lett. 2017;408:112–20. doi:10.1016/j.canlet.2017.08.030.
Visvanathan A, Patil V, Arora A, Hegde AS, Arivazhagan A, Santosh V, Somasundaram K. Essential role of METTL3-mediated m(6)A modification in glioma stem-like cells maintenance and radioresistance. Oncogene. 2018;37:522–33. doi:10.1038/onc.2017.351.
Wang X, Feng J, Xue Y, Guan Z, Zhang D, Liu Z, Gong Z, Wang Q, Huang J, Tang C, et al. Structural basis of N(6)-adenosine methylation by the METTL3-METTL14 complex. Nature. 2016;534:575–8. doi:10.1038/nature18298.
Bansal H, Yihua Q, Iyer SP, Ganapathy S, Proia DA, Penalva LO, Uren PJ, Suresh U, Carew JS, Karnad AB, et al. WTAP is a novel oncogenic protein in acute myeloid leukemia. Leukemia. 2014;28:1171–4. doi:10.1038/leu.2014.16.
Su R, Dong L, Li C, Nachtergaele S, Wunderlich M, Qing Y, Deng X, Wang Y, Weng X, Hu C, et al. R-2HG Exhibits Anti-tumor Activity by Targeting FTO/m(6)A/MYC/CEBPA Signaling. Cell. 2018;172:90–105.e23. doi:10.1016/j.cell.2017.11.031.
Hanahan D, Weinberg RA. Hallmarks of cancer: the next generation. Cell. 2011;144:646–74. doi:10.1016/j.cell.2011.02.013.
Hoenicke L, Zender L. Immune surveillance of senescent cells–biological significance in cancer- and non-cancer pathologies. Carcinogenesis. 2012;33:1123–6. doi:10.1093/carcin/bgs124.
Topalian SL, Hodi FS, Brahmer JR, Gettinger SN, Smith DC, McDermott DF, Powderly JD, Carvajal RD, Sosman JA, Atkins MB, et al. Safety, activity, and immune correlates of anti-PD-1 antibody in cancer. N Engl J Med. 2012;366:2443–54. doi:10.1056/NEJMoa1200690.
Gourd E. EMA restricts use of anti-PD-1 drugs for bladder cancer. Lancet Oncol. 2018;19:e341. i: 10.1016/s1470-2045(18)30433-9.. ; ; do.
Gómez de Liaño Lista A, van Dijk N, de V Oria de Rueda, Necchi G, Lavaud A, Morales-Barrera P, Alonso Gordoa R, Maroto T, Ravaud P, Durán A I, et al. Clinical outcome after progressing to frontline and second-line Anti-PD-1/PD-L1 in advanced urothelial cancer. Eur Urol. 2020;77:269–276; doi: 10.1016/j.eururo.2019.10.004.
Giraldo NA, Sanchez-Salas R, Peske JD, Vano Y, Becht E, Petitprez F, Validire P, Ingels A, Cathelineau X, Fridman WH, et al. The clinical role of the TME in solid cancer. Br J Cancer. 2019;120:45–53. doi:10.1038/s41416-018-0327-z.
Mantovani A, Marchesi F, Malesci A, Laghi L, Allavena P. Tumour-associated macrophages as treatment targets in oncology. Nat Rev Clin Oncol. 2017;14:399–416. doi:10.1038/nrclinonc.2016.217.
Siegel RL, Miller KD, Jemal A. Cancer statistics. 2020. CA Cancer J Clin. 2020;70:7–30; doi: 10.3322/caac.21590.
Sanli O, Dobruch J, Knowles MA, Burger M, Alemozaffar M, Nielsen ME, Lotan Y. Bladder cancer. Nat Rev Dis Primers. 2017;3:17022. doi:10.1038/nrdp.2017.22.
van Dijk N, Funt SA, Blank CU, Powles T, Rosenberg JE, van der Heijden MS. The Cancer Immunogram as a Framework for Personalized Immunotherapy in Urothelial Cancer. Eur Urol. 2019;75:435–44. doi:10.1016/j.eururo.2018.09.022.
Kamat AM, Li R, O'Donnell MA, Black PC, Roupret M, Catto JW, Comperat E, Ingersoll MA, Witjes WP, McConkey DJ, et al. Predicting Response to Intravesical Bacillus Calmette-Guérin Immunotherapy: Are We There Yet? A Systematic Review. Eur Urol. 2018;73:738–48. doi:10.1016/j.eururo.2017.10.003.
Jiang Y, Wan Y, Gong M, Zhou S, Qiu J, Cheng W. RNA demethylase ALKBH5 promotes ovarian carcinogenesis in a simulated tumour microenvironment through stimulating NF-κB pathway. J Cell Mol Med. 2020. doi:10.1111/jcmm.15228.
Wang Y, Wang R, Yao B, Hu T, Li Z, Liu Y, Cui X, Cheng L, Song W, Huang S, et al. TNF-α suppresses sweat gland differentiation of MSCs by reducing FTO-mediated m(6)A-demethylation of Nanog mRNA. Sci China Life Sci. 2020;63:80–91. doi:10.1007/s11427-019-9826-7.
Han J, Wang JZ, Yang X, Yu H, Zhou R, Lu HC, Yuan WB, Lu JC, Zhou ZJ, Lu Q, et al. METTL3 promote tumor proliferation of bladder cancer by accelerating pri-miR221/222 maturation in m6A-dependent manner. Mol Cancer. 2019;18:110. doi:10.1186/s12943-019-1036-9.
Cheng M, Sheng L, Gao Q, Xiong Q, Zhang H, Wu M, Liang Y, Zhu F, Zhang Y, Zhang X, et al. The m(6)A methyltransferase METTL3 promotes bladder cancer progression via AFF4/NF-κB/MYC signaling network. Oncogene. 2019;38:3667–80. doi:10.1038/s41388-019-0683-z.
Xie H, Li J, Ying Y, Yan H, Jin K, Ma X, He L, Xu X, Liu B, Wang X, et al. METTL3/YTHDF2 m(6) A axis promotes tumorigenesis by degrading SETD7 and KLF4 mRNAs in bladder cancer. J Cell Mol Med. 2020;24:4092–104. doi:10.1111/jcmm.15063.
Yang F, Jin H, Que B, Chao Y, Zhang H, Ying X, Zhou Z, Yuan Z, Su J, Wu B, et al. Dynamic m(6)A mRNA methylation reveals the role of METTL3-m(6)A-CDCP1 signaling axis in chemical carcinogenesis. Oncogene. 2019;38:4755–72. doi:10.1038/s41388-019-0755-0.
Battle A, Brown CD, Engelhardt BE, Montgomery SB. Genetic effects on gene expression across human tissues. Nature. 2017;550:204–13. doi:10.1038/nature24277.
Wilkerson MD, Hayes DN. ConsensusClusterPlus: a class discovery tool with confidence assessments and item tracking. Bioinformatics. 2010;26:1572–3. doi:10.1093/bioinformatics/btq170.
Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, Smyth GK. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 2015;43:e47. doi:10.1093/nar/gkv007.
Hänzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinformatics. 2013;14:7. doi:10.1186/1471-2105-14-7.
He Y, Jiang Z, Chen C, Wang X. Classification of triple-negative breast cancers based on Immunogenomic profiling. J Exp Clin Cancer Res. 2018;37:327. doi:10.1186/s13046-018-1002-1.
Yoshihara K, Shahmoradgoli M, Martínez E, Vegesna R, Kim H, Torres-Garcia W, Treviño V, Shen H, Laird PW, Levine DA, et al. Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun. 2013;4:2612. doi:10.1038/ncomms3612.
Li YM, Li ZL, Chen F, Liu Q, Peng Y, Chen M. A LASSO-derived risk model for long-term mortality in Chinese patients with acute coronary syndrome. J Transl Med. 2020;18:157. doi:10.1186/s12967-020-02319-7.
Li T, Fan J, Wang B, Traugh N, Chen Q, Liu JS, Li B, Liu XS. TIMER: A Web Server for Comprehensive Analysis of Tumor-Infiltrating Immune Cells. Cancer Res. 2017;77:e108–10. doi:10.1158/0008-5472.Can-17-0307.
Li B, Severson E, Pignon JC, Zhao H, Li T, Novak J, Jiang P, Shen H, Aster JC, Rodig S, et al. Comprehensive analyses of tumor immunity: implications for cancer immunotherapy. Genome Biol. 2016;17:174. doi:10.1186/s13059-016-1028-7.
Cui J, Chen S, Bo Q, Wang S, Zhang N, Yu M, Wang W, Han J, Zhu Y, Shi B. Preoperative prognostic nutritional index and nomogram predicting recurrence-free survival in patients with primary non-muscle-invasive bladder cancer without carcinoma in situ. Onco Targets Ther. 2017;10:5541–50. doi:10.2147/ott.S146990.
Liu X, Cui J, Gong L, Tian F, Shen Y, Chen L, Wang Y, Xia Y, Liu L, Ye X, et al. The CUL4B-miR-372/373-PIK3CA-AKT axis regulates metastasis in bladder cancer. Oncogene. 2020;39:3588–603. doi:10.1038/s41388-020-1236-1.
Song Y, Li P, Qin L, Xu Z, Jiang B, Ma C, Shao C, Gong Y. CUL4B negatively regulates Toll-like receptor-triggered proinflammatory responses by repressing Pten transcription. Cell Mol Immunol. 2019. doi:10.1038/s41423-019-0323-0.
Jia G, Fu Y, He C. Reversible RNA adenosine methylation in biological regulation. Trends Genet. 2013;29:108–15. doi:10.1016/j.tig.2012.11.003.
Xiang Y, Laurent B, Hsu CH, Nachtergaele S, Lu Z, Sheng W, Xu C, Chen H, Ouyang J, Wang S, et al. RNA m(6)A methylation regulates the ultraviolet-induced DNA damage response. Nature. 2017;543:573–6. doi:10.1038/nature21671.
Lin Z, Hsu PJ, Xing X, Fang J, Lu Z, Zou Q, Zhang KJ, Zhang X, Zhou Y, Zhang T, et al. Mettl3-/Mettl14-mediated mRNA N(6)-methyladenosine modulates murine spermatogenesis. Cell Res. 2017;27:1216–30. doi:10.1038/cr.2017.117.
Fustin JM, Doi M, Yamaguchi Y, Hida H, Nishimura S, Yoshida M, Isagawa T, Morioka MS, Kakeya H, Manabe I, et al. RNA-methylation-dependent RNA processing controls the speed of the circadian clock. Cell. 2013;155:793–806. doi:10.1016/j.cell.2013.10.026.
Pan Y, Ma P, Liu Y, Li W, Shu Y. Multiple functions of m(6)A RNA methylation in cancer. J Hematol Oncol. 2018;11:48. doi:10.1186/s13045-018-0590-8.
Taddei ML, Giannoni E, Comito G, Chiarugi P. Microenvironment and tumor cell plasticity: an easy way out. Cancer Lett. 2013;341:80–96. doi:10.1016/j.canlet.2013.01.042.
Szabados B, van Dijk N, Tang YZ, van der Heijden MS, Wimalasingham A, Gomez de Liano A, Chowdhury S, Hughes S, Rudman S, Linch M, et al. Response Rate to Chemotherapy After Immune Checkpoint Inhibition in Metastatic Urothelial Cancer. Eur Urol. 2018;73:149–52. doi:10.1016/j.eururo.2017.08.022.
Fuchs CS, Doi T, Jang RW, Muro K, Satoh T, Machado M, Sun W, Jalal SI, Shah MA, Metges JP, et al. Safety and Efficacy of Pembrolizumab Monotherapy in Patients With Previously Treated Advanced Gastric and Gastroesophageal Junction Cancer: Phase 2 Clinical KEYNOTE-059 Trial. JAMA Oncol. 2018;4:e180013. doi:10.1001/jamaoncol.2018.0013.
Roh W, Chen PL, Reuben A, Spencer CN, Prieto PA, Miller JP, Gopalakrishnan V, Wang F, Cooper ZA, Reddy SM, et al. Integrated molecular analysis of tumor biopsies on sequential CTLA-4 and PD-1 blockade reveals markers of response and resistance. Sci Transl Med. 2017;9; doi:10.1126/scitranslmed.aah3560.
Yang S, Wei J, Cui YH, Park G, Shah P, Deng Y, Aplin AE, Lu Z, Hwang S, He C, et al. m(6)A mRNA demethylase FTO regulates melanoma tumorigenicity and response to anti-PD-1 blockade. Nat Commun. 2019;10:2782. doi:10.1038/s41467-019-10669-0.
Han D, Liu J, Chen C, Dong L, Liu Y, Chang R, Huang X, Liu Y, Wang J, Dougherty U, et al. Anti-tumour immunity controlled through mRNA m(6)A methylation and YTHDF1 in dendritic cells. Nature. 2019;566:270–4. doi:10.1038/s41586-019-0916-x.
Frye M, Harada BT, Behm M, He C. RNA modifications modulate gene expression during development. Science. 2018;361:1346–9. doi:10.1126/science.aau1646.
Wang X, Zhao BS, Roundtree IA, Lu Z, Han D, Ma H, Weng X, Chen K, Shi H, He C. N(6)-methyladenosine Modulates Messenger RNA Translation Efficiency. Cell. 2015;161:1388–99. doi:10.1016/j.cell.2015.05.014.
Li A, Chen YS, Ping XL, Yang X, Xiao W, Yang Y, Sun HY, Zhu Q, Baidya P, Wang X, et al. Cytoplasmic m(6)A reader YTHDF3 promotes mRNA translation. Cell Res. 2017;27:444–7. doi:10.1038/cr.2017.10.
Huang H, Weng H, Sun W, Qin X, Shi H, Wu H, Zhao BS, Mesquita A, Liu C, Yuan CL, et al. Recognition of RNA N(6)-methyladenosine by IGF2BP proteins enhances mRNA stability and translation. Nat Cell Biol. 2018;20:285–95. doi:10.1038/s41556-018-0045-z.
Lobo J, Barros-Silva D, Henrique R, Jerónimo C. The Emerging Role of Epitranscriptomics in Cancer: Focus on Urological Tumors. Genes (Basel). 2018;9; doi:10.3390/genes9110552.
Kim SY, Yoon MJ, Park YI, Kim MJ, Nam BH, Park SR. Nomograms predicting survival of patients with unresectable or metastatic gastric cancer who receive combination cytotoxic chemotherapy as first-line treatment. Gastric Cancer. 2018;21:453–63. doi:10.1007/s10120-017-0756-z.
Yang F, Zhou Q, Meng L, Xing N. IMP3 is a biomarker for non-muscle-invasive urothelial carcinoma of the bladder associated with an aggressive phenotype. Medicine. 2019;98:e16009. doi:10.1097/md.0000000000016009.
Burdelski C, Jakani-Karimi N, Jacobsen F, Möller-Koop C, Minner S, Simon R, Sauter G, Steurer S, Clauditz TS, Wilczak W. IMP3 overexpression occurs in various important cancer types and is linked to aggressive tumor features: A tissue microarray study on 8,877 human cancers and normal tissues. Oncol Rep. 2018;39:3–12. doi:10.3892/or.2017.6072.
Schmiedel D, Tai J, Yamin R, Berhani O, Bauman Y, Mandelboim O. The RNA binding protein IMP3 facilitates tumor immune escape by downregulating the stress-induced ligands ULPB2 and MICB. Elife 2016;5; doi:10.7554/eLife.13426.

Supplementaldata.pdf

Download PDF

Journal Publication

published 11 Feb, 2022

Read the published version in Cancer Cell International →

Review #2 received at journal
14 Nov, 2021
Reviewer #2 agreed at journal
27 Oct, 2021
Review #1 received at journal
14 Sep, 2021
Reviews received at journal
03 Sep, 2021
Reviewers invited by journal
03 Sep, 2021
Editor invited by journal
03 Sep, 2021
Reviewer #1 agreed at journal
02 Sep, 2021
Editor assigned by journal
29 Aug, 2021
Submission checks completed at journal
29 Aug, 2021
First submitted to journal
28 Aug, 2021

You are reading this latest preprint version

Comprehensive analysis of N6-methyladenosine regulators with the tumor immune landscape and correlation between the insulin-like growth factor 2 mRNA-binding protein 3 and programmed death ligand 1 in bladder cancer

Status:

Journal Publication

Version 1

Abstract

Background

Methods

Results

Conclusions

Figures

Background

Methods

Dataset source acquisition

Selection of m6A regulators

Consensus clustering analysis

Identification of differentially expressed genes (DEGs) and construction of the protein-protein interaction (PPI) network

Gene set variation analysis (GSVA)

Estimation of immune signatures, TME cell infiltration level and tumor purity in BCa

Cox regression analyses

The correlation between the expression levels of m6A regulators and mutation with immune cell infiltration

Construction and validation of a predictive nomogram

Tissue specimens

Cell cultures and manipulation

Western blotting and antibodies

RNA extraction and reverse transcription-polymerase chain reaction (RT-PCR)

Flow cytometric assay

Statistical analysis

Results

Landscape of somatic mutation and CNV mutation of m6A regulators in BCa

Profiles of mRNA and protein expression level of m6A methylation regulators in BCa

Identification of m6A regulators in two subgroups using consensus clustering

Characteristics of risk score pattern based on m6A regulators

Characteristics of immune landscape with risk score pattern

Construction and validation of nomogram

Characteristics of IGF2BP3 expression in cancers

The correlation between IGF2BP3 and PD-L1 in BCa cells and tumor specimen

Discussion

Declarations

References

Supplementary Files

Status:

Journal Publication

Version 1

Comprehensive analysis of N⁶-methyladenosine regulators with the tumor immune landscape and correlation between the insulin-like growth factor 2 mRNA-binding protein 3 and programmed death ligand 1 in bladder cancer

Selection of m⁶A regulators

The correlation between the expression levels of m⁶A regulators and mutation with immune cell infiltration

Landscape of somatic mutation and CNV mutation of m⁶A regulators in BCa

Profiles of mRNA and protein expression level of m⁶A methylation regulators in BCa

Identification of m⁶A regulators in two subgroups using consensus clustering

Characteristics of risk score pattern based on m⁶A regulators