ATP2C2 Has Potential to Define Tumor Microenvironment in Breast Cancer

Tumor microenvironment (TME) is vital for the occurrence and development of breast cancer (BRCA). However, it remains challenging to understand the dynamic modulation of the stromal and immune components comprehensively in TME. Herein, we used ESTIMATE and CIBERSORT algorithm to estimate the number of stromal and immune components and the abundance of tumor-infiltrating immune cells (TICs) in 582 BRCA cases from gene expression omnibus (GEO) database. We employed three regression models including univariable Cox proportion, LASSO regression model and multivariate Cox regression, and identified 7 immune-specific genes related to BRCA survival. Of 7 genes, ATPase Secretory Pathway Ca²⁺ Transporting 2 (ATP2C2) attracts our attention for significantly predicting prognosis of BRCA patients. Further analysis indicated that ATP2C2 expression was closely related to the clinicopathological features (age, T- and N-staging) and negatively correlated with patients’ survival in BRCA. Gene Set Enrichment Analysis (GSEA) was performed to reveal pathway enrichment between ATP2C2^high and ATP2C2^low groups. The low ATP2C2 expression groups’ genes were mainly enriched for immune-related activities, while those in the ATP2C2 high-expression group were largely enriched in metabolic-related pathways. Notably, Pearson’s correlation analysis identified that ATP2C2 expression was positively correlated with T follicular helper (Tfh) cells, and negatively correlated with gamma delta (γδ) T cell, suggesting that ATP2C2 might be accountable for the maintenance of immune-dominant status for TME. To sum up, this study comprehensively analyzed the TME and shed light on prognostic immune-related biomarkers for BRCA. In particular, ATP2C2 might be helpful for predicting the prognosis of BRCA patients, which provided an extra insight for BRCA treatment.

Introduction

Breast cancer (BRCA) with high aggressiveness is one of the most common malignant tumors that seriously endanger women’s health (1). However, the mechanism of its occurrence and development is not yet fully understood. Therefore, in-depth study of the molecular mechanisms of BRCA development and finding effective therapeutic targets are of great significance for improving the prognosis of BRCA patients.

The “seed-soil” theory of tumors believes that the growth of tumor cells requires the surrounding normal cells and extracellular matrix to provide a permissive environment, that is, the tumor microenvironment (TME) (2). TME is the local environment for tumor cells to survive. In addition to tumor cells, it also contains resident stromal cells and recruited immune cells. These components are crucial in the occurrence, development, and immune evasion of tumors (3). Although stromal cells can promote tumor angiogenesis, cancer cell proliferation, invasion and metastasis, the detailed molecular mechanism to promote cancer progression has not been fully elucidated. Meanwhile, increasing attention has been paid to the influence of the immune cells in TME on tumor development. It was revealed that high levels of immune cell infiltration into cancer tissues were correlated with favorable outcomes, suggesting that valuing TME heterogeneity and remodeling the immune microenvironment may hold promise for cancer treatment (4). Immune cells in the TME can affect the host immune response by secreting chemokines, cytokines and other factors that directly or indirectly suppress or support tumor progression (5). Therefore, to better understand the immune status of the TME and investigating the distribution pattern and function of immune cells are essential for improving the effectiveness of cancer immunotherapy.

In the present study, we used ESTIMATE and CIBERSORT algorithm to quantify the tumor-infiltrating immune cell (TIC) proportion and the ratio of immune and stromal components of BRCA samples from the gene expression omnibus (GEO) database and identified a novel predictive biomarker, ATPase Secretory Pathway Ca²⁺ Transporting 2 (ATP2C2). Gene dysregulation of ATP2C2 has been linked to breast cancers: ATP2C2 (also known as SPCA2) promotes tumor growth by increasing Ca²⁺ entry through activation of the Orai1 calcium channel (6). Recently, the literature has demonstrated a novel association among ATP2C2, Kv10.1 and Orai1 involved in mediating transduction signals from TME to the BRCA cells, suggesting that ATP2C2 may play an important role in TME (7). Here we set out to compare differentially expressed genes (DEGs) produced by comparison between stromal components and immune components in BRCA cases, revealing that the ATP2C2 may be a potential indicator for the alteration of TME status in BRCA.

Materials and Methods

BRCA Datasets and Samples

