Identification of m6A-Associated Gene DST as a Prognostic and Immune-Associated Biomarker in Breast Cancer Patients.

BACKGROUND: N6-methyladenosine (m6A) is most common internal RNA modification in eukaryotic cells. Existing evidence shows that m6A is closely related to pathogenesis and progression in breast cancer (BRCA). Therefore, it is critical to investigate the key role of m6A target genes in BRCA.
METHODS: M6A target genes in BRCA are acquired using RMVar online database. The differentially expressed genes (DEGs) from three microarray datasets (GSE5764, GSE22358, GSE9014) is processed by GEO2R. Oncomine, GEPIA, UALCAN and TNMplot were applied to validate the expression of DST. Survival analyses were performed via DRUGSURV and Kaplan-Meier Plotter database. Univariable survival and multivariate Cox analysis were completed to assess the prognostic value of DST and receiver operating characteristic (ROC) curve was performed to evaluate the diagnostic value of DST. We also investigated the correlation between DST and cancer immune infiltration via using CIBERSORT, TIMER and TISIDB.
RESULTS: DST and COL11A1 were significantly expressed in both DEGs and m6A target genes set. COL11A1 show no significance on the patients' survival. However, high expression of DST was related to the favorable prognosis. Multivariate analysis revealed that the DST dysregulation is an independent prognostic factor and ROC indicated that the great diagnostic value of DST with AUC of 0.948. Subsequently, immunological analyses showed that DST was significantly associated with various immune infiltration cells, including NK cells, T helper cells and Mast cells. Furthermore, DST was also related with multiple immune checkpoints and chemokines, including LAG3, LMTK3 CD24, CXCL12, KDR and CX3CR1. These results indicated the potential roles of DST in the development of BRCA via altering the immune response.
CONCLUSION: DST can influence the development and progression of BRCA by altering the immune microenvironment.

Introduction

In women, breast cancer (BRCA) is the most common cancer and the first cause of cancer death surpassing lung cancer. BRCA is a heterogeneous disease composed of various biological subtypes.1,2 The administration of therapeutic targets based on molecular mechanism can improve clinical outcomes such as hormonal therapy and HER2-targeted therapy.3 Therefore, the discovery of novel therapeutic targets is critical for BRCA patients.

N6-methyladenosine (m6A) is a dynamic modification and its removal and installation have been demonstrated to be an important component of gene expression alteration.4,5 M6A can regulate RNA function and metabolism ranging from RNA splicing, nuclear export and translation to RNA degradation.6,7 The occurrence and disappearance of m6A is distributed to three regulators: “writers”; “erasers”; “readers”. To decide the fate of modified RNA, the “writers” generate m6A mark, “erasers” show demethylation activity and “readers” recognize m6A modification.8 Accumulating evidence indicates that m6A methylation plays a significant role in cancer pathogenesis via different mechanisms. YTHDF1 increases the translation of EIF3C via binding to m6A-modified EIF3C mRNA and promotes the tumorigenesis and metastasis of ovarian cancer.9 METTL3 interacts with the microprocessor protein DGCR8 and positively regulates the pri-miR221/222 process in an m6A-dependent manner, which plays an oncogenic role in bladder cancer.10 IGF2BP can augment the expression of SRF in m6A-dependent manner by blocking the miRNA-directed decay of the SRF mRNA, which promotes tumor cell invasion and growth.11 Thus, understanding the detailed roles of m6A in BRCA may provide us with more possibilities for treatments.

Cancer cells can functionally remodel tumor microenvironment (TME) via secreting diverse cytokines and chemokines.12 Immune cells are important constituents of the tumor stroma.13 Immune cells are composed of innate immune cells (neutrophils, macrophages, innate lymphoid cells, natural killer cells, dendritic cells, suppressor cells and myeloid-derived) and adaptive immune cells (T cells and B cells). Immune cells within TME play tumor-promoting and tumor-antagonizing properties of tumors.14 Accordingly, immunotherapy has become a promising field and therapeutic strategy.15–17

In this study, we identified a differentially expressed m6A target gene, DST, and analyzed the clinical significance of DST. Furthermore, the immunoregulatory role of DST was also explored in BRCA patients.

Methods

Data Acquisition

Three microarray datasets related with BRCA (GSE5764, GSE22358, GSE9014) were download from Gene Expression omnibus (GEO) database. The differentially expressed genes is obtained by GEO2R tool18 (). The adj. P value <0.05 and |logFC| > 1 were set as a cutoff criterion. M6A target genes in BRCA are acquired by RMVar online database19(). The co-DEGs among three dataset and m6A target genes were intersected by Venn diagram. Gene expression and clinical data were downloaded from the Cancer Genome Atlas (TCGA) website (). Clinical profiles, including age, sex, T stage, N stage, M stage, and follow-up data, were collected. Patients with incomplete date were deleted. Finally, a total of 1069 patients were included.

