Introduction

Triple negative breast cancer (TNBC) is an aggressive form of breast cancer that affects 10–20% of all breast cancer patients and is characterized by its lack of expression of estrogen, progesterone and HER2 receptors1. The standard of care for TNBC patients primarily relies on conventional anthracycline and taxane-based chemotherapy regimens, and few next-generation therapies have shown efficacy in patients with this disease2. Paclitaxel, a taxane-based chemotherapeutic commonly used in TNBC treatment3, targets microtubules to disrupt the formation of the mitotic spindle, resulting in cell cycle arrest and apoptosis. Paclitaxel treated cells often fail to achieve symmetric division of nuclear content, potentially resulting in apoptosis4. Cells that persist through paclitaxel treatment may exhibit nuclear damage in the form of either disrupted nuclear structures and/or acquired polyploidy5,6. While 22% of TNBC patients treated with paclitaxel achieve pathological complete response, the outcome for those with residual disease is relatively poor7,8. Moreover, paclitaxel monotherapy only achieves a median 5.5 month progression free survival in patients with locally advanced or metastatic disease9. Therefore, there is a need to better understand the molecular basis of paclitaxel response and mechanisms of resistance that may be targeted for therapeutic benefit.

Phenotypic plasticity enables malignant cells to rapidly adapt to therapeutic challenge10 and can also drive acquired drug resistance11. Adaptive responses often involve activation of new transcription factors that in turn upregulate programs that repress immune activation12, grant tolerance to DNA replication stress13, or enable evasion of apoptosis 14. Single-cell RNA sequencing (scRNA-seq) is a powerful approach to investigate the subtle but critical differences in transcriptional landscape that distinguish cellular phenotypic states and to identify molecular programs associated with different therapeutic sensitivities15,16,17.

To elucidate adaptive responses of TNBC cells to paclitaxel, we performed deep single-cell RNA sequencing of HCC1143 TNBC cells before and after paclitaxel treatment. Paclitaxel induced a range of phenotypic changes, including altered cell cycle phase distribution, increased proportion of multinucleated cells, increased expression of senescence and DNA damage associated biomarkers, and upregulation of interferon-related gene programs. Comparison of gene expression profiles from paclitaxel treated versus IFNB or IFNG treated cells enabled identification of genes that were uniquely upregulated after paclitaxel treatment, including a suite of transcription factors. Functional assessment with siRNA knockdown confirmed that many of these TFs are critical for mediating resistance to paclitaxel. Using live-cell imaging, we probed the temporal dynamics of these functional responses, which demonstrated that knockdown of ELF3, FOSL1 and IRF9 synergize with paclitaxel to slow cell cycle progression. Together, these analyses identify upregulation of ELF3, FOSL1 and IRF9 as important regulators of cell cycle progression that mediate paclitaxel response, and which may serve as biomarkers of response.

Results

Paclitaxel modulates multiple cancer-associated phenotypes

We identified changes in growth and proliferation induced by paclitaxel by performing an EdU incorporation assay for three TNBC cell lines (HCC1143, HCC1806 and MDA-MB-468) treated for 72 h with dose escalation of paclitaxel (0.01 nM-81 nM). Higher paclitaxel doses (1 nM-81 nM) induced increasingly frequent nuclear aberrations, including multi-lobular and multinucleated phenotypes (Fig. 1A). Comparison of cell counts at 72 h to 0 h and DMSO control revealed a sigmoidal relationship between paclitaxel dose and growth rate for all three cell lines (Fig. 1B), where the highest concentrations were cytostatic (growth rate plateau >  = 0) for two of the three cell lines18. The total DAPI intensity and total EdU intensity for each cell was calculated and thresholds were used to classify cells into 2N/EdU-, EdU + or 4N/EdU- cell cycle states (Supplemental Fig. 1A). There was a significant reduction in the fraction cells staining positive for EdU (Fig. 1C), indicating that paclitaxel treatment reduces the proportion of cells that are actively synthesizing new DNA, as expected6,19.

Fig. 1
figure 1

Paclitaxel modulates multiple cancer associated phenotypes. (A) Representative fluorescent images showing EdU incorporation for HCC1143 cells treated with DMSO or Paclitaxel at the indicated doses for 72 h. Merged image shows DAPI (blue channel, marker for DNA), EdU (red channel, marker for active DNA synthesis) and HCS Cellmask (green channel, marker used as cytoplasmic segmentation aid). Monochrome images below show DAPI and EdU channels separately, and magenta outline indicates cytoplasmic boundary from segmentation. (B) Growth rates calculated for 72 h paclitaxel treatment for three TNBC cell lines: HCC1143, HCC1806, MDA-MB-468. Dots indicate individual measurements, line indicates loess fit of the points. Growth rate was calculated from phase confluency using an Incucyte S3. (C) Barplot indicating proportion of EdU + cells for each cell line at each paclitaxel concentration. Significance assessed with Dunnett’s test. (D) Representative fluorescent images showing HCC1143 cells treated with either 0.1% DMSO or 1 nM Paclitaxel for 72 h and then stained with DAPI (blue) and TUBB3 (red). Monochrome images show DAPI and TUBB3 separately, where magenta outline indicates cytoplasmic segmentation boundary. (E–F) Barplots showing mean TUBB3 (E) and p16/p15 (F) staining intensity for triplicate wells of HCC1143 treated with Paclitaxel dose escalation, after normalization to paired DMSO control (horizontal line). Mean intensity was calculated for the entire cell segmentation area (including nuclear and cytoplasmic regions). Significance assessed with Dunnett’s test. (G) Barplot comparing the fraction of cPARP positive cells for mononucleated (magenta) versus multinucleated (cyan) cells within the same treatment condition. cPARP positive threshold was set to the 99th quantile of DMSO treated cells total cPARP nuclear intensity (Supplemental Fig. 1C). Significance assessed with proportions test. For all statistics: * = p < 0.05, ** = p < 0.01, *** = p < 0.001.

Full size image

We further assessed adaptive cellular responses by analyzing biomarkers associated with senescence (p16/p15), DNA damage induced apoptosis (cPARP), and microtubule components (TUBB3) for HCC1143 cells treated with paclitaxel dose escalation for 72 h (Fig. 1D). TUBB3 overexpression has been associated with resistance to multiple microtubule targeting drugs, and consistent with this, we found a dose-dependent relationship between TUBB3 expression and paclitaxel concentration (Fig. 1E)20,21. There was also a positive association between both cytoplasmic and nuclear p16/p15 staining and paclitaxel dose (Fig. 1F). Additionally, we observed a strong correlation between p16/p15 and TUBB3 expression at the single cell level across paclitaxel concentrations, suggesting that the TUBB3 highly expressing cells represent a senescent subpopulation of cells (Supplemental Figs. 1B,C, Pearson correlation = 0.70, r^2 = 0.5). Increasing doses of paclitaxel induced a corresponding increase in the fraction of cPARP positive cells (DMSO: 6%, 81 nM Paclitaxel: 28% cPARP positive), indicating induction of DNA damage driven apoptosis. Higher paclitaxel doses resulted in a significantly higher proportion of mononucleated cells staining positive for cPARP as compared to multinucleated (19.4% mononucleated cells and 7.1% multinucleated cells cPARP positive at 81 nM paclitaxel, proportions test p = 0.017), suggesting that multinucleated cells are less likely to undergo DNA damage-induced apoptosis (Fig. 1G). Comparison between these findings and the EdU results showed that across a range of paclitaxel doses, multinucleated cells tended to belong to the 4N/EdU- cell cycle state (Supplemental Fig. 1D). Together this suggests that multinucleated cells that survive paclitaxel treatment are cell cycle arrested and less likely to undergo DNA damage-induced apoptosis than are mononucleated cells.

