N6-Methyladenosine Regulators Are Involved in the Progression of and Have Clinical Impact on Breast Cancer

Background

Breast cancer is one of the most common and lethal malignancies in women worldwide. It is estimated that about 1 806 590 new cases and 606 520 deaths will occur in 2020 [1], and therefore immediate action is needed. Although enormous progress has been achieved in the past decades, with the development of surgical skills, diagnosis, and treatment modalities, many questions remain to be answered, such as specific risk factors of breast cancer, limits of immunotherapy, and identification of patients in need of treatment [2].

Epigenetics, mainly including histone modification, chromosome remodeling, acetylation, and methylation, could influence genes’ transcription and translation without changing sequences of DNA [3]. Epigenetics is also involved in the carcinogenesis, development, and prognosis of breast cancer [4]. As a notable methylation form, N6-methyladenosine (m6A) was first discovered and studied in the 1970s and since then has been viewed as the main methylation modification in eukaryotic cells [5]. The m6A regulators are a group of molecules that regulate m6A by adding, erasing, or recognizing the methyl group at the N6 position of adenosine. The m6A process depends on these regulators: methyltransferases (“writers”), demethylases (“erasers”), and binding proteins (“readers”). “Writers” install m6A units to RNA and are mainly composed of Wilms tumor 1- associated protein (WTAP) [6], methyltransferase-like 3 (METTL3) [7], and METTL14 [7]. “Erasers” can remove the m6A unit from molecules, and consist of α-ketoglutarate-dependent dioxygenase alkB homolog 5 (ALKBH5) [8] and fat-mass and obesity-associated protein (FTO) [9]. The modification of m6A is recognized by the “readers”, which contain heterogeneous nuclear ribonucleoproteins (HNRNPs), YT521-B homology (YTH) family members, and insulin-like growth factor 2 mRNA-binding proteins (IGF2BPs) [10]. In addition to the above, METTL16, RBM15, KIAA1429, RBM15B, EIF3, ZC3H13, and CBLL1 were included in our analysis. In all, 20 m6A regulators were selected in our study, based on previous research and reviews [11,12].

A growing number of studies are focusing on the role of m6A in breast cancer. Data have shown that m6A modification regulates breast cancer progression [13] and is correlated with prognosis [14]. The knockdown of METTL3 could decrease the m6A methylation level, and therefore inhibit breast cancer proliferation and tumor growth [13]. It has been demonstrated that upregulated YTHDF1, YTHDF3, and KIAA1429 can predict poor overall survival (OS). Overexpression of YTHDF3 is an independent factor for OS in patients with breast cancer [14]. To fully understand the role of m6A in breast cancer, we used data from The Cancer Genome Atlas (TCGA) to comprehensively analyze RNA expression patterns and their potential relevance to clinical characteristics in breast cancer. Our goal was to provide new evidence about the potential mechanisms of m6A regulators and their involvement in the development and progression of breast cancer and to clarify how this evidence could be used as a tool to predict prognosis of breast cancer patients.

Material and Methods

DATA COLLECTION AND ANALYSIS:

The R package getFirehoseData was used to acquire breast cancer data from the TCGA dataset (Breast Cancer, METABRIC, Nature 2012 & Nature Commun 2016). A total of 1904 cases with gene expression and clinical data were included in our analysis. Clinical information and RNA expression of breast cancer patients were downloaded and processed with R version 3.6.1.

