Identification of a Gene Panel Predictive of Triple-Negative Breast Cancer Response to Neoadjuvant Chemotherapy Employing Transcriptomic and Functional Validation

Triple-negative breast cancer (TNBC) patients exhibiting pathological complete response (pCR) have better clinical outcomes compared to those with residual disease (RD). Therefore, robust biomarkers that can predict pCR may help with triage and resource prioritization in patients with TNBC. Herein, we identified a gene panel predictive of RD and pCR in TNBC from the discovery (n = 90) treatment-naive tumor transcriptomic data. Eight RD-derived genes were identified as TNBC-essential genes, which were highly predicative of overall survival (OS) and relapse-free survival (RFS) in an additional cohort of basal breast cancer (n = 442). Mechanistically, targeted depletion of the eight genes reduced the proliferation potential of TNBC cell models, while most remarkable effects were for combined SLC39A7, TIMM13, BANF1, and MVD knockdown in conjunction with doxorubicin. Orthogonal partial least squares-discriminant analysis (OPLS-DA) and receiver operating characteristic curve (ROC) analyses revealed significant predictive power for the identified gene panels with an area under the curve (AUC) of 0.75 for the validation cohort (n = 50) to discriminate RD from pCR. Protein–Protein Interaction (PPI) network analysis of the pCR-derived gene signature identified an 87-immune gene signature highly predictive of pCR, which correlated with better OS, RFS, and distant-metastasis-free survival (DMFS) in an independent cohort of basal and, to a lesser extent, HER2+ breast cancer. Our data have identified gene signatures predicative of RD and pCR in TNBC with potential clinical implications.


Introduction
Triple-negative breast cancer (TNBC) represents 15-20% of all breast cancers, which is characterized by the lack of hormone receptors (HR), including estrogen receptor (ER) and progesterone receptor (PR) expression, in addition to the lack of epidermal growth factor receptor-2 (HER-2) amplification [1]. Despite the lower incidence of TNBC compared to other molecular subtypes, TNBC is oftentimes presented with an advanced tumor stage and those patients have worse overall survival (OS) [2]. Neoadjuvant chemotherapy (NAC) followed by surgery is the standard-of-care clinical practice for TNBC with larger tumors [3]. Longitudinal patterns of TNBC response to NAC are oftentimes used to assess treatment response and likelihood for OS. Pathological complete response (pCR) is associated with better prognosis, while residual disease (RD) predicts worse clinical outcome [4]. Therefore, the timely identification of patients who will not benefit from NAC is critical to avoid exposure of this patient group to toxic drugs and unnecessary delays in surgical intervention.
While tumor heterogeneity and clonal evolution have been linked to variable responses to NAC in TNBC employing single-cell analysis [5,6], currently there are no predictive signatures for RD and pCR in the context of NAC in clinical practice. We previously reported the presence of immune gene expression as the hallmark of pCR in TNBC employing single-cell analysis [7]. In the current study, we utilized publicly available datasets and characterized the transcriptome from a discovery cohort (n = 90) transcriptomic data and subsequently constructed an orthogonal partial least squares-discriminant analysis (OPLS-DA) and validated it on a second cohort (n = 50) of TNBC and identified the gene signatures predicative of RD and pCR. Interestingly, we observed genes associated with pCR to have better prediction (higher AUC) compared to genes associated with RD, which implied inherent heterogeneity in the genes driving RD, while the genes predictive of pCR were mostly indicative of immune infiltration. Employing CRISPR-Cas9 in vitro screen data, we identified eight RD-derived genes essential for TNBC, which upon targeted depletion reduced TNBC colony formation and enhanced their sensitivity to doxorubicin and paclitaxel.