Cells surviving paclitaxel treatment halt cycling and upregulate interferon response genes

We assessed paclitaxel-induced molecular programs with 10X Genomics single-cell whole transcriptome sequencing of HCC1143 cells treated with either DMSO vehicle control or 1 nM paclitaxel for 24 h or 72 h (Fig. 2A). After quality control filtering that required cells to have a minimum of 3000 unique genes and a maximum of 25% mitochondrial counts, we recovered 3194 total cells (513 – 1106 cells per condition) with a mean UMI count of 63,668 (Supplemental Fig. 2A).

Fig. 2
figure 2

Cells surviving paclitaxel treatment have halted cycling and upregulated interferon response genes. (A) UMAP color coded by treatment condition. DMSO_24 = 0.1% DMSO for 24 h, DMSO_72 = 0.1% DMSO for 72 h, PTX_24 = 1 nM Paclitaxel for 24 h, PTX_72 = 1 nM Paclitaxel for 72 h. (B) Barplot showing proportion of each condition assigned to G1, S, or G2M cell cycle state based on transcriptional profile. (C) Volcano plot of differentially expressed genes for Paclitaxel treatment versus DMSO at 24 h. (E) Significantly differentially expressed genes were determined with cutoffs of Benjamini Hochberg corrected p < 0.05 and absolute Log2FoldChange > 0.5. (D) Reactome pathway enrichment results for genes significantly upregulated after paclitaxel treatment at 24 h. Size indicates the number of genes upregulated within the pathway, color indicates significance. (E) Volcano plot of differentially expressed genes for Paclitaxel treatment versus DMSO at 72 h. Differentially expressed genes determined as in (C). (F) Reactome pathway enrichment results for genes significantly upregulated after paclitaxel treatment at 72 h. Size indicates the number of genes upregulated within the pathway, color indicates significance.

Full size image

We examined drug-induced changes in cell cycle distribution by assigning cell cycle status to each individual cell using aggregate expression of canonical gene programs for S and G2/M18,19. In agreement with the results of the EdU incorporation assay (Fig. 1C), we observed a reduction in the fraction of S phase cells after paclitaxel treatment as compared to time-matched vehicle control (Fig. 2B). UMAP visualization showed that the cells organized overall by treatment condition and cell cycle phase (Supplemental Fig. 2B).

We analyzed time-matched conditions to identify significantly differentially expressed genes induced by paclitaxel treatment (Wilcoxon rank sum test, absolute log2 fold-change > 0.5, Benjamini Hochberg FDR < 0.01). This revealed a time-dependent change in molecular programs with 66 significantly upregulated and 57 significantly downregulated genes after 24 h of paclitaxel treatment, and 256 significantly upregulated genes and 58 significantly downregulated genes after 72 h (Fig. 2C, Supplemental Fig. 2C). Reactome pathway enrichment analysis revealed that the significantly upregulated genes from the 24 h paclitaxel treated sample were enriched for multiple programs related to Interferon Signaling and Translation (Fig. 2D, Supplemental Data 2). Programs uniquely upregulated after 72 h of paclitaxel treatment include Response to Chemical Stress, Cell Cycle Progression, and Antigen Processing-Cross presentation (Fig. 2E-F). The ontologies enriched after 72-h paclitaxel treatment had low overlap with those at 24 h, suggesting that paclitaxel response is a dynamic process with transcriptional changes that continue to emerge even after the first 24 h of treatment (Jaccard Index = 0.023, Supplemental Fig. 2D). Notably, the Neutrophil Degranulation pathway was significantly enriched at both time points, with upregulated genes related to antigen presentation (HLA-B, HLA-C, B2M) and differentiation (CD47, CD55, CD59, CD63). Paclitaxel treatment also induced significant upregulation of the pro-tumorigenic chemokines CXCL1 and CXCL8 (Supplemental Fig. 2E)22,23,24,25. Together this shows that TNBC cells that survive paclitaxel treatment have altered surface marker expression and produce tumor supportive chemokines.

Paclitaxel response activates canonical interferon response genes

Despite gene enrichment consistent with interferon response, the paclitaxel treated cells showed no evidence of autocrine interferon signaling, indicating that paclitaxel induces interferon response pathways in a non-canonical manner (Supplemental Fig. 3A). To disentangle the paclitaxel response signature from an interferon-induced response, we performed a second scRNA-seq experiment with HCC1143 cells after treatment for 72 h with 7 perturbations that target ligand-receptor pairs known to play an important role in normal and pathological breast tissue26,27: Interferon Beta (IFNB), Interferon Gamma (IFNG), Transforming Growth Factor Beta (TGFB), Oncostatin M (OSM), Lymphotoxin Alpha (LTA), Notch Inhibitor (NOTCHi) and combination of Notch Inhibitor and Interferon Beta (NOTCHi_IFNB). Cells were treated for 72 h and then harvested and sequenced with the 10X Genomics scRNA-seq pipeline. After quality control filtering, we recovered 4231 total cells (295–725 cells per condition, Supplemental Fig. 3B).

Fig. 3
figure 3

Paclitaxel response activates canonical interferon response genes. (A) UMAP showing the scRNA-seq landscape after ligand perturbations. IFNB = Interferon-Beta, OSM = Oncostatin-M, NOTCHi_IFNB = Notch inhibitor + Interferon-Beta, NOTCHi = Notch inhibitor, TGFB = Transforming Growth Factor Beta, IFNG = Interferon-Gamma, LTA = Lymphotoxin-Alpha, PBS = Phosphate Buffered Saline (control). (B) Heatmap showing the Pearson correlation for log2 fold-change values for each perturbation versus time-matched control. (C,D) Gene enrichment map for Paclitaxel uniquely upregulated (3C) and Paclitaxel + Interferon shared upregulated (3D) genes. Color indicates significance, size indicates number of upregulated genes, and lines connect ontologies with shared elements. (E) ChEA3 transcription factor enrichment ranks computed from 140 Paclitaxel uniquely upregulated genes (x axis) versus 120 Paclitaxel-Interferon shared upregulated genes (y axis). Lower rank indicates higher imputed activity, and named TFs with red dots were within the top 15 ranks for either geneset. TFs to the lower right of the diagonal have higher imputed activity within the PTX + IFN shared upregulated gene set, and TFs to the upper left of the diagonal have higher imputed activity within the PTX uniquely upregulated gene set. (F) Bar plot showing Average Log2FC from paclitaxel treated scRNA-seq data for the 24 top ranked transcription factors (intersect of top 15 ranked for PTX unique or PTX shared individually). Red font indicates TFs significantly upregulated (average log2 fold-change > 0.25, FDR < 0.01) at either 24 or 72 h of paclitaxel treatment compared to vehicle control.