Protein interaction data were downloaded from the Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) database [15] (https://string-db.org/) and then analyzed with Cytoscape (https://cytoscape.org/).

The Tumor Immune Estimation Resource (TIMER) [16] (https://cistrome.shinyapps.io/timer/) was adopted to obtainthe correlation among different m6A regulators after adjusting by tumor purity.

GENE SET ENRICHMENT ANALYSIS:

Gene set enrichment analysis (GSEA) (http://software.broadinstitute.org/gsea/index.jsp) was performed to determine whether a series of a previously defined set of genes were enriched in different biological groups [17]. Termswith P<0.05 and FDR<0.25 were identified. The hallmark gene set “h.all.v6.0.symbols.gmt” was used in this study.

LEAST ABSOLUTE SHRINKAGE AND SELECTION OPERATOR:

Least absolute shrinkage and selection operator (LASSO) is a regression analysis method that performs both variable selection and regularization to enhance the prediction accuracy and interpretability of the statistical model it produces. In our study, LASSO was utilized to evaluate the weight of each m6A regulator in the prognosis of patients [18]. The risk score of each patient was calculated with the following equation: risk score ∑i=1nCoef(i)*χ(i), where coef(i) was the value of each gene from LASSO analysis and x(i) was the expression value of each gene.

CONSENSUS CLUSTERING ANALYSIS:

The ConsensusClusterPlus package in R was used to perform consensus cluster analysis in our study. The algorithm started with subsetting a part of the characteristics and items in the data file, and each subset was divided into k groups according to k-means, the clustering algorithm used. This procedure was performed for a customized number of repetitions to yield a way to represent the consensus among several runs of the clustering and evaluate the solidity of the discovered clustering. Pairwise consensus values, the ratio of clustering runs in which 2 items are grouped together, were computed and kept in a consensus matrix for each k [19]. The clusters of patients with breast cancer were obtained according to the expression of m6A regulators. The prognosis of the different clusters was then compared.

ETHICS STATEMENT:

The study was approved by the institutional review board of our institution.

STATISTICS:

Chi-square tests or one-way ANOVA were applied to compare clinical and pathological data. Univariate and multivariate Cox regression analyses were conducted to evaluate the prognostic value of multiple parameters and characteristics. The Kaplan-Meier method was used to compare the OS of patients in different groups. All statistical analyses were conducted using R v3.6.1 (https://www.r-project.org/) or SPSS v19.0 (SPSS Inc., Chicago, IL, USA).

Results

RNA EXPRESSION AND CLINICOPATHOLOGICAL CHARACTERISTICS:

The m6A regulators used in our study were collected from published studies, and 20 molecules were included. The Nottingham prognostic index (NPI) was used to predict the outcomes of breast cancers following surgery [20]. First, we evaluated the potential relationships between the RNA expression of m6A regulators and NPI, cancer cellularity, and tumor staging. As shown in Figure 1A, about half of the 20 m6A regulators were statistically significantly related with tumor staging. In particular, EIF3 expression increased as tumor staging increased from stage I to stage III (Figure 1B). By contrast, the gradually decreased expression of RBM15B was observed as tumor staging increased (Figure 1C). Also, the link of NPI to m6A regulators was determined. Interestingly, most of the genes were associated with the NPI value (Figure 1D). EIF3 and IGF2BP3 expression were further invalidated in different NPI groups. Results showed a growing trend of EIF3 expression as NPI value became greater (Figure 1E), and IGF2BP3 showed the same expression pattern as EIF3 (Figure 1F). In addition, the role of m6A regulators in cellularity was assessed. Those results indicated that more than half of these genes were involved in cellularity distribution at a statistically significant level (Figure 1G). EIF3 (Figure 1H) and RBM15 (Figure 1I) showed similar expression patterns, in which both decreased as the cellularity proportion shifted from high to low. We also displayed the expression of m6A regulators across different tumor stages, cellularity, and NPI scores in Supplementary Table 1.

:

We also performed consensus clustering to yield an appropriate k value that could perfectly separate patients with breast cancer. k=4 was selected as an adequate choice, with the stability of clustering rising from k = 2 to 10, according to the expression likeness of the m6A regulators (Figure 2A, 2B). From the k value of 4, patients were divided into 4 clusters based on the consensus clustering results. Survival analysis was then conducted to compare whether the 4 groups of patients showed distinct 10-year OS. As shown in Figure 2C, patients in different clusters had significantly different outcomes, indicating that the clustering was powerful enough to predict the prognosis of patients with breast cancer.

:

Because the regulators were all involved in m6A modification, they were expected to be closely connected with each other. We first searched the STRING database to evaluate the genes’ protein interaction with each other. As shown in Figure 3A, the identified genes in STRING were closely correlated, indicating that these m6A regulators were functionally interrelated at the protein level. Then, the correlation between the regulators’ RNA in the TCGA database was further evaluated. The results indicated that most of the RNAs were related with each other, although with a small coefficient (Figure 3B), which was not consistent with protein interactions or with the results from a previously published study [11]. This might be due to the tumor heterogeneity across different types of cancers or tumor purity. Because the TCGA dataset was not adjusted by tumor purity, we analyzed data from the TIMER database, in which gene expression is adjusted by tumor purity, to thoroughly measure the link among the regulators. With the TIMER database, we found some stronger associations among these regulators than with raw data from the TCGA database. Four pairs of regulators that were closely related to each other are shown in Figure 3C.

:

Based on the above findings, we believed that these regulators contributed to patients’ outcomes. Cox regression analysis was then used to estimate the role of these genes in the OS of breast cancer patients. The results showed that most of the modification regulators were related with poor OS and 5 of them could predict OS (P<0.05) (Figure 4A). LASSO analysis was then used to additionally assess the effect of the 5 genes on prognosis, and the risk coefficient was acquired for each of the 5 regulators (Figure 4B). The results were consistent with the Cox regression analysis, showing that the LASSO risk coefficient was not zero. Subsequently, patients were divided into 2 groups according to the median risk score, which was calculated with the following equation: risk score ∑i=1nCoef(i)*χ(i). The survival analysis indicated that the risk score could effectively separate patients’ survival (P<0.0001) (Figure 4C). Multivariate analysis also demonstrated that the risk score was an independent prognosis prediction factor (hazard ratio 1.23, 95% confidence interval [CI] 1.08–1.39; P=0.002) (Table 1). These results showed the strong power of m6A regulators in predicting prognosis of breast cancer. In addition, other independent prognosis prediction factors in our study were patient age, radiotherapy, HER2 status, positive lymph node number, and tumor size (Table 1). GSEA analysis was then conducted to explore the potential underlying mechanisms, pathways, and cancer hallmarks. As shown in Figure 4D, 4 common cancer-related pathways, AKT1 signaling, cytokinesis, RNA transcription, and APC-MYC pathway, were enriched in patients with a high-risk score. Next, we explored whether the risk score influenced the effect of radiotherapy, chemotherapy, and hormone therapy on breast cancer. First, we identified patients who had received chemotherapy, radiotherapy, or hormone therapy. Then, survival analysis was performed in these patients. The results suggested that patients with a high risk score had relatively worse survival compared with those with a low risk score (Figure 4E), indicating that high-risk cases were not as sensitive to chemoradiotherapy or hormone therapies as were low-risk cases. Therefore, we recommend that patients with a high risk should receive closer follow-up to aid in early diagnosis and treatment of recurrence.

Discussion

Prior work has documented the role of m6A regulators in breast carcinogenesis, development, and drug resistance; for example, researchers found that FTO can promote the progression of breast tumors [21] and HNRNPA2/B1 is altered in endocrine-resistant breast cancer [22]. However, these studies have either included a single m6A regulator or were not focused on clinical relevance. In the present study, we studied the m6A regulators as a whole in breast cancer to test their relevance with clinicopathological features.

We identified that the RNA expression of m6A regulators is closely related to prognosis and other clinical features such as the NPI and cellularity of breast cancer. Here, we first divided breast cancer patients into 4 clusters by consensus clustering and found that the clusters significantly affected the prognosis of patients with breast cancer. Also, patients were further split into 2 subgroups according to the risk score value. Risk was correlated with prognosis and influenced cancer crucial pathways and biological processes of breast cancer.

Twenty m6A regulators have been included in this study. These regulators were searched from previously published data and have been validated to be involved in m6A modification in various cancer types besides breast cancer. According to the theory of tumor heterogeneity [23], which states that different tumor cells and tissues of distinct morphological and phenotypic profiles exist in all cancers, the function of some m6A regulators in breast cancer might be different from that in other cancers, where the role of the m6A regulators have been validated. This could lead to some inconsistencies between our results and those of other studies. For example, WTAP was reported to promote genes in various cancers, including hepatocellular cancer [24], bladder cancer [25], and ovarian cancer [26]. However, in the present study, WTAP was found to be a protective factor of OS. The inconsistency could be the result of the limited sample size and gap between clinical and experimental differences in the present study. CBLL1 is mainly studied in lung cancer and has been recognized to promote cancer progression in lung cancer [27]. It was also related to a favorable prognosis in our study, which warrants more studies to explore the potential effect of CBLL1 in breast cancer and the tumor heterogeneity between lung cancer and breast cancer. Our findings also indicated that dysregulation of certain m6A regulators resulted in alterations of RNAs in cancers and validated that the same regulators might exert divergent functions in different tumors.

In our study, the relationship between m6A regulators and clinicopathological characteristics was also analyzed. As a “reader”, EIF3 could be directly recruited by m6A methylation, which is followed by a cap-independent translation in the cellular stress response [28]. EIF3 is closely linked to the tumor stage, NPI, and cellularity of breast cancer. Also, EIF3 was associated with poor survival of breast cancer in the present study, which is consistent with previous research showing that EIF3 contributes to poor survival of patients with HER2 (+) breast cancer [29] and triple-negative breast cancer [30]. These pieces of evidence indicate that EIF3 has a pivotal role of in breast cancer. In addition, most of the regulators were found to be associated with breast cancer staging, NPI, and cellularity in our study, which is consistent with previous studies that found m6A is widely involved in cancer progression and prognosis [31].

It has been demonstrated that FTO can regulate cell migration and invasion in breast cancer via the miR-181b-3p/ARL5B signaling pathway, and high FTO expression is associated with advanced breast cancer staging [32]. A recent study also showed that tumor status and stage were relevant to the expression level of m6A RNA methylation regulators [33]. The m6A in peripheral blood RNA is closely associated with breast cancer stage [34]. The regulators promote the invasiveness of breast cancer, and therefore are correlated with an advanced tumor stage. NPI is calculated based on 3 pathological criteria: tumor size, the number of involved lymph nodes, and tumor grade [20]. The m6A regulators are involved in lymph node metastasis [32] and are therefore associated with NPI. However, low expression of some m6A regulators, such as METTL14 and ZC3H13, is associated with tumor progression, including NPI [35]. This could be explained by the potential dual function of m6A regulators in breast cancer. We also searched data on the relationship between m6A and cellularity in breast cancer but found little evidence about that relationship. This indicates that our findings are novel and might provide some clues for future studies.

Chemotherapy, radiation, and hormone therapy have been shown to greatly improve the prognosis of breast cancer. In addition to predicting OS, the risk score also performs well when predicting the response of patients with breast cancer after receiving the above treatment modalities. Therefore, patients with a high risk score should receive a higher dose of adjuvant therapies and close surveillance after treatment. Preclinical studies have indicated that the activation of PI3K/AKT/mTOR contributes to an acquired resistance to hormone therapy [36,37]. Data from other randomized trials also demonstrated that mTOR inhibition could improve hormone therapy resistance [38,39]. It has been reported that MYC is frequently increased in chemotherapy-resistant breast cancer patients and that MYC is able to maintain breast cancer stem cells in chemotherapy-resistant patients [40]. High-risk breast cancer patients are harboring the enrichment of the above pathways, which probably leads to radiation, chemotherapy, and hormone therapy resistance. These molecules and hallmarks have been targeted to develop new drugs and overcome resistance to adjuvant therapies. Further studies targeting m6A methylation regulators are warranted.

Solid tumor tissues consist of cancerous and noncancerous components, which contain epithelial, stromal, and endothelial cells [41]. Tumor purity, characterized as the proportion of cancerous cells in a solid tumor sample, is a key feature in cancer transcriptomics and metabolomics data analysis. Noncancerous components are viewed as prevailing pollutants in commixture and constitute a massive part of tumor masses. Studies have found that m6A modulators are closely related to each other, and our results from the STRING and TIMER databases also invalidated the findings. However, the correlation result of RNA data from the TCGA database without tumor purity adjusted was not consistent with the results of STRING or TIMER, an inconsistency that could have resulted from tumor purity. Typically, harvested surgical samples could be of a purity of lower than 70%, giving rise to systematic biases in genomic analysis [42]. Tumor purity is an important confounder and therefore should be considered to increase data accuracy and reliability.

Conclusions

In conclusion, our data demonstrated that m6A regulators were related to clinicopathological characteristics and showed prognostic significance in breast cancer. Future studies to evaluate the role of m6A regulators in other types of cancer are needed. The development of targeted therapeutics based on m6A regulators is also promising.