1. Introduction
The stratification of breast cancer (BC) into morphologically and molecularly discrete subtypes has led to an overall improvement in treatment of patients with this disease [
1,
2]. TNBC, a subtype of BC with the poor prognosis is characterized by deficiency in expressions of
estrogen receptor (ER), progesterone receptor (PR), and human epidermal growth factor receptor 2 (
HER2). TNBC accounts for approximately 10-15% of all BCs and is mostly common in women younger than 40 years, women of color or those with a
BRCA1 mutation. TNBC differs from other types of invasive BC in that it tends to be more proliferative, highly metastatic, and has fewer therapeutic alternatives, therefore patients with this disease are inclined to have a worse prognosis [
2,
3,
4].
Enrichment of breast cancer stem cells (CSCs), a small population of tumor cells marked by high CD44 expression, low CD24 expression, and increased ALDH1A1/3 activity and the ability to self-renew and differentiate to repopulate the bulk tumor [
5], is believed to contribute to poor prognosis of TNBC [
6]. TNBC tumor heterogeneity driven by CSCs and lack of effective therapeutic arsenals other than chemotherapy often led to poor outcomes in patients with TNBC. In recent years, innovations in omics technologies have contributed immense knowledge about TNBC microenvironment heterogeneity [
7]. Moreover, identification of the key molecular determinants associated with the tumor heterogeneity and tumor progression in TNBC may help in development of novel strategies for treatment of this malignancy. In this study, we sought to find the critical genes associated with TNBC stemness through analysis of Gene Expression Omnibus (GEO) datasets, The Cancer Genome Atlas (TCGA), and Clinical Proteomics Tumor Analysis Consortium (CPTAC). We also employed the PGSEA algorithm, performed Gene Ontology (GO) analysis, and evaluated protein-protein interaction (PPI) networks to elucidate co-expressed genes with shared biological process terms (BP) and their linkage to stemness pathways in TNBC. Our analyses led to the identification of
CDK1, CCNA2,
CCNB1, AURKA and
EZH2 as signature genes involved in cancer stemness and progression. Moreover, we found that expressions of these Signature genes were affected by NAC1, a transcriptional co-factor belonging to the BTB gene family and possessing oncologic role [
8,
9]. The findings reported here may have potential and important implications in the prognosis and treatment of TNBC.
2. Materials and Methods
2.1. Identification and evaluation of differentially expressed genes
The following keywords: triple-negative breast cancer or TNBC, neoplasm, triple-negative breast cancer, RNA sequencing, TNBC patient samples, or breast cancer patient samples, were used to search in the Gene Expression Omnibus (GEO) database. We identified six TNBC datasets with adjacent normal tissues. The datasets were randomly grouped into training and validation groups. The training group comprised GSE65194, GSE36295, and GSE38959, while in the validation group, we had GSE21653, GSE20711, and GSE61725 datasets. GEO2R, an online platform using Limma R packages and GEOquery to deal with diverse experimental designs, was used to determine differentially expressed genes in tumors and normal samples. GEO2R analysis [
10] was performed with Limma-voom Benjamin and Hochberg and p≤0.05 as the settings for all datasets. All the other settings defaulted to GEO2R. We then characterized genes as differentially expressed between TNBC and normal if log2FC was ≥1.5 and the p-value was ≤ 0.05. Common genes in both the training and validation groups were obtained and used for downstream analysis.
2.2. GO analysis
GO and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways analysis was performed using Funrich software to determine enriched BP terms and oncogenic processes associated with the Differentially expressed genes [
11]. Integrated Differential Expression and Pathway analysis (iDEP), a platform that allows for selection of multiple pathways analysis methods and databases [
12,
13], was employed to perform Parametric Analysis of Gene Set Enrichment (PGSEA) [
14] to retrieve the genes involved in similar BP and differentially upregulated in TNBC.
2.3. Protein-protein interaction (PPI) network analysis and hub gene confirmation
To evaluate the possible PPI of the differentially expressed genes, Searching Tool for Recurring Instances of Neighboring Genes (STRING) [
15], Metascape [
16], and Cytoscape [
17] were used to construct and visualize the networks. UALCAN, a TCGA-based online platform for analysis of gene expression profiles in tumor and normal samples [
18,
19], was used to determine signature gene expressions in BC subtypes. TNM plot online platform [
20] was employed to compare the identified signature gene expressions in normal, primary, and metastatic samples. Assessment of mRNA expression and copy number alterations was executed using cBio cancer genomics portal (cbioportal) in TCGA and PanCancer Atlas samples [
21]. Immune cells infiltration in correlation with gene expression was analyzed using Timer 2.0, a comprehensive platform for systematical analysis of immune infiltrates in tumor (
https://cistrome.shinyapps.io/timer/) [
22]. Analysis of protein expression was performed using Human Protein Atlas database, a repository of tumor and normal IHC stained tissues [
23] (
https://www.proteinatlas.org), as well as the cBioPortal-Clinical Proteomic Tumor Analysis Consortium (CPTAC ) proteomics-pipeline (
cBioPortal for Cancer Genomics).
2.4. Analysis of the identified signature genes clinical relevance
The KM plotter [
24] was employed to evaluate the clinical significance of the identified signature genes. Using PAM50 a 50-gene signature classification system [
25], we analyzed the implication of the identified genes in TNBC (basal) samples.
2.5. In vitro expression of signature genes in TNBC cell lines
To determine protein expression in TNBC cells, we chose a wide range of cell lines across established TNBC subtypes; basal 1 (MDA-MB-468), basal 2 (HCC1806), mesenchymal (BT-459) and mesenchymal stem-like (MDA-MB-231). As controls, MCF10A and MCF7 were used as normal epithelial-like and non-TNBC cell lines, respectively. Cells were cultured following the ATCC recommendation. DMEM was used to culture MDA-MB-231, MCF7, and MDA-MB-468; DMEM-F12 was used for MDA-MB-453, and RPMI-1640 for HCC1806 and BT549 cell lines. MCF10A cells were cultured using mammary epithelial cell growth kit (promo-cell cat#C-21010). The media was supplemented with 10% fetal bovine serum (FBS) and 1% Penicillin-Streptomycin antibiotic mix. The cells were incubated at 37oC with 5% CO2 and 95% humidity.
2.6. Western Blot analysis
For western blot analysis, cells were harvested into a cold PBS buffer. Cells were spined down, and protein extraction was performed using Laemmli buffer supplemented with protease inhibitor. Protein quantification was done using a BCA assay kit (Thermo Scientific, #23228) and samples run in sodium dodecyl sulfate–polyacrylamide gel electrophoresis (SDS-PAGE) at 100 V. Protein blotting was done using immun-Blot® PVDF membrane in transfer buffer (100V for 75 minutes). The membranes were blocked with 5% defatted milk for 1 hour, and blots were incubated overnight with appropriate primary antibody at 4oC. Secondary antibody incubation was done at room temperature for 2 hours. TBST buffer was used to wash the membranes for 15 minutes each before and after adding horseradish peroxidase (HRP)-secondary antibody. Protein expression was visualized using the ChemiDOCTM MP imaging system (BioRAD).
3. Results
3.1. Identification of highly altered genes in TNBC
We retrieved six datasets from the GEO database and divided them randomly into two groups: training group and validation group. The training group comprised GSE65194, GSE36295, and GSE38959 datasets, which contained 96 TNBC, 29 normal, and 125 non-TNBC samples; the validation group consisted of GSE21653, GSE61725 and GSE20711 datasets, which had 140 TNBC,48 normal and 233 non-TNBC samples (
Table 1).
These datasets were subjected to a systematic analysis as depicted in
Figure 1. Differentially expressed genes in each dataset were analyzed using GEO2R tool (p≤0.05 and log2FC≥1.5). The common down- and up-regulated genes in each group are shown in
Figure 2a and S. Figure
1a-b. In the training group, we identified 112 common upregulated genes (
Figure 2a) and 74 downregulated genes (S. Figure
1c); in the validation group we observed 67 common upregulated and 23 downregulated genes (
Figure 2b and S. Figure
1d). Next, we compared the altered genes in the training group with those in the validation group and found that 64 (55 upregulated and 9 downregulated) genes were common between the training and validation groups (
Figure 2c, S. Figure
1e). A heatmap displaying the expression profiles of the common 64 genes in all 6 datasets is shown in
Figure 2d. In addition, utilizing the TCGA pan-cancer data, we explored the expression of these 64 genes among TNBC patient samples. Out of the 1084 TCGA BC samples, we analyzed 171 TNBC (basal-like) and 114 normal samples and found that these 64 genes were highly expressed in TNBC compared to normal samples (
Figure 2e).
3.2. Correlation of the highly expressed genes with cancer progression
To explore the implication of the highly expressed genes in the development and progression of TNBC, we performed a PPI analysis using STRING and Metascape [
16]. This analysis revealed possible interactions between the common 55 upregulated genes. The genes interaction network showed 55 nodes, 1359 edges, and an average node degree of 49.4. The average local clustering coefficient of the PPI network was 0.946, with a PPI enrichment
p-value of 1.0 x10
- 16 (
Figure 3a). The common 55 upregulated genes were further evaluated for their enrichment in the GO-biological process (GO-BP) terms. Intriguingly, we found that these genes were significantly enriched in the pathways associated with cell cycle, PLK1 and Aurora B kinase signaling, DNA replication, FOXM1 transcription factor (TF) network, G2/M checkpoints and p73 transcription factor network, which are highly dysregulated biological processes in many types of cancer (
Figure 3b). Additionally, we found that the regulator of TNBC progression, nuclear transcription factor Y subunit alpha (NFYA), was involved in the regulation of more than 50% of the 55 common upregulated genes (
Figure 3c). These results suggest that the upregulation of these genes may contribute to the aggressive phenotype of TNBC.
3.3. PGSEA analysis of the TNBC altered genes identifies a five gene-signature with a putative role in tumor stemness
We next performed a series of bioinformatics assessments and experimental validations to determine the effects of those altered genes on the stemness features of TNBC. Firstly, we performed the Parametric Gene Set Enrichment Analysis (PGSEA) to identify the genes with similar expression patterns and shared common GO-BP terms [
14]. This analysis pinpointed several relevant pathways, one of which is the regeneration term that includes
CDK1,
EZH2,
CCNB1,
CCNA2, and
AURKA (
Figure 3d). We then retrieved the background genes from the regeneration list (157 genes) and delved into the genes likely associated with cancer stemness/ cell self-renewal (S. Table
2). Among these genes, POU domain-class 5-transcription factor 3 (OCT4), SRY-box transcription factor 2 (SOX2), SRY-box transcription factor 9 (SOX9), Wnt signaling pathway proteins, catenin signaling genes, NOTCH signaling, matrix metallopeptidase 9 (MMP9), Krueppel-like factor 5 and MYC signaling genes, are all implicated in CSCs and self-renewal. In addition, the regeneration term contained several cytokines, such as interleukin-6 (IL-6) and interleukin 10 (IL-10). These analyses suggest that
CDK1,
EZH2,
CCNB1,
CCNA2, and
AURKA genes in regeneration term may have roles in driving the stemness of TNBC. Thus, we selected these five signature gene for further analysis.
Search tool for recurring instances of neighboring genes (STRING) and Cytoscape-cytoHubba validation analyses revealed a highly functional and physical association between these signature genes with a paired score degree values greater than 0.64 and a statistically significant enrichment (p<6.17x10-5). Interaction scores obtained from PPI analysis demonstrated that CDK1 has the highest connectivity scores with the other proteins (score >0.9) (S.Table 1). Additionally, CDK1 showed a high interaction with EZH2 that had a slightly lower interaction with the other three proteins (EZH2-AURKA, EZH2-CCNB1, and EZH2-CCNA2). To determine whether these 5 genes are significantly enriched in cancer-associated GO-biological processes, we performed a functional enrichment assessment using STRING. We found that these genes were not only enriched in cell regeneration, but also were significantly associated with other tumor-related biological processes, including histone phosphorylation, histone modification, and G2/M transition of mitotic cell cycle (FDR<0.0001) (S.Table 2). Furthermore, functional exploration of the 5 signature genes in the Wikipathways database affirms their roles in tumor-associated pathways, including retinoblastoma tumor regulation, DNA damage response, DNA repair, Ataxia telangiectasia Mutated (ATM) signaling pathway, miRNA regulation of DNA damage response, 5' adenosine monophosphate-activated protein kinase (AMPK) signaling, cell cycle regulation and DNA damage response (FDR<0.02) (S.Table 3).
To determine whether the functions of the signature genes were unique to the TNBC phenotype, we compared their expression in TNBC with other BC subtypes at both transcript and protein levels using the TCGA dataset samples. Our analysis found that except CCNB1, whose increase in TNBC did not attain a significant level compared to HER2 positive samples, all the other genes were highly expressed in TNBC samples as compared to the normal or other BC subtypes (p<0.05) (
Figure 4a). Also, using the CPTAC Platform samples, we analyzed the protein expression of the signature genes in TNBC samples (
Figure 4b). A significant upregulation of AURORA-A, CCNB1, CDK1, and EZH2 was observed in TNBC as compared to other types of BCs (O-BCs) (p<0.05). CCNA2 expression in TNBC was not significantly higher than that in O-BCs. Analysis of the signature genes protein expression using the Human Protein Atlas also showed high expressions of these proteins in tumor samples compared to normal samples (S.
Figure 2). Moreover, our western blot analysis verified the higher expressions of EZH2 and AURORA-A in MDA-MB-231 and MDA-MB-468 TNBC cells than in non-malignant breast cells (
Figure 4c and S. Figure
3a). Highest expression of CCNB1 was observed in mesenchymal TNBC cell line BT549 (
Figure 4c). HCC1806, a lesser aggressive TNBC cell line, had a lower expression of the signature genes than other TNBC cell lines (
Figure 4c).
To explore the possible roles of the identified signature genes in tumor metastasis, we analyzed the TCGA data utilizing the tumor-normal-metastasis (TNM) plotter platform. The Pearson correlation analysis demonstrated a significantly increased expression of the signature genes in primary and metastatic samples compared to normal samples: CDK1, p<4.98 x10
-56; CCNA2, p<4.42x10
-16; CCNB1, p<1.5x10
-33; AURKA, p<1.5x10
-18; and EZH2, p<6.94x10
-7 (
Figure 4d). This analysis not only confirms the increased expression of these genes in TNBC but also supports our hypothesis that these genes have important roles in TNBC tumor progression. Also, the Receiver Operating Characteristic (ROC) analysis displayed the specificity of these genes to TNBC versus non-TNBC, as illustrated by the high area under the curve (AUC) values: CCNA2: 0.857 (p<0.00001); EZH2: 0.844 (p<0.00001); AURKA: 0.784 (P<0.00001); CCNB1: 0.762 (P<0.00001), and CDK1: 0.722 (P<0.00001) (
Figure 4e). Furthermore, using the KM-Plotter platform, we assessed the clinical implication of the identified signature genes in prognosis of TNBC patients. Among the 5 signature genes, the expression of EZH2 and CCNB1 appeared to highly affect the clinical prognosis of TNBC patients (
Figure 4f and S. Figure
3a)
. Taken together, these analyses suggest that the high expression of the identified signature genes in TNBC may contribute substantially to the aggressive phenotype of TNBC including tumor stemness genes.
3.4. Signature genes expression is positively associated with immunosuppressive cells infiltration and hypoxia status in TNBC
Recent studies have revealed an association between stemness and immune signatures [
27,
28,
29]. It is believed that CSCs in highly hypoxic regions of the tumor microenvironment (TME) can interact with immune cells, including tumor-associated macrophages, myeloid-derived suppressor cells (MDSCs), and T cells [
30]. Reciprocally, the immune cells such as MDSCs and Tregs can nurture stemness niche to support CSC survival [
29]. We, therefore, ascertain whether the expressions of the signature genes have an association with the infiltration of Tregs and MDSCs in TNBC samples. TNBC-TCGA samples from the Tumor Immune Estimation Resource (TIMER)2.0. were utilized, and the published guideline for choosing transcriptome quantification methods was followed to determine the influence of gene expression on immunity [
31]. We selected the Tumor Immune Dysfunction and Exclusion (TIDE) platform [
32] for analysis of MDSCs infiltration and computational pipeline for quantification of the Tumor Immune contexture from human RNA-seq data (QUANTISEQ) for Tregs. Notably, there was a significant positive correlation between the mRNA expressions of AURKA, CCNA2, EZH2, and CDK1 with infiltration of MDSCs (P<0.05). Additionally, there was a significant positive correlation between enhanced AURKA, CCNA2, EZH2 and CDK1 mRNA expression with MDSCs and Tregs infiltration (P<0.05) (
Figure 5a). In the analysis of cell types in TME, we observed that the stromal and immune scores are positively correlated with the samples of low tumor purity [
33] and that tumor purity is associated with the expression of signature genes, suggesting that the source of signature genes mRNA transcriptome was from tumor cells but not from infiltrated immune cells (
Figure 5a). Additionally, explorative analysis of the role of these genes in predicting patients’ response to immune-checkpoint inhibitors (ICIs) including anti-PD1, anti-PD-L1, or anti-CTLA4 antibodies, we found that EZH2 (P<0.05) and AURKA (P<0.05) could better predict outcomes of patient receiving these treatments (
Figure 5b-d). Retrieving the immune-TME and TNBC aggression-associated marker genes from the Pan-cancer TCGA-dataset, we found a highly positive correlation of the expression of S100A9 and S100A8, two immune suppression marker genes, with AURKA and EZH2 in the TNBC samples. Similarly, correlation existed between Arginase 1 (ARG1), which is known to dampen T-cell proliferation, with EZH2 and AURKA. Additionally, high CCNB1 expression was correlated with elevated expression of KDM5A (a drug resistance marker) while AURKA positively correlated with HAVCR2 (TIM3) (a marker for highly exhausted T-cells) (
Figure 5e). These analyses imply that the signature genes identified here may have potential roles in modulating immune-TME and the interaction between immune cells and CSCs, contributing to the immunosuppressive phenotype of TNBC.
One of the features of TME in solid tumors is oxygen deficiency, which favors CSCs enrichment and MDSCs accumulation, contributing to an immune-suppressive environment [
34,
35]. Thus, we queried whether these signature genes are involved in the modulation of hypoxia. We retrieved the transcriptomic and copy number alterations (CNAs) data from TCGA and quantified tumor hypoxia in TNBC and non-TNBC samples, using three independent hypoxia gene signature-based scoring tools: Winter [
36], Ragnum [
37], and Buffa [
38]. The samples were stratified into TNBC (basal-like) and non-TNBC and evaluated for the differences in hypoxic status. We found that the hypoxic level was significantly higher in TNBC samples than that in non-TNBC (p<0.0001) (
Figure 6a). Sorting the samples based on Ragnum score hypoxia intensity and examining the CNAs and mRNA expression of the signature genes uncovered a correlation between mRNA amplification and CNAs with higher hypoxia intensity (S. Figure
4a).
To determine whether the expression patterns of CNAs and mRNA were also exhibited at the protein level, we explored the CPTAC pan-cancer data. A positive correlation between the expression of the 5 proteins and hypoxia scores status was observed (S. Figure
4b), and all three methods demonstrated similar changes and correlations (correlation coefficient >0.2) (S. Figure
4c). The proteins that had the highest correlations with the hypoxia scores were CDK1 (Buffa score: 0.60; Winter score: 0.51; Ragnum score: 0.52) and EZH2 (Buffa score: 0.52; Winter score: 0.45; Ragnum score: 0.44). Furthermore, EZH2-CDK1 pair demonstrated a high correlation (correlation coefficient; 0.62), while AURORA-A expression was highly correlated with CDK1 protein levels (correlation coefficient; 0.53) under hypoxia conditions (S. Figure
4c). To determine whether these 5 signature genes can predict hypoxia status in BC, we stratified the TCGA samples into two groups based on the Buffa score. The samples with positive Buffa scores were grouped as positive responders, while those with zero or negative scores were grouped as non-responders. ROC analysis demonstrated high specificity of these genes in predicting the samples with high hypoxia (
Figure 6b). Intriguingly, AUC values of the signature genes matched those of Ragnum (0.88, p<0.00001) and Winter (0.951, p<0.00001) hypoxia prediction scores: AURKA, 0.847 (P<0.00001); CCNB1, 0.827 (P<0.00001); CDK1, 0.806 (P<0.00001); CCNA2; 0.835 (P<0.00001); and EZH2; 0.812 (P<0.00001). Together, these analyses suggest that the signature genes identified in this study may have important roles in regulating or predicting hypoxic status in TNBC.
3.5. Depletion of NACC1 reduces the signature genes expression
We and others previously demonstrated the oncogenic role of
NACC1 in various types of cancer and in supporting cell survival under hypoxia conditions [
39]. More recently, we showed that NAC1 expression in melanoma cells promotes their evasion from immune surveillance [
40]. Thus, we asked whether this transcription co-factor has any effects on those signature genes in TNBC cells. Analysis of GSE183947 GEO dataset [
41] of the matched normal, primary, and metastatic TNBC samples found a significant upregulation of
NACC1 mRNA level in primary and metastatic tumors compared to normal tissues (
Figure 7a). It has been shown that p27 is enriched in less proliferative cells [
5]. We identified DEGs between p27-low and p27-high TNBC samples in GSE198713 dataset [
5]. Surprisingly, we found that the signature genes were significantly downregulated in p27 samples (
Figure 7b). Consistently, our experiments demonstrated that depletion of NAC1 in MDA-MB-231 TNBC cells marked reduced AURORA-A, CCNB1, and EZH2 expressions (
Figure 7c). Oncogenic driver gene signature analysis using GSEA on NACC1 depleted cells RNA-seq data showed that the downregulated genes were enriched in PRC2-EZH2, KRAS, BC-associated genes, epidermal growth factor receptor (EGFR), AKT signaling and WNT pathway (S. Figure
5a). Furthermore, our western analysis on AKT pathway-associated gene expression in the NACC1-depleted tumor cells showed decreased expressions of phospho-PTEN, β-catenin and cyclin-D1 (S. Figure
5b), which is consistent with the RNA-seq data. Collectively, these data suggest that high expression of NAC1 in TNBC may promote tumor progression and hypoxia adaption through affecting the oncogenic and stemness-associated pathways and genes, including the 5 signature genes identified in this study.
4. Discussion
Intra and inter-tumor heterogeneity, tumor relapse, metastasis, and therapy resistance constitute the major challenges in treatment of TNBC, and CSCs are an important contributor to these malignant phenotypes [
42]. To better understand these problems, we explored the genes that drive the malignant phenotype of TNBC. We mined six GEO datasets and TCGA samples and analyzed the altered gene expressions at both transcriptome and protein levels using PGSEA algorithm, STRING, Metascape, and TIMER2.0. Through these analyses, we identified 55 upregulated and 9 downregulated genes implicated in the regulation of tumor aggressive phenotype in TNBC. Further exploration led to the identification and validation of
CDK1,
CCNB1,
AURKA,
EZH2, and
CCNA2, a five-gene signature with high interactions and potentially associated with tumor stemness, hypoxia enhancement and TME regulation. Consistently, these 5 genes are among those reported to be potential prognostic biomarker signature for hepatocellular carcinoma [
43].
Protein-protein interactions, as well as gene co-expression, play critical roles in the regulation of various malignant features. Here, we found a close association of the 55 common differentially expressed genes as well as the identified signature genes at mRNA and protein levels, especially for CDK1-EZH2 pairing. Strong correlation of these genes at mRNA level may suggest a shared transcription regulation machinery. We further identified NFYA, known to promote tumor progression in TNBC [
44], as the most enriched transcription factor targeting over 50% of the upregulated genes. The short-isoform of NFYA transcriptionally amplifies the survival and tumor cell proliferation-associated genes. On the other hand, the long NFYA isoform enhances cell signaling through EMT transition and leads to a more belligerent phenotype of Claudin
low basal-like tumors [
44]. Similarly, the expression of NYFA and sex determinant region Y are independent prognostic markers in gastric cancer, and that NFYA promotes clear cell renal cell carcinoma tumor growth through regulation of cyclins activity [
45,
46]. These observations support the hypothesis that upregulation of the signature genes transcription and their protein translation in TNBC play an important role in tumor progression.
Recent studies have demonstrated that interaction and phosphorylation of SOX2 by CDK1 positively regulate stemness in lung cancer [
47], and CDK1 phosphorylation of both HIF-1α and EZH2 promotes tumor cells survival and proliferation [
48,
49]. Our analyses show that CDK1 is highly correlated with EZH2 as well as other signature genes in TNBC. Additionally, we found that both CDK1 and EZH2 are strongly correlated with hypoxia status at transcription, CNAs, and protein translation levels. Under hypoxia, tumor cells can stimulate MDSCs to establish a premetastatic niche that promotes cancer cells aggression [
35]. Although tumor cells have the capability to adapt to hypoxia, the hypoxia-related products can deter activation of tumor-infiltrating immune cells, establishing an immunosuppressive environment critical for tumor cell evasion from immune surveillance [
35]. IL-6 is a cytokine involved in immunosuppression and its expression is upregulated by CDK1 activity [
50].Under hypoxia, IL-6 can upregulate HIF1α to enhance CSCs survival under cisplatin therapy [
51]. We show here that the upregulation of signature genes (especially CDK1, EZH2, and AURKA) is correlated with hypoxia in TNBC at both mRNA and protein levels, and that these signature genes are highly expressed in primary and metastatic samples. Increased expression of EZH2 can enhance TNBC cells invasion and metastasis through suppression of
TIMP2 transcription and subsequent amplification of MMP-2 and MMP-9 activity [
52]. The role of hypoxia in regulation of immune-TME of TNBC is emerging, and studies have demonstrated the importance of hypoxia in immune-tumor cells interaction that drive tumor progression [
53]. The high correlation of the signature genes with tumor hypoxia, stemness status and increased immune cells infiltration may provide new insights into this issue. Our ROC analysis further collaborates these observations and implies the potential of the identified signature genes in predicting the hypoxia status of TNBC.
More recently, we reported that NAC1 expression in melanoma cells was critical for immune evasion [
40] and that NACC1 expression negatively regulates the suppressive activity of Tregs [
54] and CD8
+ memory T-cells [
55]. The current study suggests that NAC1, likely through regulation of the identified signature genes, has a role in immunosuppressive TME. It was reported that the expression of these genes in HCC positively correlates with infiltrations of neutrophils, macrophages, and dendritic cells and could predict HCC patients response to anti-PD-1, anti-PD-L1, and anti-CTLA-4 antibodies [
56]. We found here that high expression of signature genes is significantly associated with increased infiltration of MDSCs and Tregs. The infiltration of both Tregs and MDSCs suggests multiple roles of the signature genes in modulating immuno-TME in TNBC. Liu
et al., [
57] demonstrated that CCNB1 and EZH2 together could influence immune cells infiltration and serve as a prognostic hub in prostate cancer. However, our findings indicate that expression of EZH2 but not CCNB1 is associated with poor prognosis in patients treated with ICIs, and that a combination of EZH2 and NACC1 may better predict the survival of patients receiving anti-PDL1 treatment.
In conclusion, we have identified five NAC1-regulated genes implicated in TNBC aggressive phenotype, including tumor stemness and immune cells infiltration. The expressions of these genes could be further explored as useful prognostic markers of TNBC. The protein products of these genes may also be exploited as potential targets for the treatment of TNBC.
Author Contributions
Conceptualization, C.M.N., J.M. Y. and X.L.; formal analysis: methodology, C.M.N, H.H., A.O.A, O.O.; validation, C.M.N; investigation, C.M.N., O.O., A.O.A., F.F.O., H.H., F.O., X. X., X.R. and J.S.; resources, J.M. Y., J.S. and X.L.; writing—original draft preparation, C.M.N., J.M. Y., X.L., writing—review and editing, H.H., X.L., J.M. Y., O.O, and A.O.A; funding acquisition, J.M. Y., J.S. and X.L. All authors have read and agreed to the published version of the manuscript.
Funding
This work was supported by National Cancer Institute Grant R01CA221867 to J.Y. and J.S., Susan G. Komen Foundation CCR18548501 to X.L., and NIH grant P20GM121327 to X.L.
Institutional review board statement
Not applicable.
Informed consent statement
Informed consent was obtained from all subjects involved in the study.
Data availability
The dataset generated for this study are available on request to the corresponding author.
Declarations
The authors have no conflict of interests to declare.
References
- Masood, S. Breast Cancer Subtypes: Morphologic and Biologic Characterization. Womens Health 2016, 12, 103–119. [Google Scholar] [CrossRef] [PubMed]
- Mahmoud, R.; Ordonez-Moran, P.; Allegrucci, C. Challenges for Triple Negative Breast Cancer Treatment: Defeating Heterogeneity and Cancer Stemness. Cancers 2022, 14. [Google Scholar] [CrossRef] [PubMed]
- Chen, J.Q.; Russo, J. ERalpha-negative and triple negative breast cancer: Molecular features and potential therapeutic approaches. Biochim Biophys Acta 2009, 1796, 162–175. [Google Scholar] [PubMed]
- Fultang, N.; Chakraborty, M.; Peethambaran, B. Regulation of cancer stem cells in triple negative breast cancer. Cancer Drug Resist 2021, 4, 321–342. [Google Scholar] [CrossRef]
- Baldominos, P.; et al. Quiescent cancer cells resist T cell attack by forming an immunosuppressive niche. Cell 2022, 185, 1694–1708 e19. [Google Scholar] [CrossRef] [PubMed]
- He, L.; et al. The Role of Breast Cancer Stem Cells in Chemoresistance and Metastasis in Triple-Negative Breast Cancer. Cancers (Basel) 2021, 13. [Google Scholar] [CrossRef]
- Bianchini, G.; et al. Treatment landscape of triple-negative breast cancer—Expanded options, evolving needs. Nat Rev Clin Oncol 2022, 19, 91–113. [Google Scholar] [CrossRef]
- Nakayama, K.; et al. A BTB/POZ protein, NAC-1, is related to tumor recurrence and is essential for tumor growth and survival. Proceedings of the National Academy of Sciences of the United States of America 2006, 103, 18739–18744. [Google Scholar] [CrossRef]
- Perez-Torrado, R.; Yamada, D.; Defossez, P.A. Born to bind: The BTB protein-protein interaction domain. Bioessays 2006, 28, 1194–1202. [Google Scholar] [CrossRef]
- Sean, D.; Meltzer, P.S. GEOquery: A bridge between the gene expression omnibus (GEO) and BioConductor. Bioinformatics 2007, 23, 1846–1847. [Google Scholar]
- Pathan, M.; et al. FunRich: An open access standalone functional enrichment and interaction network analysis tool. Proteomics 2015, 15, 2597–2601. [Google Scholar] [CrossRef] [PubMed]
- Ge, X.J. iDEP Web Application for RNA-Seq Data Analysis. In Rna Bioinformatics, 2nd edition; 2021; Volume 2284, pp. 417–443. [Google Scholar]
- Ge, S.X.; Son, E.W.; Yao, R.N. iDEP: An integrated web application for differential expression and pathway analysis of RNA-Seq data. Bmc Bioinformatics 2018, 19. [Google Scholar] [CrossRef] [PubMed]
- Kim, S.Y.; Volsky, D.J. PAGE: Parametric analysis of gene set enrichment. Bmc Bioinformatics 2005, 6. [Google Scholar] [CrossRef] [PubMed]
- Szklarczyk, D.; et al. Correction to 'The STRING database in 2021: Customizable protein-protein networks, and functional characterization of user-uploaded gene/measurement sets'. Nucleic Acids Res 2021, 49, 10800. [Google Scholar] [CrossRef]
- Zhou, Y.; et al. Metascape provides a biologist-oriented resource for the analysis of systems-level datasets. Nat Commun 2019, 10, 1523. [Google Scholar] [CrossRef]
- Szklarczyk, D.; et al. The STRING database in 2021: Customizable protein-protein networks, and functional characterization of user-uploaded gene/measurement sets (vol 49, pg D605, 2021). Nucleic Acids Research 2021, 49, 10800. [Google Scholar] [CrossRef]
- Chandrashekar, D.S.; et al. UALCAN: An update to the integrated cancer data analysis platform. Neoplasia 2022, 25, 18–27. [Google Scholar] [CrossRef]
- Chandrashekar, D.S.; et al. UALCAN: A Portal for Facilitating Tumor Subgroup Gene Expression and Survival Analyses. Neoplasia 2017, 19, 649–658. [Google Scholar] [CrossRef]
- Bartha, A.; Gyorffy, B. TNMplot.com: A Web Tool for the Comparison of Gene Expression in Normal, Tumor and Metastatic Tissues. International Journal of Molecular Sciences 2021, 22. [Google Scholar] [CrossRef]
- Gao, J.J.; et al. Integrative Analysis of Complex Cancer Genomics and Clinical Profiles Using the cBioPortal. Science Signaling 2013, 6. [Google Scholar] [CrossRef]
- Li, T.W.; et al. TIMER2.0 for analysis of tumor-infiltrating immune cells. Nucleic Acids Research 2020, 48, W509–W514. [Google Scholar] [CrossRef] [PubMed]
- Uhlen, M.; et al. A pathology atlas of the human cancer transcriptome. Science 2017, 357, 660-+. [Google Scholar] [CrossRef] [PubMed]
- Lanczky, A.; Gyorffy, B. Web-Based Survival Analysis Tool Tailored for Medical Research (KMplot): Development and Implementation. Journal of Medical Internet Research 2021, 23. [Google Scholar] [CrossRef] [PubMed]
- Mathews, J.C.; et al. Robust and interpretable PAM50 reclassification exhibits survival advantage for myoepithelial and immune phenotypes. NPJ Breast Cancer 2019, 5, 30. [Google Scholar] [CrossRef]
- Wan, L.; et al. Phosphorylation of EZH2 by AMPK Suppresses PRC2 Methyltransferase Activity and Oncogenic Function. Mol Cell 2018, 69, 279–291 e5. [Google Scholar] [CrossRef]
- Chen, P.; et al. Cancer Stemness Meets Immunity: From Mechanism to Therapy. Cell Rep 2021, 34, 108597. [Google Scholar] [CrossRef]
- Yi, L.; et al. Integrative stemness characteristics associated with prognosis and the immune microenvironment in esophageal cancer. Pharmacol Res 2020, 161, 105144. [Google Scholar] [CrossRef]
- Bayik, D.; Lathia, J.D. Cancer stem cell-immune cell crosstalk in tumour progression. Nat Rev Cancer 2021, 21, 526–536. [Google Scholar] [CrossRef]
- Miranda, A.; et al. Cancer stemness, intratumoral heterogeneity, and immune response across cancers. Proc Natl Acad Sci U S A 2019, 116, 9020–9029. [Google Scholar] [CrossRef]
- Sturm, G.; et al. Comprehensive evaluation of transcriptome-based cell-type quantification methods for immuno-oncology. Bioinformatics 2019, 35, I436–I445. [Google Scholar] [CrossRef]
- Jiang, P.; et al. Signatures of T-cell dysfunction and exclusion predict cancer immunotherapy response. Cancer Immunology Research 2019, 7. [Google Scholar] [CrossRef]
- Yoshihara, K.; et al. Inferring tumour purity and stromal and immune cell admixture from expression data. Nature Communications 2013, 4. [Google Scholar] [CrossRef] [PubMed]
- Bhandari, V.; et al. Molecular landmarks of tumor hypoxia across cancer types. Nat Genet 2019, 51, 308–318. [Google Scholar] [CrossRef] [PubMed]
- Hu, J.Y.; et al. Hypoxia, a key factor in the immune microenvironment. Biomedicine & Pharmacotherapy 2022, 151. [Google Scholar]
- Winter, S.C.; et al. Relation of a hypoxia metagene derived from head and neck cancer to prognosis of multiple cancers. Cancer Res 2007, 67, 3441–3449. [Google Scholar] [CrossRef]
- Ragnum, H.B.; et al. The tumour hypoxia marker pimonidazole reflects a transcriptional programme associated with aggressive prostate cancer. Br J Cancer 2015, 112, 382–390. [Google Scholar] [CrossRef]
- Buffa, F.M.; et al. Large meta-analysis of multiple cancers reveals a common, compact and highly prognostic hypoxia metagene. Br J Cancer 2010, 102, 428–435. [Google Scholar] [CrossRef]
- Ren, Y.J.; et al. Silencing of NAC1 Expression Induces Cancer Cells Oxidative Stress in Hypoxia and Potentiates the Therapeutic Activity of Elesclomol. Frontiers in Pharmacology 2017, 8. [Google Scholar] [CrossRef]
- Ren, Y.J.; et al. Tumorous expression of NAC1 restrains antitumor immunity through the LDHA-mediated immune evasion. Journal for Immunotherapy of Cancer 2022, 10. [Google Scholar] [CrossRef]
- Zhang, Y.; et al. Identification of Five Cytotoxicity-Related Genes Involved in the Progression of Triple-Negative Breast Cancer. Frontiers in Genetics 2022, 12. [Google Scholar] [CrossRef]
- Das, P.K.; et al. Plasticity of Cancer Stem Cell: Origin and Role in Disease Progression and Therapy Resistance. Stem Cell Rev Rep 2020, 16, 397–412. [Google Scholar] [CrossRef] [PubMed]
- Ameri, M.; Salimi, H.; Eskandari, S.; Nezafat, N. Identification of potential biomarkers in hepatocellular carcinoma: A network-based approach. Informatics in Medicine Unlocked 2022, 28, 100864. [Google Scholar] [CrossRef]
- Dolfini, D.; Andrioletti, V.; Mantovani, R. Overexpression and alternative splicing of NF-YA in breast cancer. Scientific Reports 2019, 9. [Google Scholar] [CrossRef]
- Cao, B.; et al. Gene regulatory network construction identified NFYA as a diffuse subtype-specific prognostic factor in gastric cancer. International Journal of Oncology 2018, 53, 1857–1868. [Google Scholar] [CrossRef]
- Li, Y.; et al. Transcription factor NFYA promotes G1/S cell cycle transition and cell proliferation by transactivating cyclin D1 and CDK4 in clear cell renal cell carcinoma. Am J Cancer Res 2020, 10, 2446–2463. [Google Scholar] [PubMed]
- Huang, Z.; Shen, G.; Gao, J. CDK1 promotes the stemness of lung cancer cells through interacting with Sox2. Clinical & Translational Oncology 2021, 23, 1743–1751. [Google Scholar]
- Chen, G.F.; et al. Hypoxia induces an endometrial cancer stem-like cell phenotype via HIF-dependent demethylation of SOX2 mRNA. Oncogenesis 2020, 9. [Google Scholar] [CrossRef]
- Tang, B.; et al. ZMYND8 preferentially binds phosphorylated EZH2 to promote a PRC2-dependent to -independent function switch in hypoxia-inducible factor-activated cancer. Proceedings of the National Academy of Sciences of the United States of America 2021, 118. [Google Scholar] [CrossRef]
- Kuang, Y.; et al. Iron-dependent CDK1 activity promotes lung carcinogenesis via activation of the GP130/STAT3 signaling pathway. Cell Death Dis 2019, 10, 297. [Google Scholar] [CrossRef]
- Zhang, F.Q.; et al. Cisplatin treatment increases stemness through upregulation of hypoxia-inducible factors by interleukin-6 in non-small cell lung cancer. Cancer Science 2016, 107, 746–754. [Google Scholar] [CrossRef] [PubMed]
- Chien, Y.C.; et al. EZH2 promotes migration and invasion of triple-negative breast cancer cells via regulating TIMP2-MMP-2/-9 pathway. Am J Cancer Res 2018, 8, 422–434. [Google Scholar]
- Abou Khouzam, R.; et al. Tumor Hypoxia Regulates Immune Escape/Invasion: Influence on Angiogenesis and Potential Impact of Hypoxic Biomarkers on Cancer Therapies. Frontiers in Immunology 2021, 11. [Google Scholar] [CrossRef] [PubMed]
- Yang, J.M.; et al. NAC1 modulates autoimmunity by suppressing regulatory T cell-mediated tolerance. Science Advances 2022, 8. [Google Scholar] [CrossRef]
- Wang, L.Q.; et al. Expression of NAC1 Restrains the Memory Formation of CD8(+) T Cells during Viral Infection. Viruses-Basel 2022, 14. [Google Scholar] [CrossRef] [PubMed]
- Zou, Y.P.; et al. CDK1, CCNB1, and CCNB2 are Prognostic Biomarkers and Correlated with Immune Infiltration in Hepatocellular Carcinoma. Medical Science Monitor 2020, 26. [Google Scholar] [CrossRef]
- Liu, R.J.; et al. Identification of biomarkers, pathways and potential therapeutic target for docetaxel resistant prostate cancer. Bioengineered 2021, 12, 2377–2388. [Google Scholar] [CrossRef] [PubMed]
|
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content. |
© 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).