Full size image

UMAP projection of scRNAseq profiles revealed that ligand-treated cells largely grouped by perturbation (Fig. 3A) and cell cycle state (Supplemental Fig. 3C). The IFNB, IFNG, TGFB, NOTCHi and NOTCHi_IFNB conditions all had an increase in proportion of G1 cells compared to control, suggesting these ligands are cytostatic in this cell line (Supplemental Fig. 3D). Based on the observation that paclitaxel induced interferon related pathways, we next sought to evaluate the similarity in transcriptional response between paclitaxel and ligand perturbations. To that end, we computed the differential expression of all genes for each perturbation compared to time-matched vehicle control and then evaluated the pairwise Pearson correlation of log2 fold-change values (Fig. 3B). The IFNB and IFNG conditions were the most strongly correlated (Pearson correlation = 0.86), indicating a conserved impact on transcription despite acting through different receptors. The 72 h paclitaxel condition was highly correlated with both interferon treatments (IFNB Pearson correlation = 0.57, IFNG Pearson correlation = 0.47) as compared to the other single-agent perturbations (0.0, 0.11, 0.18, 0.38 Pearson correlation with OSM, LTA, NOTCHi and TGFB, respectively). Despite the strong transcriptional overlap between interferon treated and paclitaxel treated samples, we observed no multinucleation following interferon treatment (Supplemental Fig. 3E).

While type 1 and type 2 interferons primarily exhibit antitumor effects through activation of the immune system, some studies have shown that they have direct effects through induction of cell cycle arrest or apoptosis in malignant cells28,29. To further clarify the overlapping transcriptional responses of paclitaxel and interferon, we next sought to differentiate between pathways that were uniquely induced by paclitaxel response or that represent common responses induced by paclitaxel or interferon. Reactome pathway enrichment analysis revealed that the 140 genes upregulated after paclitaxel treatment but not after IFNG or IFNB (“paclitaxel-unique”) were enriched in molecular programs related to wound healing, protein folding, and intrinsic apoptotic signaling pathway (Fig. 3C), whereas the 117 genes upregulated by all three treatments (“paclitaxel-shared”) were associated with defense response to virus and antigen presentation (Fig. 3D).

As paclitaxel has been shown to activate cGAS/STING and subsequently induce upregulation of Interferon related genes in human breast cancer, we hypothesized that the strong overlap in interferon and paclitaxel transcriptional responses was driven by a shared increase in transcription factor (TF) activity5,30,31. We identified enriched TFs in our paclitaxel-unique and paclitaxel-shared gene signatures using ChEA332, which evaluates the expression of gene targets downstream from a TF of interest. DDIT3, JUN, KLF6, and ATF3 emerged as enriched transcription factors across both the paclitaxel-unique and paclitaxel-shared gene signatures (Fig. 3E). The paclitaxel-unique genes were enriched for TFs in the Immediate-Early Gene family, including JUN (JUN, JUNB, JUND) and FOS (FOS, FOSL1, FOSB)33. TFs enriched from the shared gene list were associated with activity of Interferon Regulatory Factors (IRF1/IRF7/IRF9) and Basic Leucine Zipper family (BATF2/BATF3), both related to antiviral response and regulation of antigen-presenting cells34,35,36. The high activity of IRF7 is consistent with activation of the cytosolic nucleotide sensor RIGI, suggesting that the nuclear damage induced by paclitaxel drives an increase in cytosolic RNA or DNA 37.

Inhibition of paclitaxel-induced transcription factors alters proliferation and nuclear morphology

We used siRNA knockdown in three basal-like TNBC cell lines (HCC1143, HCC1806, MDA-MB-468) to functionally assess prioritized TFs implicated in modulating response to paclitaxel. We nominated a panel of 13 TFs for functional testing, based on ChEA3 analysis and change in gene expression after 24 or 72 h of paclitaxel treatment (Fig. 3F, Log2FC > 0.25, Benjamini–Hochberg FDR < 0.01). The expression of the implicated TFs were not correlated at the single-cell level (Supplemental Fig. 3F), but several were upregulated after interferon treatment (Supplemental Fig. 4A). Most of the TFs included in this panel were either subunits of (ATF3, FOSL1, JUN, JUNB, MAFF)38 or known interactors with (ELF3, IRF7, DDIT3, NFE2L2)39,40,41,42 the AP-1 transcription factor family. Dysregulation of the AP-1 pathway is associated with multiple tumorigenic phenotypes including enhanced cellular growth, proliferation, and survival43. To functionally assess the role of these TFs in paclitaxel response, cells were transfected with siRNA for 24 h, then treated for 72 h with paclitaxel or DMSO, and subsequently fixed and stained with DAPI (nuclear marker). The resultant images were subjected to quantitative image analysis to identify nuclear and cellular masks, followed by assessment of total cell number and fraction of multinucleated cells for each condition.

Fig. 4
figure 4

Inhibition of paclitaxel-induced transcription factors alters proliferation and nuclear morphology. (A-B) Heatmaps showing percent inhibition of cell count (A) and proportion of multinucleated cells (B). Cell count is normalized to the same cell line DMSO + siNonTarget control, and color indicates mean value across triplicate. Significance for percent inhibition computed with two-tailed Student’s T-test with Bonferroni correction, and significance for fraction multinucleated computed with proportions test with Bonferroni correction. siNonTarget and siGAPD represent positive growth controls, and siPLK1 and siKIF11 represent negative growth controls. (C) Principal Component Analysis results for each siRNA knockdown where each combination of cell line (HCC1143, HCC1806, MDA-MB-468), feature (relative cell count, fraction multinucleated) and condition (DMSO, PTX) is considered a feature. (D) The Euclidean feature-distance from NonTarget control for each siRNA. For all statistics: * = p < 0.05, ** = p < 0.01, *** = p < 0.001.

Full size image

First, we analyzed the influence of TF knock-down on cell count after 72 h to evaluate their effects on cell viability. We found that knockdown of 3 of 13 TFs (NFE2L2, IRF7, MAFF) in the absence of paclitaxel significantly reduced cell numbers for at least one cell line (Student’s t-test, p < 0.05, Fig. 4A, upper). We then examined the influence of TF knock-down in the presence of paclitaxel to test our hypothesis that upregulation of these TFs mediates adaptive resistance. Knockdown of 5 of 13 TFs (NFE2L2, ELF3, IRF7, FOSL1, PLSCR1) in combination with paclitaxel significantly lowered cell count in at least one cell line compared to paclitaxel alone (Student’s t-test, p < 0.05, Fig. 4A, lower).