Differential Expression Analysis of RD Compared to pCR TNBC Tumors
To identify putative markers predictive of RD and pCR in the context of NAC therapy for TNBC, we retrieved RNA-Seq data from 90 TNBC patients who exhibited variable responses to NAC (pCR; n = 38 and RD; n = 52). The data were mapped to the gencode release v33 followed by comparative analysis, which revealed several differentially expressed genes, with most of the genes exhibiting downregulated expression in RD compared to pCR (1.5 FC, p < 0.05, Figure 1a and Table S1). Differentially expressed genes in RD and pCR are illustrated as a volcano plot (Figure 1b), with upregulated genes in RD are shown in red and downregulated genes in blue. To gain deeper insight into functional categories associated with RD vs. pCR TNBC, the list of differentially expressed genes were subjected to ingenuity pathway analysis (IPA). Canonical pathway analysis on the differentially expressed genes revealed predominant under representation of categories associated with immune response (Figure 1c). Figure 1d provides a high-level tree map of altered downstream functional categories utilizing the list of differentially expressed genes in RD vs. pCR. The major-colored rectangles indicate a group of associated biological functions or diseases, where blue indicate decreasing and orange indicates increasing, and the dimension of the rectangles indicates where associated functions are predicted to be up or down, while the color intensity corresponds to the absolute Z-scores. Using this analysis, IPA revealed suppression of various immune functions as the hallmark associated with RD, implying activation of immune functions as the most significant surrogate marker to pCR in TNBC in the context of NAC response (Figure 1d and Table S2). . TNBC patients were grouped according to pCR (n = 38) and residual disease (RD, n = 62) and were subjected to differential analysis (1.5 FC, p< 0.05). (a) Hierarchical clustering of TNBC as a function of pCR (n = 38) or RD (n = 62) based on differentially expressed mRNA genes. Each column represents one patient sample, and each row represents an mRNA. The expression level of each gene (log2) in a single sample is shown according to the color scale. (b) Volcano plot illustrating the upregulated (red) or downregulated (blue) genes in RD vs. pCR TNBC. The most differentially expressed mRNAs are indicated on the plot. (c) IPA canonical pathway analysis on the differentially expressed genes in RD vs. pCR. (d) IPA disease and function analysis of differentially expressed genes in RD vs. pCR. Orange color indicate activation while blue color indicates suppression.

Identification of mRNA Transcripts Predicative of RD and pCR
Given the limited number of differentially expressed genes using differential analysis, we subsequently performed receiver operating characteristic (ROC) analysis to identify the set of genes which can predict the response of TNBC patients to NAC (RD vs. pCR). ROC analysis identified 140 genes predictive of RD and 1490 genes predictive of pCR (asymptotic p-value < 0.05) as illustrated in the volcano plot ( Figure 2a and Table S3). C1orf116, CHST1, TP53INP2, RHOB, SLC9A3R2, EHF, BAD, FAM89B, PEBP1, and RHOC were the top 10 genes predictive of the RD phenotype. PPI network analysis of genes predictive of RD revealed a weak association ( Figure S1). Therefore, we employed an alternative strategy to identify potential essential genes for the RD phenotype by crossing the enriched genes in the RD from ROC analysis with the CRISPR-Cas9 pooled library functional screen data from the Achilles project in TNBC [8]. Data presented in Figure 2b shows the overlap between the enriched genes in RD and genes with significant effects on . TNBC patients were grouped according to pCR (n = 38) and residual disease (RD, n = 62) and were subjected to differential analysis (1.5 FC, p< 0.05). (a) Hierarchical clustering of TNBC as a function of pCR (n = 38) or RD (n = 62) based on differentially expressed mRNA genes. Each column represents one patient sample, and each row represents an mRNA. The expression level of each gene (log2) in a single sample is shown according to the color scale. (b) Volcano plot illustrating the upregulated (red) or downregulated (blue) genes in RD vs. pCR TNBC. The most differentially expressed mRNAs are indicated on the plot. (c) IPA canonical pathway analysis on the differentially expressed genes in RD vs. pCR. (d) IPA disease and function analysis of differentially expressed genes in RD vs. pCR. Orange color indicate activation while blue color indicates suppression.

