The mRNA expression profile data and DNA mutation data (VarScan2) of breast cancer samples were obtained from The Cancer Genome Atlas (TCGA) database; 1,068 samples with complete prognostic information were included. A total of 1,094 breast cancer samples with complete prognostic information were obtained in the METABRIC cohort from cBioportal (www.cbioportal.org). Transcripts per million data were used for subsequent analysis. The data for a cohort of patients with metastatic urothelial cancer treated with anti-PD-L1 therapy were obtained from the R software package IMvigor210CoreBiologies (IMvigor210, version 1.0.0) (Mariathasan et al., 2018).

The gene set file “c2.cp.kegg.v6.2” was downloaded from the Molecular Signatures database (MSigDB; https://www.gsea-msigdb.org/gsea/index.jsp). The CIBERSORT (Newman et al., 2015) algorithm was used to evaluate the infiltration levels of immune cells in the sample, and 23 immune cell signatures were also used to evaluate the infiltration state of the TME (Charoentong et al., 2017; Jia et al., 2018). ESTIMATE was used to evaluate the immune and stromal scores of each sample (Yoshihara et al., 2013). Single-sample gene set enrichment analysis (ssGSEA) was used to calculate the enrichment scores of the samples using the GSVA package (Hänzelmann et al., 2013).

We obtained a total of 24 m6A regulators from the study by Zhang et al. (2020). A total of 1,997 immune genes were collected from 184 TME-related signatures (Wang et al., 2021). We identified m6A-related immune genes based on Pearson’s correlation analysis between the expression levels of m6A regulators and immune genes in breast cancer. |Pearson’s correlation coefficient| ≥ 0.5 and p < 0.001 were set as cutoff values. A Venn diagram (http://bioinformatics.psb.ugent.be/webtools/Venn/) was used to identify the m6A-related immune genes associated with prognosis.

Univariate Cox regression analysis was used to screen the m6A-related immune genes (p < 0.05 was set as the cutoff value). Iterative LASSO (least absolute shrinkage and selection operator) was used to screen the m6A-related immune genes for subsequent analysis (Wang et al., 2021). The number of iterations was 500, and genes with a frequency greater than 50 were consensus genes for the iteration LASSO. The order of frequency represents the degree of influence of these features, and these features were then incorporated into the multivariate Cox regression model; the inclusion was stopped when the area under the curve (AUC) value of the receiver operating characteristic curve (ROC) reached its peak. The m6A-related immune genes obtained were used to construct the m6A-IS. The m6A-IS was expressed as follows: m6A-IS = (coefficient mRNA₁ × mRNA₁ expression) + (coefficient mRNA₂ ×expression of mRNA₂) + …+ (coefficient mRNA _n × expression mRNA _n ). Visualization of the prognostic value of m6A-IS was obtained using a nomogram.

We arranged the m6A-IS from low to high, starting from one low-expressing patient by setting the loop, calculating the corresponding individual p-value and hazard ratio (HR), and saving the calculation results. The minimum p-value obtained was used to determine the grouping information, and the samples were then divided into two groups. The Kaplan–Meier method estimated the overall survival (OS) curve, and the difference between survival distributions was evaluated using the two-sided log-rank test implemented in the R package survival. The Kaplan–Meier survival curve constructed via the R package survmin.

The limma package in R was used to determine differentially expressed genes in the breast cancer cohorts. The ggplot2 package and ComplexHeatmap package were used to draw heatmaps and other maps. The R package forest plot was used to plot the forest plots. Wilcoxon’s rank-sum test was used to analyze the differences between the two groups. The Kruskal–Wallis test was used to compare differences between three or more groups. The maftools package was used to map the gene mutations. The AUC was quantified using the pROC R package. All statistical p-values were two-sided, and statistical significance was set at p < 0.05. All analyses were performed using the R software (version 4.0.2).

In the aggregate, 24 m6A regulators were identified in our study. In order to enrich the dysregulated genes and their correlations, the RNA transcriptomic datasets containing the next-generation sequencing (RNA sequencing, RNA-seq) data of 1,109 breast cancer tissues and 113 non-tumor tissues from TCGA project (TCGA-BRCA) was downloaded. Then, we explored the differential expressions in breast cancer tissues and normal mammary tissues in TCGA. Compared with normal tissues, there were 17 m6A regulators that were significantly differentially expressed in breast cancer tissues (Figure 2A). A correlation analysis was performed for the 24 m6A regulators (Figure 2B). The results showed that the expressions of the 24 m6A regulators were significantly correlated. Subsequently, we summarized the frequency of somatic mutations in the 24 m6A regulators in breast cancer. Of 985 samples, only 90 had mutations in the m6A regulators; indeed, there was a very low frequency of mutations (9.14%). Interestingly, we observed few mutations in breast cancer populations (Figure 2C). By performing univariate Cox regression analysis, we identified the affected prognosis of the 24 m6A-associated genes in patients with breast cancer (Figure 2D). Among the 24 m6A regulators, HNRNPC and RBM15B were identified as significant protective factors for survival (p < 0 0.05), and YTHDF3 was a risk factor for survival (p < 0 0.05).

Related research works have shown that the expressions of m6A regulators are associated with the heterogeneity of the TME (Han et al., 2019; Zhang et al., 2020). Therefore, we further analyzed the relationship between the m6A regulators and cell infiltration in the TME. We analyzed the association between 22 types of immune cells (CIBERSORT algorithm) and m6A regulators (Figure 2E). The heatmap showed that most of the significant enrichments were found in immune cells. The infiltration levels of Tregs, activated natural killer (NK) cells, and memory B cells were significantly related to the expressions of most m6A regulators. The infiltration levels of naive B cells, M1 macrophages, and memory CD4 T cells were significantly positively correlated with most m6A regulators. The results showed that the expressions of m6A regulators were significantly correlated with the levels of immune cell infiltration (p < 0.05); indeed, they play a non-negligible role in the regulation of the TME in breast cancer.

Correlation analysis was used for preliminary screening to further analyze which TME-related immune genes were related to m6A-regulated expression. Through correlation analysis, a total of 534 TME-related genes were found to be significantly positively or negatively correlated with m6A regulators (|Pearson’s correlation coefficient| ≥ 0.5 and p < 0.001) (Figure 3A). Univariate Cox regression analysis further revealed that 71 TME-related immune genes were significantly associated with breast cancer prognosis (p < 0.05; Figure 3A). Finally, 71 TME-related immune genes were included in the iterative LASSO algorithm for analysis. The results showed that, when these 28 TME-related immune genes were included (Supplementary Table S1), the model had the highest accuracy in predicting the prognosis (Figure 3B).

Multivariate Cox regression analysis was performed to construct the m6A-IS. The m6A-IS was expressed as follows: m6A-IS = (coefficient mRNA₁ × mRNA₁ expression) + (coefficient mRNA₂ × expression of mRNA₂) + …+ (coefficient mRNA _n × expression mRNA _n ). The results of the multivariate Cox regression analysis and m6A-IS are shown in Supplementary Tables S2 and S3. The results of the survival analysis showed that m6A-IS = 1.38 was the best cutoff value, and the low-m6A-IS group had a longer survival time in breast cancer than the high-m6A-IS group (Figure 3C). We also observed the same result from the METABRIC cohort (Supplementary Figure S1A).

To identify the signatures of the low- and high-m6A-IS subgroups, we analyzed the differences in the expressions of m6A modulators in the two subgroups and discovered that a large proportion of m6A regulators were differentially expressed in the high- and low-m6A-IS groups dramatically (Figure 3D). The clinicopathological information of the breast cancer patients has been shown in Table 1. The results of the Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis showed that the low-m6A-IS group was mainly enriched in the activation of some immune pathways, indicating that patients in this group may be with immune subtypes (Figure 3E). The high-m6A-IS group was mainly enriched in glucose metabolism and lipid metabolism pathways, indicating that patients in this group may be with metabolic subtypes (Figure 3E).

We investigated the expression levels of m6A-IS from the five breast cancer subtypes previously reported to further analyze its characteristics (Parker et al., 2009). We found that the Her-2 and luminal B subtypes had the highest m6A-IS, followed by the basal and normal subtypes. The luminal A subtype had the lowest m6A-IS level (Figure 4A and Supplementary Figure S1B). This is consistent with the results of previous studies: luminal A subtype has the best prognosis, Her-2 and luminal B subtypes have the worst prognosis, and the prognosis of the basal and normal subtypes can be characterized at a level between the best and worst. We classified the breast cancer samples based on the six immune subtypes previously identified using immunogenomic features (Thorsson et al., 2018). The results showed that m6A-IS expression was the lowest in the C3 subtype; meanwhile, it was the highest in the C4 subtype. This result suggests that the C3 subtype had the best prognosis, while the C4 subtype had the worst prognosis (Figure 4B). Thorsson et al. discovered that the C4 subtype displayed high M2 macrophage domination and low lymphocytic infiltrate, which induced poor outcomes. In contrast, the C3 subtype, which showed a type I immune response and remarkable Th17 signature, had favorable prognosis. Our results are in line with the conclusions of a previous report (Thorsson et al., 2018), indicating that the m6A-IS has good robustness in different breast cancer molecular subtypes. A Sankey diagram was used to visualize the relationship between m6A-IS and the two breast cancer subtypes (Figure 4C). In addition, between-group comparisons of age were performed using the Wilcoxon test, and there was a significant difference in the distributions of age. Patients older than 65 years were more likely to have higher scores than patients aged 65 years or younger (Figure 4D and Supplementary Figure S1C). Interestingly, we found that, with increased malignancy, the m6A-IS showed a significant stepwise increase in non-metastatic breast cancer; meanwhile, no significant increase was observed in patients with metastatic breast cancer (stage IV) (Figure 4D).

These results indicate that m6A-IS can be used to characterize the existing subtypes and clinical features of patients with breast cancer. However, the relationship between m6A-IS and TME remains unknown. For this reason, we found through ssGSEA that the low-m6A-IS group had a higher expression of stromal and immune scores than the high-m6A-IS group (Figure 4E), which is consistent with previous results (Figure 3E). In addition, the CIBERSORT algorithm (Figure 4F) and the enrichment analysis of 23 immune cells (Figure 4G and Supplementary Figure S1D) showed that the infiltration levels of B cells, CD8 T cells, and M1 tumor-associated macrophages in the low-m6A-IS group were significantly higher than those in the high-m6A-IS group. In contrast, the infiltration level of M2 tumor-associated macrophages increased significantly in the high-m6A-IS group. Finally, we found that the expression levels of the immune checkpoint molecules (PD1 and PD-L1) in the low-m6A-IS group were significantly higher than those in the high-m6A-IS group (Figure 4H). These results suggest that the low-m6A-IS group had increased immune cell infiltration and, simultaneously, had increased immunosuppression.

We further compared the two groups and found that the low-m6A-IS group had a higher expression of major histocompatibility complex (MHC)-related antigen-presenting molecules than the high-m6A-IS group (all p < 0.05; Figure 5A). We analyzed some co-inhibitory molecules to compare their expressions in the high and low-m6A-IS groups, and the results showed that the expression levels of these molecules were significantly correlated with m6A-IS. The results showed that the low-m6A-IS group had higher expressions of immunosuppressive molecules than the high-m6A-IS group (all p < 0.05; Figure 5A). Further analysis revealed that the low-m6A-IS group showed higher chemokine expressions (Figure 5B). These results indicate that the low-m6A-IS group had higher immunogenicity than the high-m6A-IS group. However, at the same time, there was a significant immunosuppressive state in the low-m6A-IS group. This suggests that there is a potential immune escape mechanism in low-m6A-IS.

Finally, the top 20 genes were further analyzed with the highest mutation frequency in the high and low-m6A-IS groups (Figure 5C). The results showed that PIK3CA and CDH1 had increased mutation frequencies in the low-m6A-IS group, while TP53 had an increased mutation frequency in the high-m6A-IS group. These results suggest that PIK3CA and CDH1 mutations may be related to a high immune infiltration, while TP53 mutations may be related to immunosuppression, which may require further analysis.

Univariate and multivariate Cox analyses represented that age, tumor stage, and the m6A-IS were independent prognostic predictors for patients with breast cancer (Figures 6A, B). Subsequently, we conducted a more detailed stratification of breast cancer patients based on clinical characteristics. Univariate Cox regression analysis revealed that m6A-IS has a good prognostic value in each subgroup (Figure 6C). Furthermore, the triple-negative (TNBC) subtype has the strongest tumor immunogenicity of all BC subtypes (Liu et al., 2018); meanwhile, m6A-IS can also be used to assess prognosis in this subtype (Figure 6C). These results indicate that m6A-IS, as an independent prognostic indicator, may be useful for clinical prognosis evaluation.

An ROC curve was constructed based on the m6A-IS. We found that the AUCs of the ROC curve at 1, 3, and 5 years were 0.816, 0.825, and 0.802, respectively (Figure 6D). Next, we verified whether the m6A-IS can predict the prognosis of breast cancer patients in the METABRIC cohort (Supplementary Figure S1E). It was shown that m6A-IS has good ability to predict prognosis. The previous results found that age and tumor stage were also independent factors; therefore, we combined them with the m6A-IS to build the final predictive prognosis model (Figure 6E). The results showed that the multi-predictor ROC had excellent accuracy for 1-, 3-, and 5-year OS (AUC = 0.903, 0.836, and 0.806, respectively). Based on the results of the logistic regression analysis and the ROC curves, a nomogram was graphically depicted (Figure 6F). By calculating the total scores of each selected variable, the survival of individual breast cancer patients at 1, 3, and 5 years can be easily estimated by plotting a vertical line between the total points and each prognosis axis of the nomogram.

Recently, the connection of m6A modification with PD-L1 expression has been reported. In an anti-PD-L1 cohort (IMvigor210 cohort), the low-m6A-IS group exhibited significant clinical benefits and a markedly longer survival time (Figure 7A). Compared with those with high m6A-IS, patients with low m6A-IS showed significant survival advantages and clinical response to anti-PD-L1 therapy (Figures 7B–D). Next, we examined the differences in m6A-IS among the different immune phenotypes in the IMvigor210 cohort. Interestingly, patients in the higher m6A-IS group were remarkably relevant to the exclusion and desert immune phenotypes, and it was difficult to achieve antitumor effects of the checkpoint inhibitors in these phenotypes (Figure 7E). Based on accumulating evidence, high levels of tumor mutation burden (TMB) and PD-L1 expression have long-lasting clinical responses to immunotherapy. A high TMB is generally considered to be a preexisting adaptive immune response to the tumor, and patients with a high TMB who received PD-1 blocking immunotherapy showed an improved response and enhanced clinical efficacy compared to patients with a moderate or low TMB (Hellmann et al., 2018). Further analysis revealed that m6A-IS was significantly negatively correlated with TMB (Figure 7F). These results suggest that the quantification of the modification patterns of m6A is a potential and robust biomarker for the assessment of prognosis and clinical response to immunotherapy (Figure 7G). In conclusion, this study emphasizes that m6A methylation modification is significantly correlated with tumor immunophenotype and anti-PD-L1 immunotherapy response, and established m6A modification characteristics may help predict the response to anti-PD-L1 immunotherapy. An illustration of this study is shown in Figure 8.