We hypothesized that these 13 TFs may also be involved in cytokinesis, based the changes we observed in nuclear morphology following paclitaxel treatment (Fig. 1A). We found that siRNA knockdown alone caused a significant increase in the fraction of multinucleated cells for 6 of 13 TFs for at least one cell line: NFE2L2, ELF3, SP100, FOSL1, MAFF, and ATF3 (Proportions test, p < 0.05. Figure 4B). Additionally, knockdown for 11 of 13 TFs in the presence of paclitaxel resulted in significantly increased fraction of multinucleated cells as compared to paclitaxel alone, for at least one cell line. These findings suggest an important role for these transcription factors in maintaining nuclear structure and achieving symmetric cytokinesis in proliferating breast cancer cells.

We next sought to identify which transcription factors within this panel induced the largest overall impact on cell phenotype across treatments and cell lines. Here, we considered each siRNA an independent sample and each combination of cell line (HCC1143, HCC1806, MDA-MB-468), treatment (DMSO, paclitaxel) and phenotype (relative cell count, fraction multinucleated) as 12 independent features. Principal component analysis applied to the transformed data revealed that our positive and negative growth controls separated along Component 1 (Fig. 4C). To identify the TFs that had the greatest overall impact on phenotype, we computed the Euclidean feature-distance (distance for z-scored features) to identify TFs that induced the greatest feature-distance from siNonTarget positive growth control. We found that ELF3, FOSL1, and NFE2L2 knockdown had the largest Euclidean feature-distance, indicating that knockdown had a strong effect on both proliferation and regulation of nuclear morphology across the three TNBC cell lines (Fig. 4D). Protein quantification for ELF3, FOSL1 and NFE2L2 confirmed that each siRNA was effective at reducing protein levels for the target transcript (Supplemental Figs. 4A-B).

ELF3 mediates cell cycle progression under paclitaxel treatment

Motivated by the observation that many anti-cancer drugs act by targeting the cell cycle, we next explored the influence of prioritized TFs on cell cycle progression by leveraging HCC1143 cells genetically engineered to express a cell cycle reporter. The cell cycle state of HDHB-mClover/NLS-mCherry HCC1143 cells can be determined by quantification of relative HDHB-mClover (nuclear translocating cell cycle reporter) intensity within the cytoplasm compared to the nuclear signal marked by NLS-mCherry (stable nuclear localization)44,45. Cells in G1 cell cycle phase have near-equal nuclear and cytoplasmic HDHB-mClover expression, cells in S/G2 cell cycle phases exclude the HDHB-mClover from the nucleus, and cells in M phase concentrate the HDHB-mClover expression to the nucleus. Here we focused on NFE2L2, ELF3 and FOSL1, which induced the largest phenotypic effects; we additionally tested IRF9 which has been previously linked to anti-microtubule chemotherapy resistance 46. Reporter cells were subjected to siRNA transfection for 24 h and then treated with either 1 nM paclitaxel or DMSO. Treated cells were imaged every 15 min for 72 h. Nuclear and cytoplasmic masks were segmented with custom trained Cellpose47 models and the resultant data used to classify cells into four phase assignments based on HDHB-mClover expression and their number of nuclei (Fig. 5A). Mononucleated cells were assigned as G, S/G2 or M phase and multinucleated cells assigned to either M or multinucleated phase based on localization of the HDHB-mClover signal (Supplemental Fig. 5A).

Fig. 5
figure 5

ELF3 mediates cell cycle progression under paclitaxel treatment. (A) Representative images showing the HCC1143 cell cycle reporter line and a mitotic even occurring over 105 min. Orange text indicates assigned cell cycle. (B) Relative cell count for each cell cycle phase over time for each siRNA condition + /− paclitaxel (PTX). Cell counts normalized to total cell number at earliest time point. (C) Schematic showing the underlying structure of permitted transitions used in the Markov Model. (D) Barplot showing inferred cell cycle phase duration from Markov model transition rates for each treatment combination. (E) PTX + siRNA synergy computed as the ratio of inferred phase duration for combination (siRNA + PTX) versus Highest Single Agent (HSA, highest duration for either siRNA or PTX treatment alone). Value of 1 indicates no change in combination, values greater than 1 indicate synergy and values less than 1 indicate antagonism. (F) MSigDB Gene Set Enrichment (GSEA) results for ELF3 high versus ELF3 groups for samples from the METABRIC cohort. (G) Overall survival for the METABRIC breast cancer cohort stratified by ELF3 mRNA expression. High = top quartile of ELF3 expression, IQR = inner quartile range of ELF3 expression, and low = lowest quartile of ELF3 expression.

Full size image

As chemotherapeutic drugs often have peak efficacy during a specific cell cycle phase, we next asked whether the combination of paclitaxel treatment and siRNA knockdown altered the dynamics of cell cycle progression. To that end, we trained a Markov Model on the live-cell data, which enabled us to infer transition rates and the average time spent in each of the four phases for a given treatment condition44,48. This approach uses the change in fraction of cells in each cell cycle phase over time (Fig. 5B) to learn cell cycle-specific transition rates, which represent the fraction of cells that transition from one phase to another phase within a 1 hour timestep (Supplemental Fig. 5B). We constrained our model such that proliferating cells can either successfully complete the cell cycle or undergo mitotic failure into a permanent multinucleated phase (Fig. 5C). The resultant model can reconstruct the cell phase distribution and counts for the entire experimental duration using just the data from the initial timepoint and learned transition rates (Supplemental Figs. 5C,D). The model-predicted values for both siRNA and siRNA + Paclitaxel conditions had low error (Supplemental Fig. 6A,B), and the majority of the model predictions (11 of 12) had over 95% agreement with observed counts across all timepoints (Supplemental Fig. 6C).

We leveraged this model to assess how the combination of siRNA knockdown with paclitaxel synergized to disrupt the cell cycle and modulate transitions between phases. For each individual or combination perturbation, we computed the inferred phase duration for G1, S/G2 and M phases using the model’s homotypic transition rates, which represent the fraction of cells that remain in the same phase through the timestep (Fig. 5D). The inhibition of ELF3 alone strongly increased cell cycle duration (DMSO + siNonTarget = 49 h, DMSO + siELF3 = 115 h), with substantial increases to the time spent in G1 (DMSO + siNonTarget = 7.9 h, DMSO + siELF3 = 33.5 h) and S/G2 (DMSO + siNonTarget = 37.7 h, DMSO + siELF3 = 72.1 h) phases. The combination of paclitaxel treatment and siRNA knockdown resulted in the longest cell cycle durations for both the PTX + siELF3 (160 h) and PTX + siIRF9 (159 h) conditions compared to PTX + siNonTarget (84.4 h). To further assess the therapeutic impact of siRNA knockdown and paclitaxel combination treatments, we used a Highest Single Agent (HSA) model to compare the inferred phase duration under combination paclitaxel + siRNA to the highest inferred phase duration for a single agent (either paclitaxel alone or siRNA alone, Fig. 5E)49. In this approach, a duration ratio less than 1 indicates antagonism (siRNA + paclitaxel results in shorter inferred phase duration than highest single agent), a duration ratio of 1 means there is no benefit of combination compared to highest single agent, and a duration ratio greater than 1 indicates a positive synergistic effect (siRNA + paclitaxel results in longer inferred phase duration than highest single agent). Although knockdown of NFE2L2 alone resulted in longer cell cycle phases, these changes were not particularly synergistic with paclitaxel treatment. In contrast, IRF9 knockdown resulted in increased duration ratios of all three phases compared to HSA (G1: 1.93, S/G2: 1.87, M: 1.8), while FOSL1 knockdown resulted in increased duration ratios for G1 and M phases (1.96, 2.19 ratios respectively). ELF3 knockdown showed the greatest synergy for the G1 phase, with a G1 duration ratio of 2.55, indicating that the combination of ELF3 knockdown with paclitaxel treatment strongly inhibits cell cycle progression out of G1.