Identification of mRNA Transcripts Predicative of RD and pCR
Given the limited number of differentially expressed genes using differential analysis, we subsequently performed receiver operating characteristic (ROC) analysis to identify the set of genes which can predict the response of TNBC patients to NAC (RD vs. pCR). ROC analysis identified 140 genes predictive of RD and 1490 genes predictive of pCR (asymptotic p-value < 0.05) as illustrated in the volcano plot ( Figure 2a and Table S3). C1orf116, CHST1, TP53INP2, RHOB, SLC9A3R2, EHF, BAD, FAM89B, PEBP1, and RHOC were the top 10 genes predictive of the RD phenotype. PPI network analysis of genes predictive of RD revealed a weak association ( Figure S1). Therefore, we employed an alternative strategy to identify potential essential genes for the RD phenotype by crossing the enriched genes in the RD from ROC analysis with the CRISPR-Cas9 pooled library functional screen data from the Achilles project in TNBC [8]. Data presented in Figure 2b shows the overlap between the enriched genes in RD and genes with significant effects on TNBC survival (effect < −0.5), revealing eight common genes (ELOB, SLC39A7, TIMM13, BANF1, NDUFS1, NDUFB7, TRAPPC5, and MVD). The gene effect of the eight RD gene signature on various TNBC cell models based on the CRISPR-Cas9 screen is illustrated in Figure 2c as a violin plot, where a lower score indicates that the TNBC cell line is highly dependent on the corresponding gene for survival. The eight-gene RD signature was subsequently subjected to survival analysis using mean expression for the eight genes in each subject, which revealed significant predictive power for this signature on RFS ( Figure 2d) and OS (Figure 2e) for basal and HER2+ breast cancer, but not for those with Luminal A and Luminal B tumors. TNBC survival (effect < −0.5), revealing eight common genes (ELOB, SLC39A7, TIMM13, BANF1, NDUFS1, NDUFB7, TRAPPC5, and MVD). The gene effect of the eight RD gene signature on various TNBC cell models based on the CRISPR-Cas9 screen is illustrated in Figure 2c as a violin plot, where a lower score indicates that the TNBC cell line is highly dependent on the corresponding gene for survival. The eight-gene RD signature was subsequently subjected to survival analysis using mean expression for the eight genes in each subject, which revealed significant predictive power for this signature on RFS ( Figure 2d) and OS (Figure 2e) for basal and HER2+ breast cancer, but not for those with Luminal A and Luminal B tumors.

Targeted Depletion of RD-Essential Genes Inhibited TNBC Colony Forming Unit (CFU) Potential and Enhanced Sensitivity to NAC
To validate the potential role of the identified eight genes in mediating RD to NAC, we used siRNA to suppress the expression of ELOB, SLC39A7, TIMM13, BANF1, NDUFS1, NDUFB7, TRAPPC5, and MVD in MDA-MB-231 cells and assessed their effects on CFU as single agent or in combination with doxorubicin and paclitaxel. Data presented

