m6A-Related lncRNA to Develop Prognostic Signature and Predict the Immune Landscape in Bladder Cancer

Abnormal m6A methylation plays a significant role in cancer progression. Increasingly, researchers have focused on developing lncRNA signatures to evaluate the prognosis of cancer patients. The specific function of m6A-related lncRNAs in the prognosis of bladder cancer patients and the immune microenvironment of bladder cancer remains elusive. Herein, we performed a comprehensive analysis of m6A-related lncRNA prognostic values and their association with the immune microenvironment in bladder cancer using the TCGA dataset. A total of 9 m6A-related lncRNAs were dramatically correlated with overall survival outcomes in bladder cancer. Two molecular subtypes (cluster 1 and cluster 2) were identified by consensus clustering for 9 m6A-related prognostic lncRNAs. Cluster 1 was significantly correlated with poor prognosis, advanced clinical stage, higher PD-L1 expression, a higher ESTIMATEScore and immuneScore, and distinct immune cell infiltration. GSEA revealed the enrichment of apoptosis and the JAK-STAT signaling pathway in cluster 2. A prognostic risk score was constructed using 9 m6A-related prognostic lncRNAs, which functioned as an independent prognostic factor for bladder cancer. Moreover, bladder cancer patients in the low-risk score group had a higher pN stage, pT stage, and clinical stage and a lower tumor grade and immuneScore. The risk score was correlated with the infiltration levels of certain immune cells, including B cells, plasma cells, follicular helper T cells, regulatory T cells, resting NK cells, neutrophils, M0 macrophages, M1 macrophages, and M2 macrophages. Collectively, our study elucidated the important role of m6A-related lncRNAs in the prognosis of bladder cancer patients and in the bladder cancer immune microenvironment. The results suggest that the components of the m6A-related prognostic lncRNA signature might serve as a crucial mediator of the immune microenvironment in bladder cancer, representing promising therapeutic targets for improving immunotherapeutic efficacy.


Background
Bladder cancer (BC) is one of the most frequent urinary malignancies in China and has a high incidence rate [1]. Transitional cell carcinoma ranks approximately 90% of bladder cancer. Most of the BCs are nonmuscle invasive (MI) BCs, which readily relapse and develop into muscle invasive bladder cancers [2]. Approximately 30% of bladder cancers are MIBCs, whose gold standard for their treatment is radical cystectomy and pelvic lymph node dissection [3]. e prognosis of MIBC patients is poor [4]. Although several clinical features and molecular biomarkers have been applied for the prognosis of bladder cancer patients, these approaches are all limited to some extent. us, it is necessary to construct a new predictive model and identify new prognostic markers for bladder cancer.
N6-methylandenosine (m6A) modification, one of the most common epigenetic methylation modifications, plays an important part in many biological processes, including RNA splicing, export, stability, and translation [5]. m6A methylation is significantly related to the levels of intracellular methyltransferases ("writers") and demethylases ("erasers"), while binding proteins ("readers") interact with m6A methylation sites to perform several biological functions [6]. Abnormal m6A methylation is involved in the progression of cancer by regulating various biological processes, including cell differentiation, immunoreaction, and miRNA editing [7]. Moreover, certain regulators of m6A methylation serve as prognosis biomarkers for certain cancers. For example, YTHDF1 and HNRNPC were suggested to be prognostic biomarkers of colon cancer [8]. Another study revealed that HNRNPC is a marker for prognosis in glioblastoma multiforme and contributes to carcinogenesis [9]. Another bioinformatics study revealed that the regulators of m6A methylation could promote tumor progression and affect the overall survival (OS) of BC patients [10].
Long noncoding RNAs (lncRNAs), a category of RNAs with transcript lengths over 200 nucleotides, regulate 70% of gene expression in mammals by interacting with DNA, RNA, and proteins [11]. Aberrant lncRNA expression is associated with the regulation of the proliferation, invasion, and apoptosis of tumor cell, thus affecting the pathogenicity of bladder cancer and patient prognosis. e lncRNA LINC00641 has been shown to act as a prognostic biomarker and to inhibit bladder cancer progression [12]. Another study revealed that the lncRNA CCAT1 could promote tumor cell biological progression in bladder cancer [13]. However, the overall functions of m6A methylation-related lncRNAs in cancers remain a mystery. erefore, exploring m6A methylation-related lncRNAs and identifying prognostic biomarkers among these lncRNAs are of significance.
Based on e Cancer Genome Atlas (TCGA) database, the identification of prognostic gene signatures has become possible [14]. Using the lncRNA expression profiles of TCGA bladder cancer cohort, we identified prognostic biomarkers based on m6A methylation-related lncRNAs and constructed a prognostic model for bladder cancer by using bioinformatics methods. We also evaluated the correlation between the signature lncRNAs and immune infiltration in bladder cancer.