Motivated by these ELF3 observations, we hypothesized that ELF3 expression may be predictive of overall survival in breast cancer. To that end we assessed the METABRIC50 breast cancer cohort, using ELF3 expression to stratify patients into three categories: ‘high’ (Upper quartile of ELF3 expression), ‘mid’ (Inter quartile range of ELF3 expression), and ‘low’ (Lower quartile of ELF3 expression). We then compared the gene expression between the ELF3-high and ELF3-low groups and found that the ELF3-high tumors were significantly enriched for MSigDB hallmarks related to cell cycle progression (G2M_CHECKPOINT: NES = 1.99, Benjamini–Hochberg FDR = 5.0e−8, E2F_TARGETS: NES = 1.94, Benjamini–Hochberg FDR = 2.9e−7, Fig. 5F, Supplemental Fig. 7A). ELF3-low tumors were enriched for MSigDB hallmarks related to Allograft Rejection and Epithelial to Mesenchymal Transition (NES = -2.14, Benjamini–Hochberg FDR = 4.1e−11, and NES = − 1.95, Benjamini–Hochberg FDR = 6.1e−8 respectively). We next asked whether ELF3 expression levels were related to recurrence free survival in breast cancer patients, with the hypothesis that patients with high ELF3 expression would have worse survival. This comparison indicated that ELF3 high expressing tumors had significantly worse outcome than those with low ELF3 expression (HR = 1.57, p value = 1e−5, Fig. 5G). Additionally, this analysis found that ELF3 expression was prognostic in both directions when compared to the ELF3 mid-expressing group: ELF3-high tumors have significantly shorter recurrence free survival (HR = 1.21, p value = 0.025, Cox Proportional Hazard) and ELF3-low tumors have significantly longer overall survival (HR = 0.77, p value = 0.004, Cox Proportional Hazard) compared to the ELF3-mid tumors (Fig. 5G). These results support our experimental in vitro findings that ELF3 activity contributes to malignant cell proliferation, that high ELF3 expression is associated with cell cycle progression in human breast tumors, and finally that ELF3 may serve as a biomarker of progression free survival.

Discussion

Paclitaxel is a cornerstone therapy for TNBC and is an important component of first line neoadjuvant treatment for newly detected disease. Despite this, less than 20% of breast cancer patients treated with combination neoadjuvant therapy achieve pathological complete response (pCR), and 47% of TNBC patients without pCR have recurrent disease within 10 years51. Although long-term chemotherapy resistance is often facilitated by clonal selection of growth-permissive mutations52,53,54, newer molecular profiling techniques have revealed that short-term adaptive responses emerge through rapid epigenetic changes without acquisition of new mutations4,55. In this study, we sought to identify adaptive responses that emerge after paclitaxel treatment and that may be targeted to deepen therapeutic response. To that end, we characterized the phenotypic and transcriptional responses of TNBC cells to paclitaxel, with a focus on changes in cell number, multinucleation, and transcription factor programs. Using siRNA knockdown, live-cell imaging, and computational modeling, we identified several TFs that phenocopied key aspects of paclitaxel response, including reduced proliferation rates and an increased proportion of multinucleated cells. ELF3 knockdown in vitro was synergistic with paclitaxel treatment and suppressed G1 to S/G2 cell cycle progression. Analysis of the METABRIC breast cancer cohort revealed that high expression of ELF3 was associated with worse outcome and higher cell-cycle related pathway activity. Together, these findings reveal a new role for ELF3 as a novel therapeutic target and biomarker of progression free survival for patients with TNBC.

Many prior drug and gene manipulation studies have focused on viability or other cell count proxies at a terminal timepoint56,57,58,59. While valuable, more recent studies have demonstrated that therapies modulate multiple cancer-associated hallmarks, including cell cycle phase behavior, senescence and nuclear morphology44,60,61. Further, cellular systems are inherently dynamic, and measures that capture temporal behavior are critical for gaining mechanistic insights45,62,63,64. While our live-cell studies captured important changes in cell cycle dynamics and the population distribution of various cell cycle states, no single metric captured the complete biological response. Future studies could deploy a richer panel of reporter molecules to gain deeper insights into the timing and order of transcription factor activation, activation of specific cell cycle checkpoints, or assess critical breast cancer pathways such as ERK, senescence or apoptosis6566,67,68.

In this study we identified dual roles of the transcription factor ELF3 that contribute to paclitaxel tolerance by: 1) permitting cells to transition from G1 to S/G2, and 2) enabling successful division into two mononuclear daughter cells. These findings were enabled by a Markov Model of cell cycle progression built on population level cell count data to learn the transition rates between cell cycle phases and inferred cell cycle phase durations44,48. While the inferred cell cycle durations represent an accurate prediction of the population’s average behavior, they cannot inform whether this arises from a homogenous or heterogenous distribution of cell cycle durations. This is of particular interest in the case of cancer treatment, as a small population of cells with a fitness advantage may eventually overtake other populations to achieve therapeutic resistance69. An alternative approach could track individual cells and their progeny to build complete lineages with accompanying cell cycle timing information. Lineage based approaches tend to be relatively low throughput due to the computational and experimental requirements, but they offer the opportunity to discern between heterogenous states of differing cycling speeds48. Another limitation of our Markov Model’s implementation is the assumption that transition rates are static throughout the observation duration. While the output of the model mapped well within the 72 h measurement window, there was some divergence at the end of the experiment that may suggest a weakening of either siRNA or paclitaxel effect. An updated model with temporally dynamic transition rates could provide additionally accuracy, enabling future studies aimed at prediction of drug combination effects and optimizing the drug schedule for maximum disruption of cell cycle progression. Optimal drug scheduling would be particularly interesting to examine, as TNBC disease is commonly treated with combination of paclitaxel with additional chemotherapies3,7,51.

Paclitaxel inhibits cell growth by simultaneously promoting microtubule assembly and inhibiting microtubule depolymerization, which results in mitotic checkpoint failure and subsequent apoptosis or senescent arrest19,70. The in vitro experimentation performed in this study represents an extensive investigation into the phenotypic and molecular responses of TNBC cells to paclitaxel, however we acknowledge that tumors are comprised of diverse cell types and microenvironmental signals not captured in this study71. Importantly, the tumor microenvironment has been shown to influence therapeutic response, including through cell–cell interactions72,73,74 and extracellular matrix remodeling75,76. In line with this, we observed that paclitaxel induced changes to malignant epithelial cells that suggest an altered capability to interact with immune cells via secretion of CXCL1 and CXCL8. Tumor-derived CXCL1 can recruit immunosuppressive myeloid cells that in turn inhibit CD8+ T cell infiltration77 while CXCL8 has been shown to recruit immunosuppressive neutrophils78, promote angiogenesis79 and maintain breast cancer stem cells80. Future studies that more deeply consider the influence of stromal and immune cell signals in modulating therapeutic response will be needed to create a complete picture of paclitaxel resistance.