Expression and Clinical Significance Analysis

Gene expression was obtained from Oncomine20 () and validated by GEPIA21 (), UALCAN22 (), and TNMplot ().23 With the help of Kaplan-Meier Plotter24 () and DRUGSURV25 (), survival analysis was accomplished. The diagnostic value is assessed by the receiver operating characteristic (ROC) curve. Univariate and multivariate Cox regression analysis were performed to evaluate the prognostic power.

Correlation and Functional Enrichment Analysis

The co-expressed genes were acquired from LinkedOmics database26 (). Then, Gene Ontology biological process (GO_BP), Gene Ontology cellular component (GO_CC), Gene Ontology molecular function (GO_MF) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways were analyzed by the gene set enrichment analysis (GSEA) in the LinkInterpreter module.

Immunological Analysis

The correlations between gene expression and tumour-infiltrating immune cells from TCGA were performed. Tumor Immune Estimation Resource (TIMER) database27 () and TISIDB database28 () were applied for validation. Simultaneously, immune score of every patient is generated with ESTIMATE algorithm.

Statistical Analyses

R 4.1.0 and Statistical Package for Social Sciences 26.0 for Windows (SPSS Inc, Chicago, IL, United States) were used to perform statistical analyses. Clinicopathological parameters were compared using chi‐square or Fisher’s exact tests.

Univariate and multivariate analyses were conducted using the Cox-regression method. The pROC package29 was applied for ROC analysis. P value <0.05 was significant difference. The main tools involved in this study are summarized in Table 1.

Table 1

The Main Tools Involved in This Study

Database	Samples	URL	References
GEO2R	Tissue	https://www.ncbi.nlm.nih.gov/geo/geo2r/	[18]
RMVar	Tissue	http://www.rmvar.renlab.org/browse.html	[19]
Oncomine	Tissue	http://www.oncomine.org	[20]
GEPIA	Tissue	http://gepia.cancerpku.cn/	[21]
UALCAN	Tissue	http://ualcan.path.uab.edu/analysis.html	[22]
TNMplot	Tissue	https://www.tnmplot.com/	[23]
Kaplan-Meier Plotter	Tissue	http://kmplot.com/analysis/	[24]
DRUGSURV	Tissue	http://www.bioprofiling.de/GEO/DRUGSURV/index.html	[25]
LinkedOmics	Tissue	http://www.linkedomics.org/login.php	[26]
TIMER	Tissue	https://cistrome.shinyapps.io/timer/	[27]
TISIDB	Tissue	http://cis.hku.hk/TISIDB	[28]

Results

Identification of m6A Target Genes

The BRCA datasets of GSE5764, GSE22358 and GSE9014 were acquired to identify differentially expressed genes by GEO2R. 9308 DEGs was found in GSE9014, while 1538 DEGs in GSE22358 and 45 DEGs in GSE8977. As shown in Figure 1A, Venn plot showed 11 co-DEGS. We download m6A target geneset in BRCA from RMVar online database. Venn plot (Figure 1B) between co-DEGS and m6A target geneset depicted two m6A target genes (DST and COL11A1). Next, the prognostic values of DST and COL11A1 were explored by using survival data of BRCA patients from DRUGSURV databases. Survival curves of DST (Figure 1C and D) are statistically significant in both GSE11121 and GSE17705, but COL11A1 (Figure 1E and F) is not. Therefore, further investigations were focused on DST.

Figure 1

Identification of m6A target genes. (A) Venn plots of co-DEGs. (B) Venn plots of m6A target co-DEGs. (C–F) Kaplan–Meier curves of DST (C and D) and COL11A1 (E and F) in GSE11121 and GSE17705.

Validation of the Expression Level and Clinical Significance of DST

As depicted in Figure 2A, the expression of DST from Oncomine database was downregulated in various cancers, including BRCA. Subsequently, three online tools, UALCAN, TNMplot and GEPIA, were used to validate the downregulated DST (Figure 2B–D). Patients were divided into two groups according to the expression of DST. The association of DST expression with clinicopathological parameters is shown in Table 2. There was no significant difference in T, N, M stage and sex between the two groups. High expression of DST was significantly associated with the older age (p = 0.002). Moreover, in TCGA-BRCA, the patients with low level of DST displayed poor survival probability (p = 0.0054) and the hazard ratio (HR) is 0.63 (95% confidence interval [CI], 0.46–0.88) (Figure 2E). Univariate and multivariate Cox regression analyses (Table 3) indicate that DST is an independent prognostic factor (Figure 2F). The ROC curve showed a great diagnostic value of DST with AUC of 0.948 (95% CI: 0.925–0.970) (Figure 2G).