mRNA expression data for breast cancer (BRCA) were procured from GEO (https://www.ncbi.nlm.nih.gov/geo/). The keywords “Breast cancer” and “Survival” were used for retrieval. Finally, five gene expression microarray datasets were chosen as a training set, including (GSE42568, n=121), (GSE88770, n=117), (GSE16446, n=120), (GSE37751, n=108), and (GSE7390, n=198). 582 samples with the survival status and survival time were downloaded for DEG analysis from the five datasets. In addition, the GSE20711 cohort (n = 88) were used as a validation set.

Harmonized RNA sequencing data and related clinical information for BRCA were downloaded for clinical correlation and GSEA analysis from TCGA (https://portal.gdc.cancer.gov/, up to July 30, 2020), which included 1222 samples, 113 normal tissue samples, and 1109 tumor samples.

ESTIMATE

The stromal, immune and ESTIMATE scores were outputted by the R package “estimate” (8).

DEGs Identification Based on StromalScore and ImmuneScore

Based on the comparison with the median scores of StromalScore and ImmuneScore, 582 BRCA cases were divided into high or low score groups. Limma R package was applied to determine DEGs between high-scoring samples and low-scoring samples. Genes with p < 0.05 and |log2 FC (fold-change)| > 1 were screened DEGs. A flowchart for constructing DEG signature is displayed in Figure 1.

FIGURE 1

Figure 1 Experimental technical roadmap.

Heatmaps

Heatmaps were plotted with the R package pheatmap (9).

Functional Enrichment and Pathway Analysis

To understand the potential function of the common DEGs that stemmed from the immune scores and stromal scores, Gene Ontology (GO) and the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis were performed by “clusterProfiler” package in R. GO and KEGG terms with a p- and q-value of both <0.05 were considered significantly enriched.

Construction and Validation of the Prognostic Prediction Models

A univariate Cox model was conducted to screen potential prognostic immune-related genes (IRGs) and TICs. A LASSO Cox regression model was further applied to narrow the range of prognostic IRGs and TICs. Next, a multivariate Cox regression model was utilized to select IRGs and TICs most closely related to survival, and those IRGs and TICs were used to construct two risk models. Two formulas for the IRGs risk score and TICs risk score was established to predict patient survival: IRGs risk score = (−0.212 × expression level of ADRB2) + (0.208 × expression level of ATP2C2) + (−0.287 × expression level of CELF2) + (−0.129 × expression level of CXCL12) + (0.257 × expression level of LGMN) + (−0.204 × LIPA) + (0.214 × expression level of SLCO2B1), TICs risk score = (3.359 × abundance of naïve B cells) + (4.957 × abundance of T follicular helper (Tfh) cells) + (−4.349 × abundance of gamma delta (γδ) T cells) + (−3.966 × abundance of resting mast cells). The individual IRGs risk score and TICs risk score in both training set and validation set were calculated according to the formulas accordingly.

The Kaplan-Meier method and the ROC curves were used to analyze OS and to assess signature’s sensitivity and specificity of the models.

Construction and Validation of Nomograms

Based on seven IRGs and four TICs, we separately plotted the nomograms to predict the probability of 1-, 2-, and 3-OS of BRCA patients. Validation of the nomograms was evaluated by the discrimination and calibration.

Gene Set Enrichment Analysis (GSEA)

GSEA was performed between high- and low-ATP2C2 expression (10). The number of random sample permutations was set at 1000, and NOM p-value < 0.05, FDR q-value < 0.25, and | NES | > 1 were set as the significance threshold.

Estimation of TIC Types

CIBERSORT algorithm was utilized to evaluate the relative abundance of 22 types of TICs. After the quality filtering (p-value < 0.05), 545 BRCA samples were selected for following analysis.

Tumor Immune Dysfunction and Exclusion (TIDE) Analysis

TIDE is a computational framework construct to predict immune checkpoint inhibitors (ICIs) response (11).

Results

DEGs Based on ImmuneScore and StromalScore Were Mainly Presented as the Enrichment of IRGs

To explore the accurate alterations of gene profile in TME regarding immune and stromal components, a comparative analysis of high- and low-scoring samples was performed. A total 1091 DEGs including 692 upregulated and 399 downregulated genes were obtained from ImmuneScore (Figure 2A). Similarly, a total of 1261 DEGs including 852 upregulated and 409 downregulated genes were obtained from StromalScore (Figure 2B). By intersecting these DEGs from ImmuneScore and StromalScore, a total of 246 up-regulated genes and 47 down-regulated genes (Figures 2C, D). These DEGs were probably key factors for the status of TME. GO enrichment analysis showed that the functions of the DEGs were predominantly associated with the immune response, such as leukocyte proliferation and leukocyte migration (Figure 2E; Figure S1A). The KEGG enrichment analysis also indicated that the DEGs were involved in hematopoietic cell lineage, Th1 and Th2 cell differentiation, and cytokine–cytokine receptor interaction (Figure 2F; Figure S1B). Therefore, immune-related activities tended to represent the main functions of DEGs, implying that the involvement of immune components was a predominant feature of TME of BRCA.

FIGURE 2

Figure 2 Comparison of gene expression profile with stromal and immune scores of BRCA. (A, B) The heatmap presenting the top 20 upregulated and downregulated DEGs in ImmuneScore and StromalScore. DEGs were detected by Wilcoxon rank sum test (q < 0.05 & |log2FC| > 0.5). (C, D) Venn diagram analysis of aberrantly expressed genes based on stromal and immune scores. (E, F) DEGs-related biological functions and pathways in BRCA, terms with p and q < 0.05 were considered to be significantly enriched.

Screening of Prognosis-Specific IRGs and Construction of Prognosis Prediction Model

293 differentially expressed IRGs were further subjected to the univariate Cox regression model; we identified 67 IRGs significantly related with OS (Table S1). Then, these 67 IRGs were used in the LASSO regression for feature selection. A set of seventeen genes (ADD3, ADRB2, ALDH1A1, ATP2C2, CCR2, CELF2, CRTAM, CXCL12, EPN3, HOXC8, IL2RB, ITM2A, KLRG1, LGMN, LIPA, SLCO2B1, and TMEM243) and their coefficients were computed (Figures 3A, B). Then, the Akaike information criterion (AIC) was utilized to screen significant prognostic IRGs in multivariate Cox regression models (Figure 3C).

FIGURE 3

Figure 3 Signature-based risk score is a promising marker of survival in BRCA patients. (A, B) LASSO Cox analysis identified seven genes most correlated with OS, and 10-round cross validation was performed to prevent overfitting. (C) Multivariate Cox analysis of 7 immune-related hub genes. (D, G) Risk score distribution in training set (D) and validation set (G). (E, H) Survival overview in training set (E) and validation set (H). (F, I) Expression profile of 7 immune genes in training set (F) and validation set (I). (J, K) Kaplan–Meier estimates of OS according to the seven-immune-related gene signature in training set (J) and validation set (K). The differences between the two curves were evaluated by the two-side log-rank test. (L, M) The ROC curve analysis of the seven-immune-related gene signature for predicting OS in training set (L) and validation set (M).

The risk score of each BRCA patient was calculated, and the patients were classified into the high-risk (n = 291) or low-risk (n = 291) group by the median cut-off value (Figure 3D). Remarkably, the number of deaths was significantly higher in the high-risk group (Figure 3E). A heatmap revealed that patients in the high-risk group tended to have increased LGMN, ATP2C2, and SLCO2B1 expression levels, as well as decreased expression levels of CXCL12, LIPA, ADRB2 and CELF2 (Figure 3F). The Kaplan–Meier analysis indicated that low-risk patients had a better OS than high-risk patients (P < 0.001) (Figure 3J). The accuracy of the prognostic model was displayed in the ROC curve. The AUCs under ROC curve for predicting OS in the first, third, and fifth year were 0.831, 0.731, and 0.723, respectively (Figure 3L). To demonstrate the robustness of the prognostic signature, the predictive ability was assessed in a GSE20711 cohort. The risk score for each patient was calculated based on the same formula, and the best threshold was chosen as the cut-off for patients stratified as high- and low-risk. The risk score profiles and gene expression are displayed in Figures 3G–I. Patients in high-risk group had a shorter OS than patients in the low-risk group in the GSE20711 cohort, which was consistent with the results of the training set (Figure 3K). The AUC of the signature for OS was 0.92, 0.668 and 0.693 at 1, 3 and 5 years, respectively (Figure 3M).

Based on the multivariate Cox analysis results, seven IRGs were integrated in the nomogram to predict the OS of BRCA patients (Figures 4A, B). The calibration curves revealed acceptable accuracy. These results indicated that the IRGs risk model nomogram had very appropriate calibration.

FIGURE 4

Figure 4 Nomogram predicting OS for BRCA patients based on risk score. (A) Nomogram model for predicting the probability of 1-, 2-, and 3-year OS of BRCA. Points are assigned for seven features. The sum of these points is located on the total points axis. The total points on the bottom scales correspond to the predicted 1-, 2-, and 3-year survival. (B) Calibration plots of the nomogram for predicting the probability of OS at 1, 2, and 3 years. The X-axis represents nomogram-predicted survival, and the Y-axis represents actual survival.

ATP2C2 Is Highly Expressed in BRCA and Correlates With Unfavorable Prognosis

Ion channels have a key role in mediating TME signal (12). Previous studies had shown that ATP2C2 can regulate the localization and the activity of Kv10.1 and Orai1 channels, mediating transduction signals from TME to the BRCA cells (7). However, the clinical role of the ATP2C2 has not yet been determined. Here we divided all BRCA samples into ATP2C2 high and low expression groups based on ATP2C2 median expression. The Kaplan–Meier analysis elucidated that BRCA patients with ATP2C2 low expression had longer survival than those with ATP2C2 high expression (Figure 5A). In addition, we further found that ATP2C2 levels are associated with poor prognosis of various cancers, including thyroid carcinoma (THCA), head-neck squamous cell carcinoma (HNSC), kidney renal clear cell carcinoma (KIRC), lung squamous cell carcinoma (LUSC), and esophageal squamous cell carcinoma (ESCC) on a web platform Kaplan–Meier Plotter (http://kmplot.com/analysis/) (Figures S2A–E). After that, clinical pathologic features and ATP2C2 expression were downloaded from the TCGA database. The expression of ATP2C2 in the BRCA samples was significantly higher than that in normal breast tissues or paired breast non-tumor by Wilcoxon rank sum test (Figures 5B, C). We also observed a similar up-regulation of ATP2C2 expression in BRCA in the HPA dataset (https://www.proteinatlas.org/) (Figures S3A–F). Collectively, these outcomes clearly demonstrated that high ATP2C2 expression in TME was correlated with worse prognosis of BRCA patients. Remarkably, the levels of ATP2C2 were increased along with the progression of age, T- and N- staging (Figures 5D–I).

FIGURE 5

Figure 5 The expression of ATP2C2 and correlation with OS and clinicopathological features of BRCA patients revealed by bioinformatic analysis. (A) Survival outcome for BRCA patients with different ATP2C2 expression. The differences between the two curves were detected by the two-side log-rank test. (B) ATP2C2 mRNA is highly expressed in BRCA tissues. Statistical significance was evaluated using Wilcoxon-Mann-Whitney test. (C) ATP2C2 mRNA levels in the paired BRCA tissues were evaluated. P-values are based on the Wilcoxon Test. (D–I) The correlation of ATP2C2 expression with clinicopathological features. Statistical significance was evaluated using Wilcoxon rank sum test.

ATP2C2 Had Potential to Be an Indicator for TME Modulation

Given the expressions of ATP2C2 were negatively associated with the OS, T- and N- staging of patients with BRCA, GSEA was performed to explore the gene sets enriched in different ATP2C2 subgroups. As shown in Figure 6A, ATP2C2 low-expression group’ genes were significantly enriched in immune-related activities, such as cytokine–cytokine receptor interaction, hematopoietic cell lineage and primary immunodeficiency, while those in ATP2C2 high-expression group were mainly enriched in metabolic pathways, including citrate cycle TCA cycle, amino sugar and nucleotide sugar metabolism, and fructose and mannose metabolism (Figure 6B). Together, these findings implied that ATP2C2 might be a potential indicator for the status of TME.

FIGURE 6

Figure 6 GSEA revealed biological pathways correlated with ATP2C2 in the cohort from TCGA. (A) Gene sets enriched in ATP2C2 high-expression group. (B) Gene sets enriched in ATP2C2 low-expression group. NOM p < 0.05, FDR q < 0.25, and |NES| > 1 are set as the significance threshold.

ATP2C2 Predicts the Infiltration of Immune Cells Into BRCA Microenvironment

To further explore the indicative roles of ATP2C2 on TME, we used the CIBERSORT algorithm to detect the proportions of 22 kinds of immune cells in the BRCA microenvironment (Figure 7A). Besides, the results of the Wilcoxon-Mann-Whitney test showed that the fractions of the naïve B cells, γδ T cells, activated memory CD4⁺ T cells, naïve CD4⁺ T cells, monocytes, M0 macrophages, M2 macrophages, and activated dendritic cells in ATP2C2 high-expression group was relatively less than that in ATP2C2 low-expression group, and regulatory T cells (Tregs), activated NK cells, and M1 macrophages were relatively greater in ATP2C2 high-expression group (Figure 7B).

FIGURE 7

Figure 7 Analysis of immune cell infiltration. (A) Fraction data of the 22 types of immune cells. (B) Differential immune cell type expression was observed between the high and low-ATP2C2 groups. Significant statistical differences between the two subgroups were assessed using the Wilcoxon test.

ATP2C2 Was Correlated With Distribution Pattern of T Cell Subsets

All TICs were integrated into a univariate Cox regression model (Table 1). The significative TICs, screened by univariate Cox regression analysis were subjected to LASSO Cox regression analysis. The results of the Lasso regression suggested that the model was not overfitting (Figure S4A, B). After that, multivariate Cox regression analysis revealed that the model composed of naïve B cells, Tfh cells, γδ T cells, and resting mast cells had the smallest AIC (Figure S4C). According to the median cut-off value of TICs risk score, BRCA patients were classified into high- and low-risk groups. The distribution of risk score, survival status and expression profile of the four TICs of each patient are presented in Figures S4D–F (training set) and Figures S4G–I (validation set). The Kaplan–Meier survival analysis in the two datasets revealed significantly worse prognosis in the high-risk group (Figures S4J, K). Next, the prognostic accuracy of the TICs risk score was examined in the training set and validation set by using time‐dependent ROC curves analysis (Figures S4L, M).

TABLE 1

Table 1 Univariate analysis showing associations between different immune cell subsets and OS in BRCA. Unadjusted HRs are shown with 95 percent confidence intervals.

To further improve the accuracy of the prognostic, the nomogram based on the multivariate analysis was constructed (Figure S5A). In addition, the calibration curve and the ROC demonstrated good discrimination and concordance (Figure S5B).

Pearson analysis was applied to demonstrate the co-expression patterns among diversified immune cells (Figure 8A). Likewise, correlation relationship between immune cells and key genes was further analyzed and illustrated (Figure 8B). Here we emphatically analyzed the correlation between T cells and ATP2C2. Among them, Tfh cells were positively correlated with ATP2C2 expression (R = 0.12, P = 0.0059) (Figure 8C); γδ T cells was negatively correlated with ATP2C2 expression (R = -0.18, P = 1.9e-05) (Figure 8D). These outcomes further verified that the expressions of ATP2C2 influenced the immune activity of TME.

FIGURE 8

Figure 8 Correlation of TICs proportion with ATP2C2 expression. (A, B) Correlation analysis of different TICs and the relationships between different TICs and DEGs in tumor tissues of BRCA. Pearson’s correlation coefficient (r) was used for the significance test. (C, D) Correlation between Tfh cells, γδ T cells and ATP2C2 expression. Correlation test is conducted by the Pearson coefficient. p-value < 0.01 is the significance threshold.

ATP2C2 Promotes Immune Evasion and Resistance to ICIs

Accumulating evidence suggests that the level of cytotoxic T lymphocytes (CTL) is correlated with a better prognosis of patients. Through the analysis of the TIDE algorithm, we found that in BRCA patients with low ATP2C2 expression levels, high CTL levels indicate a better prognosis, while the above phenomenon is not observed in BRCA patients with high ATP2C2 expression levels (http://tide.dfci.harvard.edu/) (Figure 9). Similar results are observed in acute myeloid leukemia (AML), lung adenocarcinoma (LUAD), and ovarian serous cystadenocarcinoma (OV) patients (Figures S6A–C). In addition, we further used the TIDE algorithm to predict the efficacy of patients with ICIs. We found that the expression level of ATP2C2 is related to the enhanced efficacy of ICIs treatment in melanoma and KIRC (Figures S6D, E). These results suggested that ATP2C2 may be appropriate candidates for immunotherapy, especially ICIs.

FIGURE 9

Figure 9 Validation of ATP2C2 as a regulator of tumor immune escape. The association between the CTL level and OS for BRCA with different ATP2C2 levels.

Discussion

Herein, we sought to mine TME-related genes that contributed to the classification of TNM stages and the OS in BRCA patients from the GEO and TCGA-BRCA datasets. ATP2C2 was determined to be involved in immune activities. Interestingly, we found that ATP2C2 might be an indicator for the status of TME in BRCA patients through bioinformatics analysis.

TME is the key regulator of carcinogenesis and consists of tumor cells, stromal cells, and immune cells (13, 14). TICs are closely associated with tumorigenesis, angiogenesis and the growth and metastatic potential of tumor cells, which could alternately modulate the pattern of immune cells (15–18). Now, immunotherapy has been widely used to treat a broad spectrum of cancers including BRCA (19–22). However, not all patients can benefit from it. Thus, it is of crucial significance to improve treatment efficacy of BRCA and to uncover strong prognostic biomarkers for BRCA. Here, 7 prognosis-specific IRGs were identified by a series of bioinformatics analysis. The risk score calculated using the 7 IRGs was an independent prognostic factor for BRCA. Among all 7 prognosis-specific IRGs, four (e.g., ADRB2, CXCL12, LGMN, LIPA) have been reported to be involved in the immune microenvironment-associated pathogenesis of BRCA, implying that our bioinformatics analysis using GEO cohorts has prognostic value. The remaining three genes including ATP2C2, CELF2, and SLCO2B1 have not been formerly reported to be associated with BRCA patients’ prognosis and could serve as novel potential biomarkers for BRCA. Here we gave special attention to ATP2C2 and then embarked from the transcriptomic analysis of BRCA in TCGA database, which revealed that increased ATP2C2 expression was significantly correlated with the advanced clinicopathological characteristics (age, T‐ and N‐staging) and unfavorable prognosis of BRCA patients.

ATP2C2 is mainly expressed in salivary glands, gastrointestinal and respiratory tracts, and mammary gland, and participates in Ca²⁺ transport in secretion (23). ATP2C2 has been shown to be involved in many important biological functions, such as the regulation of calcium homeostasis, and modulation of phonological short-term memory in language impairment (24, 25). In addition, it was shown that ATP2C2 helps colon cancer cells adapt to hypoxia, prevents cancer cells death, increases proliferation capacity and promotes tumor growth (26). ATP2C2 can serve as an independent prognostic factor and has better prediction for the survival of thyroid cancer patients (27). Therefore, ATP2C2 seems to play an antitumor role in BRCA. Here, relationship between ATP2C2 expression and TME were further studied. GSEA results revealed that high ATP2C2 expression was associated with metabolic-related signaling pathways, such as fructose and mannose metabolism, citrate cycle TCA cycle, and amino sugar and nucleotide sugar metabolism. In the ATP2C2 low-expression group, immune pathways including cytokine–cytokine receptor interaction, hematopoietic cell lineage and primary immunodeficiency were significantly enriched. Previous studies have validated that hypoxia-reprogrammed TCA cycle promotes breast tumorigenic cells growth (28). The expression of ATP2C2 is highly related to the TCA cycle, which may affect the prognosis of BRCA patients through the TCA cycle. Changes in cytosolic Ca²⁺ trigger events critical for tumorigenesis, such as cellular motility, proliferation, and apoptosis (6). Suppression of ATP2C2 attenuates basal intracellular Ca²⁺ levels and breast tumorigenicity (6). ATP2C2 can regulate Ca²⁺ metabolism and affect the prognosis of BRCA patients. Consistent with this, our analysis verified that BRCA patients with higher ATP2C2 expression had shorter OS time. In another study, the authors concluded that ATP2C2 antagonizes epithelial mesenchymal transition and suppresses BRCA cell migration and tumor metastasis (29). This study is inconsistent with our present observations and other studies (6, 7, 30). Perhaps it emphasizes that ATP2C2 can regulate epithelial to mesenchymal transition and BRCA metastasis, rather than studying the effect of ATP2C2 on the growth of BRCA. Together, our results suggested that ATP2C2 could participate in the conversion of TME from immune-dominant to metabolic-dominant status.

Considering the importance of immune cell infiltration in tumors (31–33), CIBERSORT was further applied to evaluate the abundance ratios of 22 types of immune cells in each BRCA specimen from GEO. Mounting evidence suggests that the interaction between the tumor and the microenvironment plays an important role in the progression of BRCA and immunotherapeutic efficacy. Therefore, we evaluated the potential of ATP2C2 to reflect the infiltration of immune cells. We found out the different proportions of numerous immune cells in different ATP2C2 subgroups. ATP2C2 low-expression group had significantly higher proportions of the naïve B cells, naïve CD4⁺ T cells, activated memory CD4⁺ T cells, monocytes, γδ T cells, M0 macrophages, M2 macrophages, and activated dendritic cells and significantly lower proportions of M1 macrophages, activated NK cells, and Tregs than ATP2C2 high-expression group. Tumor-infiltrating macrophages play a crucial role in tumor behavior and clinical outcome. Classically (M1) and alternatively activated (M2) macrophages exhibit different phenotypes and functions. M1 macrophages secrete cytokines, including TNF-α, IL-6, and IL-12 (34), killing tumor cells in the TME (35). M2 macrophages can secrete anti-inflammatory factors such as TGF-β and IL-10, and promote tumor growth and metastasis (36). We found that more M1 macrophages infiltrated in ATP2C2 high-expression group compared with ATP2C2 low-expression group, implying that tumor infiltrated macrophages exert immune response functions and exhibit anti-tumor effects.

With further use of LASSO Cox regression models, as a statistical method for screening cell variables to establish the TICs risk model, the predictive accuracy could be improved significantly. Moreover, a nomogram was made based on four immune infiltrating cells to establish a more accurate prognostic prediction model for BRCA. Finally, the correlation analysis revealed that there was a significantly positive correlation between Tfh cells and ATP2C2 expression, and a negative correlation between γδ T cells and ATP2C2 expression. The main function of Tfh cells is to help B cells produce antibodies. As early as 2013, Tfh cells were first found in BRCA (37). However, the role of these Tfh cells in BRCA is not clear. It is discovered that the original Tfh cells in BRCA secrete CXCL13 to promote the accumulation of Tregs in the tumor, inhibit the body’s anti-tumor immunity, and ultimately promote the development of BRCA (38). One small component of this microenvironment in humans is γδT lymphocytes, which display both innate and adaptive functions (39). Early studies have shown that γδ T lymphocytes in the BRCA cell lines have anti-tumor activity (40). Previous studies have shown that γδ T lymphocytes inhibit angiogenic signalling pathways associated with AKT and ERK, and increase apoptosis (41, 42). Interestingly, ATP2C2 promotes BRCA cell-cycle progression and cell proliferation via RAS-ERK pathway (6). Therefore, ATP2C2 might be responsible for the preservation of immune-active status in TME.

Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

Author Contributions

HL and GR conceived and designed the study. JZL, YXW, and YSW searched and selected studies, analyzed the data, wrote and revised the paper. JL and JS prepared the figures and tables. All authors contributed to the article and approved the submitted version.

Funding

This study was supported by National Natural Science Foundation of China (#31420103915, #81472475).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

Thanks to the GEO for providing the original data.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu.2021.657950/full#supplementary-material

References

1. Krug K, Jaehnig EJ, Satpathy S, Blumenberg L, Karpova A, Anurag M, et al. Proteogenomic Landscape of Breast Cancer Tumorigenesis and Targeted Therapy. Cell (2020) 183(5):1436–56 e31. doi: 10.1016/j.cell.2020.10.036

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Kaymak I, Williams KS, Cantor JR, Jones RG. Immunometabolic Interplay in the Tumor Microenvironment. Cancer Cell (2020) 39(1):28–37. doi: 10.1016/j.ccell.2020.09.004

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Toyama T, Iwase H, Yamashita H, Hara Y, Omoto Y, Sugiura H, et al. Reduced expression of the Syk gene is correlated with poor prognosis in human breast cancer. Cancer Lett (2003) 189(1):97–102. doi: 10.1016/s0304-3835(02)00463-9

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Janes PW, Vail ME, Ernst M, Scott AM. Eph receptors in the immune-suppressive tumor microenvironment. Cancer Res (2020) 81(4):801–5. doi: 10.1158/0008-5472.CAN-20-3047

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Akhand SS, Liu Z, Purdy SC, Abdullah A, Lin H, Cresswell GM, et al. Pharmacologic Inhibition of FGFR Modulates the Metastatic Immune Microenvironment and Promotes Response to Immune Checkpoint Blockade. Cancer Immunol Res (2020) 8(12):1542–53. doi: 10.1158/2326-6066.CIR-20-0235

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Feng M, Grice DM, Faddy HM, Nguyen N, Leitch S, Wang Y, et al. Store-independent activation of Orai1 by SPCA2 in mammary tumors. Cell (2010) 143(1):84–98. doi: 10.1016/j.cell.2010.08.040

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Peretti M, Badaoui M, Girault A, Van Gulick L, Mabille MP, Tebbakha R, et al. Original association of ion transporters mediates the ECM-induced breast cancer cell survival: Kv10.1-Orai1-SPCA2 partnership. Sci Rep (2019) 9(1):1175. doi: 10.1038/s41598-018-37602-7

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Bi KW, Wei XG, Qin XX, Li B. BTK Has Potential to Be a Prognostic Factor for Lung Adenocarcinoma and an Indicator for Tumor Microenvironment Remodeling: A Study Based on TCGA Data Mining. Front Oncol (2020) 10:424. doi: 10.3389/fonc.2020.00424

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Li H, Zhao X, Wang J, Zong M, Yang H. Bioinformatics analysis of gene expression profile data to screen key genes involved in pulmonary sarcoidosis. Gene (2017) 596:98–104. doi: 10.1016/j.gene.2016.09.037

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Kuleshov MV, Jones MR, Rouillard AD, Fernandez NF, Duan Q, Wang Z, et al. Enrichr: a comprehensive gene set enrichment analysis web server 2016 update. Nucleic Acids Res (2016) 44(W1):W90–7. doi: 10.1093/nar/gkw377

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Jiang P, Gu S, Pan D, Fu J, Sahu A, Hu X, et al. Signatures of T cell dysfunction and exclusion predict cancer immunotherapy response. Nat Med (2018) 24(10):1550–8. doi: 10.1038/s41591-018-0136-1

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Arcangeli A. Ion channels and transporters in cancer. 3. Ion channels in the tumor cell-microenvironment cross talk. Am J Physiol Cell Physiol (2011) 301(4):C762–71. doi: 10.1152/ajpcell.00113.2011

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Karamitopoulou E. The Tumor Microenvironment of Pancreatic Cancer. Cancers (Basel) (2020) 12(10):3076. doi: 10.3390/cancers12103076

CrossRef Full Text | Google Scholar

14. Anderson NM, Simon MC. The tumor microenvironment. Curr Biol (2020) 30(16):R921–R5. doi: 10.1016/j.cub.2020.06.081

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Liu Y, Li X, Zhang Y, Wang H, Rong X, Peng J, et al. An miR-340-5p-macrophage feedback loop modulates the progression and tumor microenvironment of glioblastoma multiforme. Oncogene (2019) 38(49):7399–415. doi: 10.1038/s41388-019-0952-x

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Kumagai S, Togashi Y, Sakai C, Kawazoe A, Kawazu M, Ueno T, et al. An Oncogenic Alteration Creates a Microenvironment that Promotes Tumor Progression by Conferring a Metabolic Advantage to Regulatory T Cells. Immunity (2020) 53(1):187–203 e8. doi: 10.1016/j.immuni.2020.06.016

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Bustos SO, Antunes F, Rangel MC, Chammas R. Emerging Autophagy Functions Shape the Tumor Microenvironment and Play a Role in Cancer Progression - Implications for Cancer Therapy. Front Oncol (2020) 10:606436. doi: 10.3389/fonc.2020.606436

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Fabris L, Sato K, Alpini G, Strazzabosco M. The Tumor Microenvironment in Cholangiocarcinoma Progression. Hepatology (2020) 18(S1):75–85. doi: 10.1002/hep.31410

CrossRef Full Text | Google Scholar

19. Di Cosimo S. Advancing immunotherapy for early-stage triple-negative breast cancer. Lancet (2020) 396(10257):1046–8. doi: 10.1016/S0140-6736(20)31962-0

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Lieben L. Immunotherapy: Keeping breast cancer in check. Nat Rev Cancer (2017) 17(8):454–5. doi: 10.1038/nrc.2017.55

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Nathan MR, Schmid P. The emerging world of breast cancer immunotherapy. Breast (2018) 37:200–6. doi: 10.1016/j.breast.2017.05.013

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Salatino M, Girotti MR, Rabinovich GA. Glycans Pave the Way for Immunotherapy in Triple-Negative Breast Cancer. Cancer Cell (2018) 33(2):155–7. doi: 10.1016/j.ccell.2018.01.015

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Vanoevelen J, Dode L, Van Baelen K, Fairclough RJ, Missiaen L, Raeymaekers L, et al. The secretory pathway Ca2+/Mn2+-ATPase 2 is a Golgi-localized pump with high affinity for Ca2+ ions. J Biol Chem (2005) 280(24):22800–8. doi: 10.1074/jbc.M501026200

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Fenech MA, Carter MM, Stathopulos PB, Pin CL. The pancreas-specific form of secretory pathway calcium ATPase 2 regulates multiple pathways involved in calcium homeostasis. Biochim Biophys Acta Mol Cell Res (2020) 1867(1):118567. doi: 10.1016/j.bbamcr.2019.118567

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Newbury DF, Winchester L, Addis L, Paracchini S, Buckingham LL, Clark A, et al. CMIP and ATP2C2 modulate phonological short-term memory in language impairment. Am J Hum Genet (2009) 85(2):264–72. doi: 10.1016/j.ajhg.2009.07.004

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Jenkins J, Papkovsky DB, Dmitriev RI. The Ca2+/Mn2+-transporting SPCA2 pump is regulated by oxygen and cell density in colon cancer cells. Biochem J (2016) 473(16):2507–18. doi: 10.1042/BCJ20160477

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Zhao H, Zhang S, Shao S, Fang H. Identification of a Prognostic 3-Gene Risk Prediction Model for Thyroid Cancer. Front Endocrinol (Lausanne) (2020) 11:510. doi: 10.3389/fendo.2020.00510

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Tang K, Yu Y, Zhu L, Xu P, Chen J, Ma J, et al. Hypoxia-reprogrammed tricarboxylic acid cycle promotes the growth of human breast tumorigenic cells. Oncogene (2019) 38(44):6970–84. doi: 10.1038/s41388-019-0932-1

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Dang DK, Makena MR, Llongueras JP, Prasad H, Ko M, Bandral M, et al. A Ca(2+)-ATPase Regulates E-cadherin Biogenesis and Epithelial-Mesenchymal Transition in Breast Cancer Cells. Mol Cancer Res (2019) 17(8):1735–47. doi: 10.1158/1541-7786.MCR-19-0070

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Dang D, Prasad H, Rao R. Secretory pathway Ca(2+) -ATPases promote in vitro microcalcifications in breast cancer cells. Mol Carcinog (2017) 56(11):2474–85. doi: 10.1002/mc.22695

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Sui S, An X, Xu C, Li Z, Hua Y, Huang G, et al. An immune cell infiltration-based immune score model predicts prognosis and chemotherapy effects in breast cancer. Theranostics (2020) 10(26):11938–49. doi: 10.7150/thno.49451

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Grabovska Y, Mackay A, O’Hare P, Crosier S, Finetti M, Schwalbe EC, et al. Pediatric pan-central nervous system tumor analysis of immune-cell infiltration identifies correlates of antitumor immunity. Nat Commun (2020) 11(1):4324. doi: 10.1038/s41467-020-18070-y

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Braun DA, Hou Y, Bakouny Z, Ficial M, Sant’ Angelo M, Forman J, et al. Interplay of somatic alterations and immune infiltration modulates response to PD-1 blockade in advanced clear cell renal cell carcinoma. Nat Med (2020) 26(6):909–18. doi: 10.1038/s41591-020-0839-y

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Shapouri-Moghaddam A, Mohammadian S, Vazini H, Taghadosi M, Esmaeili SA, Mardani F, et al. Macrophage plasticity, polarization, and function in health and disease. J Cell Physiol (2018) 233(9):6425–40. doi: 10.1002/jcp.26429

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Hedbrant A, Wijkander J, Seidal T, Delbro D, Erlandsson A. Macrophages of M1 phenotype have properties that influence lung cancer cell progression. Tumour Biol (2015) 36(11):8715–25. doi: 10.1007/s13277-015-3630-9

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Valeta-Magara A, Gadi A, Volta V, Walters B, Arju R, Giashuddin S, et al. Inflammatory Breast Cancer Promotes Development of M2 Tumor-Associated Macrophages and Cancer Mesenchymal Cells through a Complex Chemokine Network. Cancer Res (2019) 79(13):3360–71. doi: 10.1158/0008-5472.CAN-17-2158

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Gu-Trantien C, Loi S, Garaud S, Equeter C, Libin M, de Wind A, et al. CD4(+) follicular helper T cell infiltration predicts breast cancer survival. J Clin Invest (2013) 123(7):2873–92. doi: 10.1172/JCI67428

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Gu-Trantien C, Migliori E, Buisseret L, de Wind A, Brohee S, Garaud S, et al. CXCL13-producing TFH cells link immune suppression and adaptive memory in human breast cancer. JCI Insight (2017) 2(11):e91487. doi: 10.1172/jci.insight.91487

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Davey MS, Willcox CR, Baker AT, Hunter S, Willcox BE. Recasting Human Vdelta1 Lymphocytes in an Adaptive Role. Trends Immunol (2018) 39(6):446–59. doi: 10.1016/j.it.2018.03.003

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Bank I, Book M, Huszar M, Baram Y, Schnirer I, Brenner H. V delta 2+ gamma delta T lymphocytes are cytotoxic to the MCF 7 breast carcinoma cell line and can be detected among the T cells that infiltrate breast tumors. Clin Immunol Immunopathol (1993) 67(1):17–24. doi: 10.1006/clin.1993.1040

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Morrow ES, Roseweir A, Edwards J. The role of gamma delta T lymphocytes in breast cancer: a review. Transl Res (2019) 203:88–96. doi: 10.1016/j.trsl.2018.08.005

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Aggarwal R, Lu J, Kanji S, Das M, Joseph M, Lustberg MB, et al. Human Vgamma2Vdelta2 T cells limit breast cancer growth by modulating cell survival-, apoptosis-related molecules and microenvironment in tumors. Int J Cancer (2013) 133(9):2133–44. doi: 10.1002/ijc.28217

PubMed Abstract | CrossRef Full Text | Google Scholar