As key regulators of multiple molecular programs, many transcription factors contribute to cancer-associated phenotypes81 and therapeutic response82,83. We found that the ETS family transcription factor ELF3 was upregulated during early response to paclitaxel treatment, and siRNA knockdown of ELF3 synergized with paclitaxel treatment to slow cell growth. Other studies have found that high ELF3 activity is associated with inhibition of epithelial to mesenchymal transition 84. Furthermore, ELF3 inhibition reduces proliferation across multiple cancer models including lung adenocarcinoma85, neuroendocrine carcinoma86 and prostate cancer87. Conserved dysregulation of ELF3 across cancer types may be related to its genomic ___location (loci 1q32) which is commonly amplified across cancers89,90 and also encodes for a number the cancer related genes including MDM4 (p53 suppressor)91,92. Our results support the idea that ELF3 is not only commonly dysregulated across these epithelial malignancies, but can also act as a transcriptional mechanism of acquired drug resistance. Future works that employ similar methodology to probe the function of ELF3 in these other malignancies, as well as systematically perturb combinations of the implicated transcription factors, may provide additional insight into this complex response phenotype and how it can be interrupted for therapeutic benefit.

Taken together, this work has identified ELF3 upregulation as an acquired mechanism of paclitaxel resistance. These findings support the development of pharmacological agents that inhibit ELF3 activity and could be used in combination with paclitaxel to further improve patient outcomes. While it has been historically difficult to develop targeted transcription factor inhibitors due to their lack of enzymatic activity, recent advances, such as targeted siRNA nanoparticles and indirect inhibition through targeting multiple interacting proteins, have made pharmacomodulation of transcription factors more tenable93,94,95. Until such therapies are developed, ELF3 may serve as a useful biomarker which predicts the development of paclitaxel resistance and continued malignant proliferation.

Methods

Cell culture