Figure 2

Validation of the expression and clinical significance of DST. (A–D) The expression of DST in breast cancer and normal validated by oncomine (A), UCLCAN (B), TNM plot (C), GEPIA (D). (E) The Kaplan–Meier curve of DST on overall survival of TCGA-BRCA. (F) Forest plot of risk factor. (G) The ROC curve of DST on diagnosing breast cancer. *P<0.05.

Table 2

The Characteristics of BRCA Patients from TCGA Database

Characteristic	Low Expression of DST (n=309)	High Expression of DST (n=760)	P value
Age
≤60	150	448	0.002
>60	159	312
Gender
Female	302	755	0.052
Male	7	5
T stage
T1-T2	257	639	0.616
T3-T4	52	118
M stage
M0	247	643	0.6743
M1	7	15
N stage
N0-N1	244	615	0.458
N2-N3	60	133
AJCC stage
StageI-II	215	530	0.560
StageIII-IV	80	180

Table 3

Univariate and Multivariate Analysis with OS in BRCA Patients

Clinicopathologic Parameter	Univariate Analysis			Multivariate Analysis
	HR	95% CI	P value	HR	95% CI	P value
Age (>60 vs ≤60)	1.97	1.42–2.7	<0.001	2.00	1.38–2.91	<0.001
Gender (male vs female)	0.83	0.12–5.95	0.854	0.53	0.07–3.87	0.532
T stage (T3-T4 vs T1-T2)	1.80	1.25–2.59	0.002	0.96	0.56–1.62	0.864
M stage (M1 vs M0)	4.90	2.93–8.21	<0.001	2.27	1.17–4.40	0.016
N stage (n1-n2 vs n0-n1)	2.30	1.58–3.37	<0.001	0.97	0.53–1.76	0.911
AJCC Stage (III-IVvsI-II)	2.67	1.92–3.72	<0.001	2.63	1.36–5.07	0.004
DST expression (high vs low)	0.65	0.47–0.91	0.12	0.68	0.46–0.99	0.045

Abbreviations: HR, hazard ratio; CI, confidence interval; OS, overall survival; AJCC, American Joint Committee on Cancer.

DST Co-Expression and Enrichment Analyses

To explore the biological function of DST in BRCA, the LinkedOmics tool was used to obtain the co-expression pattern of DST in TCGA-BRCA. 5785 genes showed a negative correlation with DST, while 9003 genes showed a positive correlation (Figure 3A). Figure 3B and C showed the top 50 genes positively and negatively related to DST, respectively. Enrichment analysis was completed to explore the potential biological functions of DST. The enrichment results of GO terms and KEGG pathway are shown in Figure 3D–G. Mitochondrial respiratory chain complex assembly, mitochondrial protein complex and structural constituent of ribosome were the top enriched terms in biological process, cellular component, and molecular functions, respectively. Proteasome, ribosome and oxidative phosphorylation were the top three enriched KEGG pathways.

Figure 3

The co-expression genes with DST from the LinkedOmics database. (A) Genes significantly related with DST by Pearson test. (B and C) Heat maps of top 50 genes positively and negatively correlated to DST. Red represents genes positively correlated and blue represents genes negatively associated. (D–G) GO terms for molecular functions (D), biological process (E), cellular component (F) and KEGG pathways (G) of DST by GSEA analyses.

The Role of DST in the Immune Microenvironment

The TME is increasingly recognized as an important factor in the development and progression of BRCA.30 Immune score of TCGA-BRCA, calculated using the ESTIMATE package,31 showed that patients with higher immune score had a better OS (Figure 4A). To investigate the role of DST in the immune microenvironment of BRCA, we further explored associations between DST expression and immune infiltration cells using ssGSEA algorithm. NK cells, T helper cells, Mast cells, eosinophils, Tgd, Dendritic cells (DC), macrophages, B cells, TFH and neutrophils were found to be positively correlated with DST in BRCA (Figure 4B). The top three cell types (NK cells, Mast cells, T helper cells) were validated using the TISIDB database (Figure 4C–E) and TIMER (Figure 4F–H) database, indicating that DST might influence the immune microenvironment in BRCA patients. Additionally, we also observed the negative associations between DST levels and multiple immune checkpoint molecules, including LAG3 (Figure 4I), LMTK3 (Figure 4J) and CD24 (Figure 4K).

Figure 4

The associations between immune infiltration and DST. (A) The Kaplan–Meier curve of immune score on overall survival of TCGA-BRCA. (B) Immune infiltration cells associated with DST using ssGESA algorithm. Validations of top three cells by TIMER database. (C–E) Validations of top three cells by TISIDB database, such as NK cells (C), T helper cells (D), Mast cells (E). (F–H) Validations of top three cells by TIMER database, such as NK cells (F), T helper cells (G), Mast cells (H). (I–K) The correlation between DST with LAG3 (I), LMTK3 (J) and CD24 (K).

