M6A expression in BC
In this study, three public cohort datasets and one clinical dataset were included, which were as follows: 1091 BC cases and 113 non-tumor normal cases from the TCGA date, 1985 cases from the METABRIC date, BC cases from the KM Plotter website, and 134 BC cases from one clinical dataset.
In this study, the following 28 m6A regulators were studied: 10 writers were METTL3, METTL14, METTL16, WTAP, ZC3H13, RBM15B, RBM15, CBLL1, KIAA1429, and NSUN2; 1 eraser: ALKBH5; and 17 readers: YTHDC1, YTHDC2, YTHDF1, YTHDF2, YTHDF3, FXR2, FXR1, EIF3A, EIF4G2, IGF2BP1, IGF2BP2, IGF2BP3, ELAVL1, G3BP1, HNRNPA2B1, LRPPRC, and ABCF1.
Genetic differences between BC and adjacent normal tissues
Based on the TCGA data cohort, the analysis results demonstrated that 21 identified genes were differentially expressed between BC tissues and adjacent normal tissues, including 14 up-regulated genes and 7 down-regulated genes (Fig. 2a). These up-regulated genes were IGF2BP1, IGF2BP3, RBM15, CBLL1, KIAA1429, FXR1, ELAVL1, NSUN2, ABCF1, LRPPRC, YTHDF2, YTHDF1, HNRNPA2B1, and EIF4G2. These down-regulated genes were IGF2BP2, METTL14, METTL16, ZC3H13, YTHDC1, WTAP, and EIF3A. In addition, no significant differences were observed for the other seven genes (Fig. 2a).
Genetic differences between different types of BC
Based on the TCGA cohort, the ANOVA showed that the expression of 28 m6A regulators in various BC types was significantly different (Fig. 2b). The results showed that the expression of m6A regulators was significantly correlated with different molecular typing of BC, except EIF3A and YTHDF2 (Fig. 2b, supplementary Figure 1). Specifically, IGF2BP1, IGF2BP2, and IGF2BP3 were highly expressed in TNBC, while they were expressed at low levels in HER2+, luminal A, and luminal B BC (Fig. 2b).
Interaction of m6A regulators
By analyzing the TCGA cohort, the results indicated that all enrolled 28 studied genes exhibited gene-gene interactions. The results also showed negative gene-gene expression interaction among IGF2BP1, IGF2BP2, and IGF2BP3, and strong positive gene-gene expression interactions among the remaining 25 genes (Fig. 2c).
Genetic mutation and copy number variations (CNVs) of the m6A regulators
According to the TCGA data cohort, the results showed that these enrolled 28 studied genes rarely had genetic mutations, and the overall mutation frequency was only 8.45% (71/840). The statistical analysis indicated that the mutation frequency of KIAA1429, which had the highest frequency, was only approximately 1%, and as for the other genes, the frequencies were below 1% (Fig. 3a). We also used COSMIC (the Catalogue Of Somatic Mutations In Cancer, https://cancer.sanger.ac.uk/cosmic/) to explore the variants of m6A regulators in BC. The results demonstrated that the mutation frequency of IGF2BP2 was 6.3%, IGF2BP3 was 6.0% and LRPPRC was 4.2%. The mutation frequencies of the remaining m6A regulators remained between 0.5% and 2.9%.
In this study, we also analyzed the following CNVs of m6a regulators: homozygous deletion (HOMD), heterozygous deletion (HETD), neural (NEUT), gain, and amplification (AMP). The top three genes with CNV gains were KIAA1429, YTHDF3, and YTHDF1; the top three genes with HOMD mutations were ZC3H13 and FXR2, and WTAP; the top three genes with HETD mutations were FXR2, METT16, and ALKBH5; and the top three genes with AMP mutations were KIAA1429, YTHDF3, and IGF2BP1 (Fig. 3b). Moreover, the results indicated a clear relationship between m6A regulators and CNVs. (Supplementary Figure 2) Supplementary Figure 2 also showes that the expression trend of some genes was from low to high. Furthermore, the expression of ABCF1, ALKBH5, ELAVL1 and FXR1 expression was gradually increased from deletion, loss, neural, gain, to amplification (Supplementary Figure 2a-d).
Prognostic significance of m6A regulators in BC
Prognostic significance from TCGA data and METABRIC data
In the TCGA cohort, 816 eligible samples were selected according to the following inclusion criteria: 1. overall survival (OS) > 6 months; and 2. enrolled sample has complete stage, molecular classification and relapse-free survival (RFS) data. The Cox univariate analysis showed that BC patients with CBLL1 high expression had a better RFS than those with low CBLL1 expression. (HR95%CI = 0.51(0.26–0.97)) (Fig. 4a).
In the METABRIC cohort, the inclusion criteria of the population were as follows: 1. OS > 12 months; 2. The sample has complete grade, TNM stage and molecular classification data. The Cox univariate analysis showed that the expression of the following 4 genes was related to better OS: CBLL1 (HR95%CI = 0.70 (0.54–0.91)), FXR2 (HR95%CI = 0.75 (0.58–0.98)), RBM15B (HR95%CI = 0.76 (0.61–0.94)), and KIAA1429 (HR95%CI = 0.76 (0.61–0.94)). In contrast, YTHDF1was related to poor OS. (HR95%CI = 1.38 (1.13–1.70)) (Fig. 4b).
In summary, higher CBLL1 expression in BC was associated with better RFS in the TCGA cohort and better OS in the METABRIC cohort. That is, higher CBLL1 expression is correlated with better prognosis.
Relationship between CBLL1 and clinicopathological characteristics of BC
Based on the TCGA cohort, the results showed that CBLL1 expression was not significantly different in different stages of BC (P = 0.16, Fig. 5a). Additionally, we analysed the relationship between CBLL1 and the PAM50 molecular subtypes of BC. The expression of CBLL1 in luminal BC was significantly higher than that in HER2 positive BC and TNBC (P < 0.05, Fig. 5b). However, CBLL1 expression in luminal A BC was not significantly different from that in luminal B and PAM50-normal BC. (P > 0.05, Fig. 5b).
Prognostic significance of CBLL1 in BC
From the TCGA cohort, a total of 816 samples were included to analyse the relationship between CBLL1 expression and the prognosis of BC patients. The results showed that patients with high CBLL1 expression had better RFS (P = 0.005, Fig. 6a). In the METABRIC cohort, the results of the analysis with the quartile method showed that patients with high CBLL1 expression had better OS with the quartile method involved (P = 0.02, Fig. 6b). Moreover, based on that cohort in the Kaplan-Meier Plotter website, we confirmed that high CBLL1(affymetrix ID: 227187_at) expression is correlated with better RFS (HR = 0.64 (0.55–0.75), log-rank p = 2.7e− 08) (See Fig. 6c). In addition, the analysis of 134 selected clinical samples showed that patients with high CBLL1 expression in the cytoplasm also had better OS (HR = 0.46 (95%CI 0.24–0.89), log-rank p = 0.046) (See Fig. 6d). In this study, CBLL1 −/+ was defined as low expression, and CBLL1++/+++ was defined as high expression (See Fig. 6e).
Functional analysis of CBLL1 in luminal BC
In a previous study, Hakai, as a coregulator of oestrogen receptor alpha, was found to play a negative role in the development and progression of BC cells [23]. Thus, we performed functional analysis of CBLL1 in ER- or PR-positive BC. Based on the TCGA cohort, a total of 597 ER- or PR-positive and HER2-negative BC patients patients were included according to the inclusion criteria. Altogether, 521 genes with differential gene expression (log FC > 0.8 or < − 0.3) were included in the analysis.
A total of 597 patients was were divided into the CBLL1-high (CBLL1-H) and CBLL1-low (CBLL1-L) groups according to the median value of CBLL1 expression. The heat map compared the pathological features and signalling pathway features of BC in the CBLL1 high-expression group with those in the CBLL1-low expression group. The following parameters were observed in these two groups: pTNM staging, PAM50, PIK3CA, GATA3, apoptosis-related pathways, ESR1-related pathways, and immune-related pathways (Fig. 7a). The expression of the following apoptosis regulators was significantly different in the CBLL1-H and CBLL1-L group, as follows: POMK, NRIP1, GTF2I, SEMA3C, VPS13C, RAB27B, PIK3C2A, ITPR2, HIPK2, and ADAM9. In addition, the expression of ESR1 regulators,including ZNF770, AFF3, PRLR, SLC7A2, and CLSTN2, was also significantly different in CBLL1-H and CBLL1-L group. Furthermore, immune regulators, including CALML5, ISG15, S100A8, LTB, CYBA, S100A9, TNFRSF4, SCT, and IFI27, were significantly different between the CBLL1-H and CBLL1-L groups (Fig. 7a).
Gene set enrichment analysis (GSEA) demonstrated that the CBLL1-related differentially expressed genes were significantly enriched in the ESR1-related signalling pathway, cell apoptosis related pathway and immune-related pathway. In particular, CBLL1 expression promoted the upregulation of genes related to ESR1 and cell apoptosis-related pathways (Fig. 7b). Moreover, single-sample gene set enrichment analysis (ssGSEA) demonstrated that low CBLL1 expression increased the occurrence of tamoxifen resistance and tumor-associated hypoxia. High CBLL1 expression induced apoptosis, mitotic cell death and the upregulation of ESR1 (Fig. 7c).