Datasets.
e mRNA expression profiles of BC patients were isolated from the TCGA database, and we obtained data on 412 BC tissues and 19 normal bladder tissues. Moreover, the clinical characteristics of the BC patients, including age, grade, and survival status, were obtained from the TCGA.

Identification of Prognosis-Related lncRNA.
We then explored the prognostic value of m6A-related lncRNAs using the "survival" package with a p value of 0.0001. We then generated a forest plot according to the data of the univariate Cox analysis. Moreover, we also constructed a heatmap using the "pheatmap" package to show the expression of those lncRNAs with significant prognosis value in BC tissues and normal tissues.

Bioinformatics
Analysis. m6A methylation-related lncRNAs with significant prognostic value were selected for further analysis. We then classified bladder cancer into different subtypes using the "ConsensusClusterPlus" package (1,000 iterations and resampling rate of 80%). e gene expression patterns between each subtype were evaluated and visualized with the "pheatmap" package. Gene set enrichment analysis (GSEA) was conducted to detect the function of each bladder cancer subtype with a simulation of 500 and FDR of 0.05. e ESTIMATE algorithm was utilized to calculate the immuneScore, stromalScore, and ESTIMATEScore for each bladder cancer patient using the "estimate" package. e infiltration abundance of 22 immune cell types in each BC subtype was visualized with the "vioplot" package. e prognostic signature of m6A methylation-related lncRNAs was developed by LASSO regression analysis. By LASSO regression analysis, the coefficients of each bladder cancer case were calculated with the following computational equation: risk score � sum of coefficients × the lncRNA expression. e risk score of all the BC patients was computed in the training and test cohorts. Subsequently, the patients were separated into high-and low-risk groups with the cutoff point set as the median value of the risk score. Moreover, the correlation between the risk score and the abundance of immune cells was also calculated and "ggplot2" was used to visualize the result.

Statistical Analysis.
All statistical tests were conducted with R version 4.0.1. e Mann-Whitney U test was applied to explore the mRNA levels of m6A methylation-related lncRNAs.
e differences between two subgroups were evaluated with Student's t-test.
e chi-square test was performed to compare categorical variables in the training and test cohorts. Survival curves were drawn using the Kaplan-Meier method. Pearson correlation tests were performed to explore the correlation among subtypes, clinicopathological features, risk scores, immune checkpoints expression, and immune infiltration levels. Univariate and multivariate analyses were performed with Cox regression models to explore the independent prognostic value of the risk scores integrated with other clinical features. p < 0.05 indicates statistical significance. Figure 1 displays the workflow of the current study. e clinical characteristics of bladder cancer patients were downloaded from the TCGA database and rearranged. A total of 412 cases of TCGA bladder cancer were obtained, and the clinical characteristics are shown in Table 1. Based on the lncRNA annotation file in GENCODE, we identified 14087 lncRNAs in the TCGA bladder cancer dataset. After extracting the expression profiles of 23 m6A methylation regulators in TCGA bladder cancer dataset, we evaluated the association between 23 m6A methylation regulators and 14087 lncRNAs. lncRNAs that were associated with one or more of the 23 m6A methylation regulators (|Pearson R| > 0.3 and p < 0.01) were defined as m6A-related lncRNAs. As a result, we obtained 762 m6A-related lncRNAs.

Identification of m6A-Related lncRNAs in Bladder Cancer.
e coexpression network between 23 m6A methylation regulators and 762 m6A-related lncRNAs is presented in Figure 2(a). Based on the prognostic value of these m6A-related lncRNAs, a univariate Cox regression analysis was performed to identify m6A-related prognostic lncRNAs with a p value of 0.0001. e data revealed that 9 m6A-related lncRNAs were markedly associated with OS in bladder cancer patients (Figure 2(b) and Table 2). We then analyzed the mRNA level of the 9 m6A-related prognostic lncRNAs in bladder cancer, which indicated the mRNA level of 8 m6A-related prognostic lncRNAs (PTOV1-AS2, AC116914.2, EHMT2-AS1, AC004148.1, AL136295.2, KCNQ1OT1, AC104564.3, and AC073534.2) was upregulated and the expression of 1 m6A-related prognostic lncRNA (ATP1B3-AS1) was decreased in BC tissues compared with their expression levels in bladder tissues (Figure 2(c)).