Targeted Depletion of RD-Essential Genes Inhibited TNBC Colony Forming Unit (CFU) Potential and Enhanced Sensitivity to NAC
To validate the potential role of the identified eight genes in mediating RD to NAC, we used siRNA to suppress the expression of ELOB, SLC39A7, TIMM13, BANF1, NDUFS1, NDUFB7, TRAPPC5, and MVD in MDA-MB-231 cells and assessed their effects on CFU as single agent or in combination with doxorubicin and paclitaxel. Data presented in Figure 3a illustrated significant inhibition of CFU following gene knockdown, which was most remarkable in BANF1, TIMM13, SLC39A7, MVD, NDUFB1, ECO13, TRAPPC5, and NDUFS1, respectively. Interestingly, the effects were further enhanced when BANF1, MVD, TIMM13, SLC39A7, and ECO13 were combined with doxorubicin. Similarly, the effects were most remarkable when BANF1, SLC39A7, TIMM13, and MVD were combined with paclitaxel. Similar results were also obtained when using the BT-549 TNBC models ( Figure S2). in Figure 3a illustrated significant inhibition of CFU following gene knockdown, which was most remarkable in BANF1, TIMM13, SLC39A7, MVD, NDUFB1, ECO13, TRAPPC5, and NDUFS1, respectively. Interestingly, the effects were further enhanced when BANF1, MVD, TIMM13, SLC39A7, and ECO13 were combined with doxorubicin. Similarly, the effects were most remarkable when BANF1, SLC39A7, TIMM13, and MVD were combined with paclitaxel. Similar results were also obtained when using the BT-549 TNBC models ( Figure S2).  The results were further confirmed using acridine orange and ethidium bromide viability staining in the MDA-MB-231 which revealed induction of cell death in siRNA knockdown cells and was further enhanced when combined with doxorubicin ( Figure 4), therefore corroborating a role for these genes in TNBC survival and sensitive to doxorubicin.

Discriminant Analysis
To determine the potential predictive value of the identified genes, we perfo Orthogonal Projections to Latent Structures Discriminant Analysis (OPLS-DA) on th covery cohort. Data from the OPLS-DA model is presented as score plots to visuali differences between RD and pCR ( Figure 5a). The loading plot in Figure 5b indica variables that express the observed difference between RD and pCR, where each va is colored according to density score. The score contribution of each variable is sho Figure 5c. The variable influence on projection (VIP) predicted score for each vari listed in Table S4. The model was tested on the same discovery cohort which revea AUC = 0.93 to discriminate RD from pCR ( Figure S3) and was subsequently valida

Discriminant Analysis
To determine the potential predictive value of the identified genes, we performed Orthogonal Projections to Latent Structures Discriminant Analysis (OPLS-DA) on the discovery cohort. Data from the OPLS-DA model is presented as score plots to visualize the differences between RD and pCR (Figure 5a). The loading plot in Figure 5b indicates the variables that express the observed difference between RD and pCR, where each variable is colored according to density score. The score contribution of each variable is shown in Figure 5c. The variable influence on projection (VIP) predicted score for each variable is listed in Table S4. The model was tested on the same discovery cohort which revealed an AUC = 0.93 to discriminate RD from pCR ( Figure S3) and was subsequently validated on a second cohort of TNBC (n = 50) revealing an AUC = 0.75 to discriminate RD from pCR in this independent validation cohort (Figure 5d). a second cohort of TNBC (n = 50) revealing an AUC = 0.75 to discriminate RD from pCR in this independent validation cohort (Figure 5d).

PPI Network Analysis for Genes Predictive of pCR
In order to gain more insight in to the functional and network enrichment in pCR, the list of common genes predicative of pCR from the discovery and validation cohorts using univariate analysis were imported into the STRING database and were subjected to network analysis. Data presented in Figure 6a illustrates a dense network with multiple interactions. Interestingly, the highest network enrichment score was for functional categories related to immune and cell cycle regulation (Figure 6b and Table S5). Given the apparent role of immune signatures in predicting pCR, we identified an 87-immune signature base on GO analysis, which was subjected to separate STRING analysis revealing activation of various immune processes such as innate immunity, inflammatory response, myeloid leukocyte activation, adaptive immune response, activation of neutrophils and natural killer cells as the most affected functional categories (Figure 7 and Table S6).

PPI Network Analysis for Genes Predictive of pCR
In order to gain more insight in to the functional and network enrichment in pCR, the list of common genes predicative of pCR from the discovery and validation cohorts using univariate analysis were imported into the STRING database and were subjected to network analysis. Data presented in Figure 6a illustrates a dense network with multiple interactions. Interestingly, the highest network enrichment score was for functional categories related to immune and cell cycle regulation (Figure 6b and Table S5). Given the apparent role of immune signatures in predicting pCR, we identified an 87-immune signature base on GO analysis, which was subjected to separate STRING analysis revealing activation of various immune processes such as innate immunity, inflammatory response, myeloid leukocyte activation, adaptive immune response, activation of neutrophils and natural killer cells as the most affected functional categories (Figure 7 and Table S6).