HCC1143 (ATCC), HCC1806 and MDA-MB-468 cells were authenticated by STR profiling and tested negative for mycoplasma. HCC1143 and HCC1806 cells were cultured in RPMI 1640 with L-glutamine (cat. 11875119, Life Technologies Inc.) supplemented with 10% fetal bovine serum (#16000-044, Gibco). MDA-MB-468 cells were cultured in DMEM (#11965-092, Life Technologies Inc.) supplemented with 10% fetal bovine serum (#16000-044, Gibco). All lines were incubated at 37C with 5% CO2. For perturbation experiments, cells were seeded into appropriate assay vessel for 24 h prior to treatment with either vehicle control (DMSO; PBS) or perturbation (table below).

Perturbation

Shorthand

Concentration used for scRNA-seq

Source

Identifier

Vehicle

Paclitaxel

PTX

1 nM

LC labs

P-9600

0.1% DMSO

Notch inhibitor

NOTCHi

1 uM

Millipore sigma

BM0018-5MG

0.1% PBS

Interferon beta

IFNB

20 ng/mL

PBL assay science

11,410–2

0.1% PBS

Interferon gamma

IFNG

20 ng/mL

R&D systems

385-IR-100

0.1% PBS

Transforming growth factor beta

TGFB

10 ng/mL

Biotechne

7754BH005

0.1% PBS

Lymphotoxin

LT

10 ng/mL

Biotechne

8884-LY-025

0.1% PBS

Oncostatin M

OSM

10 ng/mL

Cell signaling technology

5367SC

0.1% PBS

Fixed cell assays

Cells were plated at 3000 cells in 100ul of complete media per well in a 96 well plate (#08-772-225, FisherScientific). After 24 h, an additional 100ul of either vehicle (0.1% DMSO) or paclitaxel containing complete media was added. After 72 h cells were fixed with 4% Formaldehyde (#28908, ThermoFisher Scientific) for 15 min at room temperature, then permeabilized with 0.3% Triton X-100 (#X100-100ML, Sigma Aldrich) for 10 min at room temperature, then washed twice with PBS. Fixed cells were blocked with 1% BSA (A7906-100G, Millipore Sigma) in PBS for 1 h at room temperature and then stained overnight with 1:100 anti-CDKN2A/p16INK4A + CDKN2B/p15INK4B-AF644 (#ab199756, Abcam), and 1:100 anti-cPARP-AF647 (#6987S, Cell Signaling Technology) or 1:500 anti-TUBB3-AF647 (#ab190575, Abcam) overnight at 4C. Each well was washed twice with room temp PBS then stained with 0.5ug/mL DAPI (4083S, Cell Signaling Technology) in PBS for 15 min at room temperature. Following DAPI staining, wells were washed once with PBS, then stained with 1:20,000 HCS CellMask (Orange: #H32713, Green: #H32714, Invitrogen. Used for cytoplasmic segmentation) in PBS for 15 min at room temperature. Wells were washed twice with room temperature PBS and then 4 fields of view per well imaged on an InCell 6000 (GE Healthcare). Images were segmented with two custom Cellpose47 models to segment the nucleus (from DAPI channel) and cytoplasm (from HCS Cellmask channel). Image quantification was performed in R (v4.3.1) using EBImage (v4.42.0), and cells were annotated based on the number of distinct nuclei segmented within each cytoplasmic mask. EdU incorporation assays (#C10640, Invitrogen) were performed following manufacturer instruction with default concentrations, with an initial incorporation incubation time of 1 h at 37C, and additional staining with DAPI (#4083S, Cell Signaling Technology) and HCS CellMask (#H32714, Invitrogen) to aid in segmentation. Segmentation was performed using custom trained Cellpose models and feature quantification was performed with EBImage (v4.42.0).

scRNA-seq library preparation and sequencing

Experiment 1 (DMSO 24 h, DMSO 72 h, Paclitaxel 24 h, Paclitaxel 72 h): Each condition had a single-cell RNA library prepared using 10X Genomics Single Cell 3’ v2 kits and sequenced on an Illumina NextSeq 500 for 500e6 reads per library.

Experiment 2: All conditions were multiplexed using Hashtag Oligonucleotide barcoding technology (TotalSeq-B, Biolegend) following manufacturer standard protocol. A paired feature-barcode library and mRNA library were generated using the Single Cell 3’ v3 kit (10X Genomics) following manufacturer instructions and then sequenced on an Illumina NovaSeq for 800e6 reads.

scRNA-seq data processing

For both experiments; raw base call files were converted to FASTQ format with bcl2fastq (Illumina). Cellranger count (v6.0.2) was used to align reads to the GRCh38 transcriptome (GRCh38-2020-A, accessed from 10X Genomics) and count UMI reads. The R package Seurat96,97 (4.0.5) was used to perform variable feature identification, linear and nonlinear dimensionality reduction, unsupervised clustering and differential gene expression.

Variance Stabilizing Transformation was used to identify the top 2000 variable genes and Principal Component Analysis (PCA) was used to reduce these 2000 genes to 10 components for UMAP embedding and unsupervised clustering. Differential expression analysis was performed using the FindMarkers function of Seurat with default parameters. Geneset enrichment analysis was performed with the R package clusterProfiler98 (v4.8.2) using significantly upregulated genes compared to time-matched vehicle control (abs(avg_log2FC) > 0.5, Benjamini Hochberg FDR < 0.05).

Transcription factor enrichment analysis

Significantly upregulated genes (avg_log2FC > 0.5, Benjamini Hochberg FDR < 0.05) were computed for paclitaxel, IFNB and IFNG treated samples compared to time-matched vehicle treated cells. ChEA3 enrichment analysis was performed with default settings using R code from the CHEA3 API documentation (https://maayanlab.cloud/chea3/) to perform an online query using either the genes uniquely upregulated in paclitaxel treated cells, or those shared between paclitaxel and either of the interferon responses. The top 15 ranked transcription factors from both the paclitaxel unique and paclitaxel-interferon shared TF enrichment lists were considered when nominating siRNA knockdown targets. Any TF that also had at least 0.25 log2 fold change for paclitaxel at either 24 or 72 h compared to vehicle control was included in the siRNA knockdown panel.

siRNA knockdown

Cells were plated in 90ul of serum free media per well of a 96 well plate. 24 h later, siRNA knockdown mixture was prepared using a cell-line optimized concentration of Lipofectamine RNAiMAX (cat 13,778,075–075, Invitrogen) and siRNA (Horizon Discovery ON-TARGETplus) following RNAiMAX recommended protocol. The final concentration of siRNA per well was 1 pmol and the final volume of RNAiMAX per well was 75nL for HCC1143, and 37.5nL for HCC1806 or MDA-MB-468 in 100uL of cell containing volume. 24 h after siRNA transfection cells were treated with an addition of 100uL complete media containing either DMSO vehicle control or paclitaxel. The Z’-factor was positive for each cell-line by drug combination, indicating a valid screening assay with sufficient dynamic range between positive and negative controls (Supplemental Data 5)99.

Protein isolation

Protein isolation: siRNA knockdown of HCC1143 cells was performed using the siNonTarget, siELF3, siFOSL1, and siNFE2L2 pools as described above. After 24 h of knockdown, perturbation containing media was added such that media volume doubled and had a final concentration of either 0.1% DMSO (vehicle control) or 1 nM Paclitaxel. After 72 h of perturbation, cells were washed with 4C PBS then lysed by 5 min incubation at 4C with RIPA buffer (R0278, Sigma) supplemented with 1X Halt Protease and Phosphatase Inhibitor Cocktail (1861281, Thermo Scientific). Remaining cells were scraped from the plate and lysate was snap frozen in liquid nitrogen then stored at − 80C overnight. The following day lysate was clarified by centrifugation at ×21,130g for 10 min at 4C. The supernatant was collected and the protein concentration was immediately quantified. Remaining protein was stored at − 80C.

Western blot

Protein quantification was performed using the Western Simple protocol on the Jess capillary western machine using the 12–230 kDa cartridge and following manufacturer instructions (SM-W002, Biotechne). Primary antibodies targeting the protein products of ELF3 (anti-ESE1, NBP3-32320, Novus), FOSL1 (anti-FRA1, #5281S, Cell Signaling Technology), and NFE2L2 (anti-NRF2, #20733S, Cell Signaling Technology) were used at 1:50 dilution. Lysates were loaded at a concentration of 2 mg/mL and volume of 5uL per capillary well, and the Anti-rabbit detection kit (DM-001, Biotechne). was used to quantify primary antibody levels. The Total Protein Detection Module (DM-TP01) was used to normalize input protein across all lanes. Quantification was performed using the included Compass software with default settings (v6.3.0, Biotechne), and the resultant peak areas was plotted using R (v4.3.1).

HDHB reporter live-cell assays

siRNA knockdown and drug treatment was performed as described above, and then the plate was loaded on an Incucyte S3 (Sartorious) and cells imaged every 15 min for 72 h post drug treatment. At each timepoint 4 fields of view were captured at 20 × magnification in each well using the phase, red and green channels. A cytoplasmic mask was computed from the mean of normalized red/green channel, and a nuclear mask was computed from the red channel using custom trained Cellpose47 models. Image quantification was performed in R (v4.3.1) using EBImage (v4.42.0). An additional perinuclear ring mask was computed as the 11 pixel dilation from the nuclear mask, but still bound by the cytoplasmic mask. To determine mClover localization thresholds for cell cycle assignment, 250 cell images were randomly selected and manually assigned to the G1, S/G2 or M cell cycle state based on mClover localization. The mClover intensity ratios were then used to determine thresholds for automated cell cycle phase calling which was applied to the rest of the data set (Supplemental Fig. 5A). Mononuclear cells with a Perinuclear:Nuclear mean intensity ratio greater than 0.8 and Nuclear:Cytoplasmic total intensity less than 0.5 were assigned to the S/G2 phase. Mononuclear and Multinuclear cells with a Nuclear:Cytoplasmic total intensity ratio greater than 0.8 and Perinuclear:Nuclear mean intensity ratio less than 0.8 were assigned to the ‘M’ phase. The remainder of mononuclear cells were assigned ‘G1’, and the remainder of multinucleated cells were assigned ‘Multinucleated’.

Markov modeling

The 5-frame moving average of cell count per cell cycle phase was down sampled to one value per hour and used to train a Markov model for each unique siRNA (NonTarget, ELF3, FOSL1, NFE2L2, IRF9, PLK1) +/− paclitaxel condition. The transition matrix of the model was constrained such that cells could remain in their current phase, progress through the cell cycle (G1—> S/G2, S/G2—> M, M—> G1 with replication) or transition from M phase to an absorbing (permanent) multinucleated phase. Models were trained for 15 epochs, and the first epoch was seeded with an identity transition matrix. 3000 random transition matrices were generated each epoch, and the 5 with lowest error were used as seeds for the following epoch. The prior best performing matrices were updated with randomly generated matrices at a learning rate of 0.1 for the first epoch, halving every 2 epochs.

The prediction for counts for each future state (Sn+1) is calculated as the product of the counts at the prior state (Sn) by the transition matrix (P) and the replication matrix (RM).

$$RM= \begin{array}{ccccc} & G1& S/G2& M& Multi.\\ G1& 1& 1& 0& 0\\ S/G2& 0& 1& 1& 0\\ M& 2& 0& 1& 1\\ Multi.& 0& 0& 0& 1\end{array}$$
$$P = \begin{array}{ccccc} & G1& S/G2& M& Multi.\\ G1& ?& ?& 0& 0\\ S/G2& 0& ?& ?& 0\\ M& ?& 0& ?& ?\\ Multi.& 0& 0& 0& ?\end{array}$$
$${S}_{n+1}={S}_{n}*P*RM$$

The error of the Markov predicted cell counts (cexp) compared to observed counts (cobs) was computed as the arithmetic mean of the Root Mean Squared Relative Error (RMSRE) of each cell cycle phase across all predicted timepoints. The noise floor of RMSRE was estimated with a second-order loess fit with span of 0.75 (loess function from R package ‘stats’, v4.3.1).

$$RMSRE= \sqrt{\frac{1}{n}*\sum \frac{{({c}_{exp}-{c}_{obs})}^{2}}{{{c}_{obs}}^{2}}}$$

The mitotic success rate (MSR) of each condition was computed as the ratio of M-to-G1 transition (PM,G1) to the sum of the transition rates for M-to-G1 (PM,G1) and M-to-multinucleated (PM,Multi):

$$MSR=\frac{{P}_{M,G1}}{{P}_{M,G1}+{P}_{M,Multi}}$$

The expected duration of G1, S/G2 and M cell cycle phases was calculated from the homotypic transition rates as100:

$$For (i==j): {ExpectedDuration}_{i,j}= \frac{1}{1-{P}_{i,j}}$$

METABRIC survival and microarray analysis

The METABRIC50 microarray and patient metadata was accessed through cbioportal101,102,103 and analyzed using R (v4.3.2) and the ‘survival’ package (v3.5.7). The z-scored microarray expression data was used to categorize patients into ‘high’ (highest expressing quartile), ‘mid’ (first to third expressing quartile) or ‘low’ (lowest expressing quartile) based on expression of ELF3. For survival analysis, patients were filtered to those with microarray data and then Kaplan–meier survival curves were generated with the ‘ggsurvfit’ package (v1.0.0). Cox proportional hazard statistics were calculated with the ‘coxph’ function of the ‘survival’ package (v3.5.7). Differential expression was calculated from the log normalized microarray data using the ‘wilcoxauc’ function from the ‘presto’ package (v1.0.0). Significantly differentially expressed genes (abs(logFC) > 0.5 and adjusted p < 0.05) where used to compute MSigDB hallmark GSEA using the ‘clusterprofiler’ (v4.10.1) and ‘msigdbr’ (v7.5.1) packages.

Full reagent list

Reagent

Shorthand

Type

Source

Identifier

HCC1143

–

Cell line

ATCC

CRL-2321

HCC1806

–

Cell line

ATCC

CRL-2335

MDA-MB-468

–

Cell line

ATCC

HTB-132

RPMI 1640

RPMI

Reagent

Life Technologies

11875119

DMEM

DMEM

Reagent

Life Technologies

11965-092

Fetal Bovine Serum

FBS

Reagent

Gibco

16000-044

Dimethyl Sulfoxide

DMSO

Reagent

Millipore Sigma

D8418-250ML

Phosphate Buffered Saline

PBS

Reagent

Gibco

14190235

Paclitaxel

PTX

Reagent

LC Labs

P-9600

Interferon-Beta

IFNB

Reagent

PBL Assay Science

11410-2

Interferon-Gamma

IFNG

Reagent

R&D Systems

385-IR-100

Human Oncostatin M

OSM

reagent

Cell Signaling Technology

5367SC

Recombinant Human Lymphotoxin alpha1/beta2 protein

LT

Reagent

Biotechne

8884-LY-025

Recombinant Human TGFB-Beta 1

TGFB

Reagent

Biotechne

7754BH005

BMS-906024

NOTCHi

Reagent

Millipore Sigma

BM0018-5MG

16% Formaldehyde (w/v)

–

Reagent

ThermoFisher Scientific

28,908

Triton X-100

–

Reagent

Millipore Sigma

X100-100ML

Normal Goat Serum Blocking Solution

–

Reagent

MP Biomedicals

#0219135680

Lipofectamine RNAiMAX transfection Reagent

–

Reagent

ThermoFisher Scientific

13778075

Bovine Serum Albumin

BSA

Reagent

Millipore Sigma

A7906-100G

RIPA buffer

–

Reagent

Sigma

P0278

100X Halt Protease and Phosphatase Inhibitor Cocktail

–

Reagent

ThermoFisher Scientific

1861281

Click-iT Plus EdU Cell Proliferation kit (AF647)

–

Reagent

Invitrogen

C10640

anti-CDKN2A/p16INK4A + CDKN2B/p15INK4B-AF644

p16

Antibody

Abcam

ab199756

anti-cPARP-AF647

cPARP

antibody

Cell Signaling Technology

6987S

anti-TUBB3-AF647

TUBB3

Antibody

Abcam

ab190575

HCS CellMask Green

CellMask Green

Stain

Invitrogen

H32713

HCS CellMask Orange

CellMask Orange

Stain

Invitrogen

H32714

DAPI

DAPI

Stain

Cell Signaling Technology

4083S

ELF3/ESE-1 Antibody

–

Antibody

Novus

NBP3-32320

FRA1 (D80BP) Rabbit mAb

–

Antibody

Cell Signaling Technology

5281

NRF2 (E5F1A) Rabbit mAb

–

Antibody

Cell Signaling Technology

20733

12–230 kDa separation module

–

Reagent

Biotechne

SM-W002

Total Protein Detection Module

–

Reagent

Biotechne

DM-TP01

Anti-Rabbit Detection Module

–

Reagent

Biotechne

DM-001

siATF3 Smartpool

siATF3

siRNA

Hoizon Discovery

L-008663-00

siDDIT3 Smartpool

siDDIT3

siRNA

Hoizon Discovery

L-004819-00

siELF3 Smartpool

siELF3

siRNA

Hoizon Discovery

L-016080-00

siFOSL1 Smartpool

siFOSL1

siRNA

Hoizon Discovery

L-004341-00

siIRF7 Smartpool

siIRF7

siRNA

Hoizon Discovery

L-011810-00

siIRF9 Smartpool

siIRF9

siRNA

Hoizon Discovery

L-020858-00

siJUNB Smartpool

siJUNB

siRNA

Hoizon Discovery

L-003269-00

siJUN Smartpool

siJUN

siRNA

Hoizon Discovery

L-003268-00

siKIF11 Smartpool

siKIF11

siRNA

Hoizon Discovery

L-003317-00

siKLF6 Smartpool

siKLF6

siRNA

Hoizon Discovery

L-021441-00

siMAFF Smartpool

siMAFF

siRNA

Hoizon Discovery

L-003903-00

siNFE2L2 Smartpool

siNFE2L2

siRNA

Hoizon Discovery

L-003755-00

siPLK1 Smartpool

siPLK1

siRNA

Hoizon Discovery

L-003290-00

siPLSCR1 Smartpool

siPLSCR1

siRNA

Hoizon Discovery

L-003729-00

siSP100 Smartpool

siSP100

siRNA

Hoizon Discovery

L-015307-00

ON-TARGETplus Non-targeting Control

siNonTarget

siRNA

Hoizon Discovery

D-001810-10

ON-TARGETplus GAPD Control

siGAPD

siRNA

Hoizon Discovery

D-001830-10