Landscape of somatic mutation and CNV mutation of m6A regulators in BCa
A total of 24 m6A regulators in BCa were used in the present study. We first clarified the incidence of CNV and somatic mutations in m6A regulators in BCa. Among the 412 samples, 116 (28.2%) experienced somatic mutations in m6A regulators (Fig. 1A). METTL3 showed the highest mutation frequency among the m6A 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 m6A somatic mutations did not exhibit a co-occurrence relationship, except for FMR1 and YTHDF2, YTHDF1 and KIAA1429, WTAP and METTL3, ZC3H13, and RBM15 (Fig. 1B). Next, we summarized the CNV mutation frequency among the m6A regulators, and KIAA1429, YTHDF1, YTHDF3, and IGF2BP2 had a widespread frequency of CNV amplification, while METTL16 and ALKBH5 showed high CNV deletion frequency (Fig. 1C). We also explored the CNV mutation in normal tissues, and only 7 m6A regulators had a CNV mutation burden, with an extremely low frequency (Additional file 1: Fig. S2). The location of the CNV mutation in m6A regulators on different chromosomes is shown in Fig. 1D.
Profiles of mRNA and protein expression level of m6A methylation regulators in BCa
After exploring the mutation of m6A regulators, we investigated the mRNA expression levels of m6A 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 Fig. 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 m6A 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 (Fig. 2C, D). Not only did the m6A 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 (Additional file 1: Fig. S3A). Immunohistochemistry staining images of m6A regulator proteins were retrieved from the Human Protein Atlas, revealing cellular sublocalization and intensity (Fig. 2E and Additional file 1: Fig. S3B). The above results revealed that cross-talk among m6A regulators might construct important modification patterns.
Identification of m6A regulators in two subgroups using consensus clustering
The empirical cumulative distribution function was plotted to analyze the optimal k value at which the cluster model achieved maximum stability (Fig. 3A and Additional file 1: Fig. S4A). The results showed that k = 2 gained the most powerful clustering efficacy, and the samples were divided into two subclusters using unsupervised clustering (Fig. 3B and Additional file 1: Fig. S4B-F). PCA analysis was used to judge the classification, and the two clusters could gather together (Fig. 3C). Prognostic analysis for the two clusters did not show a statistically significant difference, but a trend in overall survival (OS) (Fig. 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 (Fig. 3E). The m6A expression profiles showed that all m6A regulators were upregulated in cluster 2, except for IGF2BP1 (Additional file 1: Fig. S4G). GSVA enrichment analysis was performed to explore the biological behaviors of the two clusters. As shown in Fig. 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.
Characteristics of risk score pattern based on m6A regulators
To further explore the prognostic value of m6A regulators in BCa, we first performed a univariable Cox regression analysis on the expression levels of m6A regulators. The results showed that high expression of IGF2BP3 (hazard ratio [HR] = 1.2, 95% confidence interval [CI] = 1.1–1.3) and LRPPRC (HR = 1.0, 95% CI = 1.0–1.1) had worse survival outcomes in patients with BCa, whereas YTHDC1 (HR = 0.9, 95% CI = 0.9–1.0) and WTAP (HR = 1.0, 95% CI = 0.9–1.0) were regarded as protective markers for BCa (Fig. 4A). However, high expression of ZC3H13 seemed to have worse survival outcomes in patients with BCa regardless of P value.
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 m6A 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 (Fig. 4B, Additional file 1: Fig. S5A, Table S3). The final LASSO model with the best optimal lambda included six m6A 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 (Fig. 4C). As shown in Fig. 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 risk score pattern with AUC = 0.64, indicating that the risk score could predict the OS of patients with BCa (Fig. 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 (Fig. 4G and Additional file 1: Table S4). Moreover, the expression of risk genes (IGF2BP3, LRPPRC, and FTO) was higher in high-rsik patients, while YTHDC1, YTHDF2, and WTAP tended to be expressed in the low-risk group. Univariable and multivariable Cox analyses were performed to determine the independent prognostic value of the risk score pattern. Patient age, tumor T status, N status, and risk score were independent prognostic predictors in patients with BCa (Fig. 4H).
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 (Additional file 1: Fig. S5B, 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 (Additional file 1: Fig. S5D, E). Then, we investigated the distribution differences of somatic mutation between low and high risk score by using “maftools” R package. As shown in Additional file 1: Fig. S5F–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 m6A 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 (Additional file 1: Fig. S6A, B). Next, we investigated the correlation between the immune landscape and the risk score pattern. As shown in Fig. 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 that in the low-risk group (Fig. 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 (Fig. 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 m6A regulators in tumorigenesis and progression may be associated with tumor immunity.
Next, we explored the immune characteristics of independent m6A regulators in the risk group using the TIMER database to investigate the correlation between m6A regulator expression and immune cells, including B cells, CD8 + T cells, CD4 + T cells, macrophages, neutrophils, dendritic cells, and tumor purity. As shown in Fig. 5G and Additional file 1: Fig. S6C, 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 m6A regulators. As shown in Fig. 5H and Additional file 1: Figure S6D, IGF2BP3 amplification was associated with dendritic cells.YTHDC1 deletion was related to B cells, while amplification was 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 m6A regulators in the risk score model and several immune checkpoints (Fig. 5I). The results indicate that m6A 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 m6A 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 multivariable Cox regression model to predict OS in patients with BCa (Fig. 6A). The AUCs of the nomogram for predicting the 3- and 5-year OS were 0.69 and 0.70, respectively (Fig. 6B, C). The c-index of the nomograms for OS in the training set was 0.68. As shown in Fig. 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 were 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.75 and 0.72, respectively (Additional file 1: Fig. S7A, 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 (Additional file 1: Fig. S7C, 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 Fig. 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 (Fig. 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 (Fig. 7B, C).
The Oncomine database was used to determine the expression level of IGF2BP3. And as shown in Fig. 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 (Fig. 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 (Fig. 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 (Fig. 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 (Fig. 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 (Fig. 8G–H). The correlation between IGF2BP3 and PD-L1 expression was analyzed using BCa specimens. As shown in Fig. 8I, the 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.