Survival Analysis of 87-Immune-Gene Signature in an Independent Breast Cancer Cohorts
The identified 87-immune-gene signature was assessed for association with OS, RFS, and DMFS in additional cohorts from the KMplotter database, a commonly used database for cancer survival analysis. Interestingly, the 87-immune-gene signature was highly pre-

Discussion
The quest to identify potential biomarkers predictive of TNBC response to NAC have been a daunting task, that nevertheless remains important to pursue for patient triage and resource prioritization. In the current study we took multiple approaches to define gene signatures predicative of TNBC pCR and RD and the potential targeting of RD-associated genes to enhance the efficacy of standard NAC. Our data identified 140 genes to be predicative of RD, while the bulk of differentially expressed genes were found to be associated with pCR. This finding implies that genes associated with RD are likely driven by tumor cells, with each patient having their own oncogenic signatures. On the other hand, the signature predictive of pCR was mostly driven by immune infiltrating cells; therefore, it shows higher predicative power and was largely present across patients from the discovery and validation cohorts. In our quest to identify potential RD-essential genes, we integrated our computational and prediction data with CRISPR-Cas9 functional screen data from the Achilles project, leading to the identification of eight essential genes for the RD phenotype. The most catastrophic effects of gene knockdown were for BAF Nuclear Assembly Factor 1 (BANF1), Translocase of Inner Mitochondrial Membrane 13 (TIMM13), Solute Carrier Family 39 Member 7 (SLC39A7), and Mevalonate Diphosphate Decarboxylase (MVD). A recent report suggested a correlation between elevated expression of BANF1 and lymph node metastasis status and TNM staging in TNBC [9]. In addition, BANF1 has been shown to regulate PARP1-directed DNA damage response to oxidative