Consensus Clustering Categorized
Patients according to m6A-Related Prognostic lncRNAs. Consensus clustering was utilized to separate bladder cancer patients into subgroups according to the expression of m6A-related prognostic lncRNAs. k � 2 was found to be optimal clustering stability from k � 2 to 9 based on the similarity displayed by the expression levels of m6A-related prognostic lncRNAs (Figure 3(a)). e cumulative distribution function, increment in the AUC, and tracking plot of subgroups for k � 2-9 are presented in Supplementary Figures 1B and 1C, respectively. A total of 406 bladder cancer patients were separated into cluster 1 and cluster 2 (Figure 3(a)). e OS rate of bladder cancer patients in cluster 1 was worse than that of those in cluster 2 ( Figure 3(b), p � 0.022). Moreover, we found that cluster 1 was markedly related to an advanced clinical stage (Figure 3(c), p < 0.05).

Consensus Clustering Correlated with Immune
Infiltration. To explore the role of m6A-related prognostic lncRNAs in the bladder cancer immune microenvironment, we then analyzed the difference in the immuneScore and immune cell infiltration level between cluster 1 and cluster 2. e average immuneScore (Figure 4(a), p � 8.5e − 12 ), stromalScore ( Figure 4(b), p � 3.5e − 13 ), and ESTIMATE-Score ( Figure 4(c), p � 1.3e − 13 ) were higher in cluster 1 than in cluster 2. e infiltration abundance of 22 types of immune cells in each cluster is shown in Figure 4(d). As shown in Figures 4(e)-4(g), cluster 1 had a higher abundance of CD4 memory-activated T cells (p � 0.019), a lower abundance of regulatory T cells (p � 0.0012), and NK resting cells (p � 0.038) compared to cluster 2. Moreover, we also found a positive correlation between macrophage M1 and macrophage M2 (Supplementary Figure 2). We then detected the mRNA level of immune checkpoints in each subtype and their correlation with m6A-related prognostic lncRNAs. e expression of PD-L1, CTLA4, HAVCR2, LAG3, PDCD1, PDCD1LG2, and TIGIT was higher while SIGLEC15 expression was significantly lower in cluster 1 than in cluster 2 ( Figure 5(a), p < 0.001). We also found that the expression of CTLA4, TIGIT, and SIGLEC15 in BC tissues was significantly elevated ( Figure 5

Enrichment Analysis of Each Bladder Cancer Subtype.
GSEA was conducted to clarify the potential regulatory mechanisms leading to the differences between the two clusters of bladder cancer patients. Some cancer-related hallmarks, including the spliceosome and the mTOR and Notch signaling pathways, were significantly associated with cluster 1 (Supplementary Figure 5). Apoptosis and the chemokine, Toll-like receptor, and JAK-STAT signaling pathways were associated with cluster 2 (Supplementary Figure 6).

Construction of Prognostic
Signature. LASSO Cox analysis was applied to construct the prognostic signature using the 9 identified m6A-related prognostic lncRNAs. e coefficient and partial likelihood deviance of prognostic signature are shown in Figures 6(a) and 6(b). e 406 bladder cancer patients were randomly separated to training cohort and test cohort. e risk score of each patients was calculated using the following equation: . e TCGA bladder cancer patients were separated into high-risk and low-risk groups. And OS curve suggested a poor prognosis in patients in the high-risk group compared with those in the low-risk group in both the training cohort ( Figure 6(c)) and test cohort ( Figure 6(e)), with an AUC of 0.651 in the training cohort ( Figure 6(d)) and 0.737 in the test cohort ( Figure 6(f )), indicating that the signatures comprising 9 m6A-related prognostic lncRNAs had a favorable discrimination performance in predicting bladder cancer patient prognosis.
To further explore the factors affecting the prognosis of bladder cancer, we then performed univariate and multivariate Cox regression analyses. e univariate Cox regression analysis showed that age, sex, clinical stage, pT stage, pN stage, and risk score were related to the prognosis of bladder cancer patients in the training cohort ( Figure 6(g), all p < 0.05). Moreover, the multivariate Cox regression analysis revealed that the risk score (p < 0.001) was still significantly associated with the prognosis of bladder cancer patients ( Figure 6(h)). In the test cohort, the univariate and multivariate Cox regression analyses were also performed, revealing the risk score as the factor affecting the prognosis of bladder cancer patients (Figures 6(i) and 6(j)). Figure 7 shows the risk score distribution (Figures 7(a) and 7(b)) and survival status (Figures 7(c) and 7(d)) of each BC patient in the training and test cohorts. e mRNA levels of lncRNA, including KCNQ1OT1 and ATP1B3-AS1, were downregulated in the high-risk group, whereas protective m6Arelated lncRNAs, including PTOV1-AS2, AC116914.2, EHMT2-AS1, AL136295.2, AC104564.3, and AC073534.2, were expressed at low levels in the high-risk group in both the bladder cancer training and test cohorts (Figures 7(e) and 7(f)).
We then verified the prognosis of the risk score in different groups of BC patients, which found that the highrisk group had a poor prognosis compared with the low-risk group in bladder cancer patients aged >50 years (Supplementary Figure 7A

Risk Score Associated with Clinical Characteristics.
e heatmap in Figure 8(a) revealed the mRNA level of lncRNAs in the high-and low-risk groups in TCGA BC cohort. e results suggested that the expression levels of KCNQ1OT1 and ATP1B3-AS1 were higher in the high-risk group, whereas the expression levels of PTOV1-AS2, AC116914.2, EHMT2-AS1, AL136295.2, AC104564.3, and AC073534.2 were lower in the high-risk group than in the low-risk group in the bladder cancer cohort (Figure 8(a)). e heatmap also demonstrated the differences in terms of pN stage, pT stage, clinical stage, grade, immuneScore, and cluster subtype between the high-and low-risk groups (Figure 8(a), all p < 0.01). More specifically, bladder cancer in stage N1-N3 was related to a higher risk score than that in stage N0 (Figure 8  immuneScore group was significantly higher than that of the low immuneScore group (Figure 8 . e results also revealed a different risk score between the two clusters ( Figure 8(g), p � 3.9e − 16 ). ese results demonstrated that the risk score was linked to the clinical characteristics in BC.

Discussion
m6A methylation is the most common form of mRNA modification and plays a vital role in regulating gene expression at the posttranscriptional level [17]. Abnormal m6A methylation plays a vital role in the progression of cancer by regulating many biological processes, including cell differentiation, immunoreaction, and miRNA editing [7]. Increasingly, researchers have focused on developing lncRNA signatures to evaluate the prognosis of cancer patients [18]. However, limited studies have been performed to study the role of m6A-related lncRNAs in patient prognosis and the immune microenvironment of malignancies, including bladder cancer. In our study, we performed a comprehensive analysis of the expression, prognostic value, and effects on the immune microenvironment of m6A-related lncRNAs in bladder cancer.
We first identified m6A-related lncRNAs by constructing a coexpression network, and a total of 762 m6Arelated lncRNAs were obtained. is was followed by univariate Cox regression analysis for the identification of m6Arelated prognostic lncRNAs. As a result, a total of 9 m6Arelated lncRNAs were significantly related to overall survival outcomes in bladder cancer patients, and the expression of all of these lncRNAs was upregulated or downregulated in bladder cancer tissues compared with that in normal tissues. us, these 9 m6A-related prognostic lncRNAs were selected for further analysis. Based on consensus clustering for the 9 m6A-related prognostic lncRNAs, two subtypes (cluster 1 and cluster 2) of bladder cancer were identified. Interestingly, the cluster 1/2 subtype stratification showed a significant correlation with the prognosis and clinical stage of bladder cancer patients and PD-L1 expression. Moreover, we found higher immuneScores, stromalScores, and ESTI-MATEScores in cluster 1 compared with those in cluster 2. Interestingly, these results were consistent with the conclusion of the previous study, which revealed that bladder cancer patients with high immuneScore and stromalScore had a poor overall survival rate [19].
e GSEA results revealed that the spliceosome and mTOR signaling pathways were significantly associated with cluster 1. Apoptosis and the JAK-STAT signaling pathway were associated with cluster 2. e JAK-STAT signaling pathway plays a significant role in many biological processes, including cell division, apoptosis, and immune regulation [20,21]. e mTOR signaling pathway is one of the most investigated therapeutic targets in bladder cancer research [22]. Apoptosis inhibition is an important hallmark of bladder cancer. e above results revealed that cluster 2 bladder cancer patients had a better overall survival outcome than those in cluster 1.
us, we hypothesized that apoptosis and the JAK-STAT signaling pathway were more active in cluster 2 than in cluster 1.   LASSO Cox analysis was performed to construct the m6A-related lncRNA prognostic signature based on the 9 m6A-related prognostic lncRNAs. Based on the risk scores for the overall survival outcome of each patient calculated by the LASSO algorithm, bladder patients were separated into high-and low-risk groups. We found that bladder cancer patients in the high-risk group had a poor prognosis compared with those in the low-risk group.
is signature of 9 m6A-related lncRNAs had a favorable discrimination performance for predicting bladder cancer patient prognosis.
is result is consistent with the results of previous studies, which have revealed that   lncRNA-related signatures play a vital role in predicting the prognosis of bladder cancer patients. Wang et al. performed a bioinformatics analysis to identify a signature of seven immune-related lncRNAs, which could serve as prognostic biomarkers for bladder cancer [23]. Another prognostic signature based on immune-related lncRNAs can be used to predict the prognosis and immunotherapeutic response of bladder cancer patients [24]. Moreover, an eight-lncRNA signature was suggested as a candidate prognostic biomarker for bladder cancer [25]. In our study, we developed a prognostic signature based on 9 previously unstudied m6A-related prognostic lncRNAs for bladder cancer. Univariate and multivariate Cox regression analyses revealed that the risk score was an independent factor for predicting the prognosis of bladder cancer patients. We then analyzed the association between the prognostic risk scores and clinical characteristics of bladder cancer patients. We found that patients with bladder cancer in stage N1-N3 and with a low tumor grade had a higher risk score compared with those with N0 and high-grade tumors. e risk score increased as the pT stage and the clinical stage increased. Interestingly, the risk score of the high immu-neScore group was significantly higher than that of the low immuneScore group. is finding was consistent with the results of a previous study, which indicated that bladder cancer patients with a high immuneScore had a poor overall survival compared with those with a low immuneScore [26]. e tumor immune microenvironment exerts a vital function in tumorigenesis and cancer progression, and its heterogeneity can influence multiple factors, including patient prognosis and therapeutic response [27,28]. e results of a previous study indicated that immune cell infiltration can regulate tumor progression and metastasis, thus affecting patient prognosis [29,30]. Another important finding of our study is that the risk score was correlated with the infiltration levels of certain immune cells, including B cells, plasma cells, follicular helper T cells, regulatory T cells, resting NK cells, neutrophils, M0 macrophages, M1 macrophages, and M2 macrophages. As the risk score increased, the infiltration levels of neutrophils increased. is result was consistent with previous results, which suggested that tumor-infiltrating neutrophils were associated with a high risk of disease recurrence and poor overall survival outcomes [31]. A previous study revealed better overall survival rates in bladder cancer patients with high regulatory T-cell infiltration levels [32]. In our study, we found a negative correlation between regulatory T-cell infiltration levels and the risk score. Bladder cancer patients with a high risk score had poor overall survival.
is result was consistent with previous results.
ere are some limitations of the current study. First, our results were obtained by analyzing data from the TCGA, and it would be beneficial to verify our findings using the GEO database.
e regulatory mechanism of the m6A-related lncRNA prognostic signature warrants further investigation to determine methods to reshape the immune microenvironment and improve precision immunotherapy for bladder cancer.   , immuneScore (f ), and cluster (g) of bladder cancer patients. * p < 0.05, * * p < 0.01, and * * * p < 0.001.

Conclusions
In conclusion, this study systematically analyzed the expression, prognostic value, and effects on the immune microenvironment of m6A-related lncRNAs in bladder cancer. By consensus clustering for m6A-related prognostic lncRNAs and construction of a prognostic signature, our study elucidated the important role of m6Arelated lncRNAs in patient prognosis and the immune microenvironment in bladder cancer. e results suggest that the m6A-related prognostic lncRNA signature might serve as a crucial mediator of the immune microenvironment in bladder cancer, representing promising therapeutic targets for improving immunotherapeutic efficacy. e findings of our study provided potentially theoretical foundation for future animal and clinical studies about m6A-related lncRNAs as promising therapeutic targets for bladder cancer.

Data Availability
All data generated or analyzed during this study are included in this published article.

Conflicts of Interest
e authors declare that they have no conflicts of interest.

Authors' Contributions
Zuwei Li and Peiyuan Huang performed data analysis work and aided in writing the manuscript. Yuwu Li and Weizhang Zhong designed the study and assisted in writing the manuscript. All authors read and approved the final manuscript.

Supplementary Materials
Supplementary Figure