Further investigations were focused on extending the recognition of immune roles of DST. TISIDB database was applied to analyze the associations between DST with immunostimulators, immunoinhibitors, chemokines, and receptors. As shown in (Figure 5A and B), multiple immunostimulators and immunoinhibitors showed associations with the levels of DST, with CXCL12 (Spearman:ρ=0.312) as the most significant immunostimulatory marker and KDR (Spearman:ρ=0.39) as the most significant immunoinhibitory marker. Simultaneously, Figure 6A and B showed that CXCL12 (Spearman:ρ=0.312) was the most significant chemokine and CX3CR1 (Spearman:ρ=0.344) was the most significant chemokine receptors. Next, we wanted to verify the associations between these molecules in above-mentioned GEO datasets. As shown in , CXCL12 was still significantly related with DST in GSE22538 (Spearman:ρ=4.94E-13) and GSE9014 (Spearman:ρ=0.038). Interestingly, KDR was also significantly related with DST in GSE22538 (Spearman:ρ=1.93E-9) and GSE9014 (Spearman:ρ=3.72E-3). However, no significant association could be found in GSE5764, which needs to be further clarified. These results together suggested that DST might play a critical role in the progression and development of patients with BRCA by regulating the immune response.

Figure 5

Associations of the DST expression with immunomodulators from TISIDB database. (A) Correlations between immunostimulators and DST with four most significant molecules. (B) Correlations between immunoinhibitors and DST with four most significant molecules.

Figure 6

Associations of the DST expression with chemokines and receptors from TISIDB database. (A) Correlations between chemokines and DST with four most significant molecules. (B) Correlations between receptors and DST with four most significant molecules.

Discussion

In this study, we found a m6A target gene, DST, plays an important role in the pathogenesis and development of BRCA, acting as a potential tumor suppressor. Further investigations showed that the DST was involved with the immune microenvironment of BRCA, indicating that DST might make an impact on BRCA by affecting the immune microenvironment. These results suggested that DST might be reckoned as a biomarker or therapeutic target in patients with BRCA.

There are several studies introducing the importance of m6A in BRCA. M6A is the most abundant internal modification of RNA and is enriched in various cancers PMID: 32355831. M6A modification is installed by the m6A methyltransferases, removed by the demethylases and are recognized by m6A-binding proteins, which alter gene expression via regulating RNA metabolism, including translation, splicing, export and degradation.32 M6A target genes alteration in BRCA can modulate cellular processes including cell growth, metastasis, invasion and apoptosis.33–35

We obtained a list of m6A target genes from RMVar online database. These genes were intersected with DEGs from three GEO datasets. We acquired DST and COL11A1. According to survival data from DRUGSURV database, DST was selected as a candidate. Subsequently, we validated the levels of DST in BRCA with Oncomine and three online databases. Additionally, we also found that DST was downregulated in many cancers, but little mechanism was explored in human cancer.

We explored the clinical significance of DST. DST was an independent prognostic factor via multivariate Cox regression analyses and the ROC curve showed a great diagnostic value of DST with an AUC of 0.948. Higher expression of DST was accompanied by better survival outcomes, but there was no significant correlation between T, N, M stage and sex.

The above results reveal that DST might be a critical biomarker in BRCA. Recently, increasing evidence has shown that the degree of heterogeneity of the TEM in BRCA could be recognized as an important element for cancer development and progression, as well as potential therapeutic targets.36–38 Therefore, we investigated the relevance between DST and tumor immune microenvironment. DST was positively correlated with multiple immune cells and negatively correlated with many immune checkpoints. The top three cell types (NK cells, Mast cells, T helper cells) were validated using the TISIDB and TIMER databases. NK cells function as the main effector cells toward cancer in innate immunity,39 and inhibition of NK cells leads to immune escape and multi-step metastatic process of BRCA.40 Mast cells play a dual role in both anti-tumoral and pro-tumoral functions in BRCA,41 which highlights the importance of personalized treatments. T helper cells are critical regulators of immune with great therapeutic and prognostic value in BRCA.42,43 In this study, we found DST might function via altering the immune microenvironment of BRCA, and this is the first study reporting the correlation between DST and immune microenvironment in BRCA.

There were several limitations in this study. These results were mainly derived from the bioinformatics analysis but lacked adequate experimental validation. More experimental and clinical researches were needed to further validate the biological role of DST in breast cancer.

Conclusion

In conclusion, DST may alter the development and progression of BRCA via influencing the immune microenvironment and possess diagnostic and therapeutic values to some extent.