Discussion
The quest to identify potential biomarkers predictive of TNBC response to NAC have been a daunting task, that nevertheless remains important to pursue for patient triage and resource prioritization. In the current study we took multiple approaches to define gene signatures predicative of TNBC pCR and RD and the potential targeting of RDassociated genes to enhance the efficacy of standard NAC. Our data identified 140 genes to be predicative of RD, while the bulk of differentially expressed genes were found to be associated with pCR. This finding implies that genes associated with RD are likely driven by tumor cells, with each patient having their own oncogenic signatures. On the other hand, the signature predictive of pCR was mostly driven by immune infiltrating cells; therefore, it shows higher predicative power and was largely present across patients from the discovery and validation cohorts. In our quest to identify potential RD-essential genes, we integrated our computational and prediction data with CRISPR-Cas9 functional screen data from the Achilles project, leading to the identification of eight essential genes for the RD phenotype. The most catastrophic effects of gene knockdown were for BAF Nuclear Assembly Factor 1 (BANF1), Translocase of Inner Mitochondrial Membrane 13 (TIMM13), Solute Carrier Family 39 Member 7 (SLC39A7), and Mevalonate Diphosphate Decarboxylase (MVD). A recent report suggested a correlation between elevated expression of BANF1 and lymph node metastasis status and TNM staging in TNBC [9]. In addition, BANF1 has been shown to regulate PARP1-directed DNA damage response to oxidative stress [10]. Thus far, TIMM13 has not been linked to cancer, so our data suggest a plausible role for this gene in TNBC response to NAC. Knockdown of SLC39A7 was show to inhibit cell growth and induce apoptosis in colorectal cancer (CRC) [11]. Recent data suggested SLC39A7 to promote glioma cell malignancy through the TNF-α-mediated NF-κB pathways [12]. Elevated expression of SLC39A7 was also reported in tamoxifen-resistant breast cancer cells [13]. However, our data are the first to implicate SLC39A7 in the response of TNBC to NAC. Among the identified RD genes, MVD has not been linked to human cancers thus far. Elongin B (ELOB) knockdown exhibited minimal effects on TNBC viability, which was significantly enhanced when combined with doxorubicin and paclitaxel. ELOB is one subunit of the transcription factor B complex. Although there is currently no report connecting ELOB with human cancers, our data suggest its potential use to enhance TNBC sensitivity to NAC.
Interestingly, we observed the presence of an immune gene signature as a common feature in TNBC who responded to NAC in two independent cohorts. In particular, we identified an 87 gene signature that predicted better OS, RFS, and DMFS in basal and HER2+ breast cancer. Our data are consistent with recently published reports [14,15] and with our recently published work on single-cell analysis of TNBC in the context of TNBC response to NAC [7]. This interesting observation raises a key question about whether TNBC tumors with heavy immune infiltration represent different molecular subtype with different sensitivity to NAC compared to tumors lacking immune infiltration, which warrants further investigation. Alternatively Of particular interest, SLAMF7 has recently been shown to program T cells toward exhaustion phenotype via SLAMF7+ tumor-associated macrophages (TAMs) in the tumor microenvironment of renal cell carcinoma [16]. A recent study suggested expression of SLAMF7 to promote cytotoxicity of human NK cells [17]. IFNG was also among the identified genes predicative of pCR. Concordantly, our recent single-cell analysis also identified activation of the IFNG pathway in TNBC patients exhibiting pCR. Taken together, our data have highlighted activation of immune cells as the hallmark associated with pCR.
One limitation of the current study is the relatively small sample size of the validation cohort. Nonetheless, only expression data were available from the validation cohort. Interestingly and despite the differences in bioinformatics tools and reference genome used for the discovery and validation cohorts, the gene signature showed reasonable performance in the validation cohort (AUC = 0.75). However, those findings remain to be validated in additional cohorts and in TNBC patients receiving other NAC regimens.
Our results identified gene panels predictive of RD and pCR in TNBC, where the strongest prediction was for pCR. Eight genes were identified as TNBC-essential genes, which were highly predicative of OS and RFS, and were essential for TNBC proliferation. Interestingly, pCR-associated genes were enriched for immune cell functions and were predictive for better OS, RFS, and DMFS in basal and to a lesser extent in HER2+ breast cancer patients. Our findings have potential clinical implications in patient stratification and tailored therapy.

Study Cohorts
RNA-Seq data were retrieved from the PRJNA688066 [14] (discovery) and GSE192341 [15] (validation cohort): Clinical characteristics for patients from both cohorts are shown in Table 1.

Data Retrieval and Bioinformatics
The Kallisto index was generated by creating a de Bruijn graph using the reference transcriptome GENCODE release (V33) and a k-mer length of 31. FASTQ files from the PRJNA688066 dataset were mapped and aligned to the generated GENCODE release (V33) index using Kallisto v0.46 as described before [18,19]. Expression values (TPM) were subjected to differential, principal component, and clustering analysis in AltAnalyze as described before [20,21]. Expression data from the validation cohort was retrieved from the gene omnibus (GEO) dataset GSE192341.

Statistical Analysis
ROC analysis was employed to identify gene sets predictive of RD and pCR using SPSS version 26. The volcano plot was generated in GraphPad prism (v9.0) by plotting the AUC on the y-axis and −log10 p-value on the x-axis. Other plots were also generated in GraphPad prism v9.0.

Ingenuity Pathway Analysis
Differentially expressed genes in RD vs. pCR were imported into the IPA software (Ingenuity Systems; www.ingenuity.com/ accessed on 4 August 2022).) and were used for functional annotations and network analysis using canonical, upstream regulator, and downstream effector analyses The p-value is the negative log of P and represents the possibility that focus genes in the network being found together by chance [7,22].

Discriminant Analyses
To determine the ability of predicted variables (genes identified from the ROC analysis) to discriminate between RD and pCR, discriminant analyses were performed using an OPLS-DA classifier using Soft Independent Modelling by Class Analogy (SIMCA) software (version 16; Umetrics, Sweden) on the discovery dataset (n = 90). The model was then tested on the validation dataset (n = 50) and the performance was assessed by generating ROC curve and determining the AUC value. The sensitivity and specificity constants of the test were determined based on similar classification scores by OPLS-DA.

Protein-Protein Interaction (PPI)
PPI network analysis was conducted using the STRING database as described before [23]. The 500 identified common genes or the 87-immune-gene signatures predictive of pCR were subjected to PPI in the STRING database v 11.5 (https://string-db.org/ accessed on 4 August 2022) as described before [24,25].

CRISPR-Cas9 Screen Data Retrieval
Genome-wide CRISPR-Cas9 screen data were retrieved from the Achilles project and gene effect score of ≤−0.5 for TNBC cell lines were included in the analysis [8].

Gene Silencing Using SiRNA
The scrambled siRNA control and ON-TARGETplus SMARTpool siRNA targeting human ELOB, SLC39A7, TIMM13, BANF1, NDUFS1, NDUFB7, TRAPPC5, and MVD were obtained from Dharmacon (Lafayette, CO, USA). Transfection was carried out using a reverse transfection strategy as previously described [19]. Briefly, siRNA was diluted in 50 µL of Opti-MEM (GIBCO, Carlsbad, CA, USA), and 1 µL of Lipofectamine 2000 (Invitrogen) was diluted in 50 µL of OPTI-MEM. The diluted siRNA (at a final concentration of 30 nM) and Lipofectamine 2000 were mixed together and allowed to form complexes at room temperature for 20 min. For transfection, twenty microliters of transfection mixture were first placed in individual wells of the tissue culture plate, and then 10,000 of MDA-MB-231 or BT-549 TNBC cells in 60 µL transfection medium (complete DMEM without Pen-Strep) were added to each well. The transfection cocktail was then replaced with complete DMEM one day later.

Colony Forming Unity (CFU) Assay and Chemotherapy Sensitivity of TNBC Models
Transfected cells were cultured for 48 h and were then exposed to paclitaxel (10 nM). Three days later, cells were fixed for 5 min using 4% PFA and then were washed twice using PBS followed by staining with crystal violet (0.1% in 10% EtOH) at room temperature for 10 min. The images were then captured and compared to experiment controls. The plates were then allowed to air dry, followed by CFU quantification of dissolved crystal violet in 5% SDS. Absorbance was measured at 590 nm. The experiments were repeated at least twice, and data were presented as the mean ± SD from four replicas.

Detection of Cell Death Using Fluorescence Microscopy
AO/EtBr fluorescence microscopy was used to assess apoptotic cells in TNBC cells after gene knockdown as single treatment or in combination with paclitaxel in TNBC cells, as we previously described [26]. Briefly, TNBC cells were transfected in a 24-well flat-bottom tissue culture plate. Forty eight hours later, media was changed and paclitaxel was added at a final concentration of 7.5 nM. Two days later, TNBC cells were washed twice using PBS and were subsequently stained with dual fluorescent staining solution (100 µg/mL AO and 100 µg/mL EtBr) for 2 min. AO and EtBr were purchased from Sigma (Sigma Aldrich, St. Louis, MO, USA);Stained ells were then visualized under inverted microscope (Olympus IX73 fluorescence microscope, Olympus, Tokyo, Japan). AO staining was used to visualize the number of cells that had undergone apoptosis, while EtBr-positive cells indicated necrotic cells or late apoptotic cell, with the later also showing condensed chromatin.

Conclusions
Our data have identified gene signatures predicative of RD and pCR in TNBC with potential clinical implications.