Coordination under stress, for better or worse
Cells adapt to oxidative stress in part by expressing antioxidant genes, many of which are transcribed by NRF2 and p53. Using three-dimensional cultures of normal and premalignant breast epithelial cells and mathematical modeling, Pereira et al. found that NRF2 and p53 coordinated the oxidative stress response through cooperation and mutual compensation. Although critical for normal duct development, their models suggested that NRF2-mediated tolerance to oxidative stress in premalignant tissue may permit the emergence of p53 mutations that drive malignant progression and render NRF2 dispensable. The findings reveal further complexity and cell state specificity in redox signaling networks and the relevance of these networks to normal tissue development, tumor progression, and therapeutic strategies.
Abstract
Breast and mammary epithelial cells experience different local environments during tissue development and tumorigenesis. Microenvironmental heterogeneity gives rise to distinct cell regulatory states whose identity and importance are just beginning to be appreciated. Cellular states diversify when clonal three-dimensional (3D) spheroids are cultured in basement membrane, and one such state is associated with stress tolerance and poor response to anticancer therapeutics. Here, we found that this state was jointly coordinated by the NRF2 and p53 pathways, which were costabilized by spontaneous oxidative stress within 3D cultures. Inhibition of NRF2 or p53 individually disrupted some of the transcripts defining the regulatory state but did not yield a notable phenotype in nontransformed breast epithelial cells. In contrast, combined perturbation prevented 3D growth in an oxidative stress–dependent manner. By integrating systems models of NRF2 and p53 signaling in a single oxidative stress network, we recapitulated these observations and made predictions about oxidative stress profiles during 3D growth. NRF2 and p53 signaling were similarly coordinated in normal breast epithelial tissue and hormone-negative ductal carcinoma in situ lesions but were uncoupled in triple-negative breast cancer (TNBC), a subtype in which p53 is usually mutated. Using the integrated model, we correlated the extent of this uncoupling in TNBC cell lines with the importance of NRF2 in the 3D growth of these cell lines and their predicted handling of oxidative stress. Our results point to an oxidative stress tolerance network that is important for single cells during glandular development and the early stages of breast cancer.
INTRODUCTION
Among glandular tissues, the breast-mammary epithelium is unique because of the marked expansion and reorganization that occur after birth (1). During puberty, a branched network of epithelial ducts is pioneered by terminal end buds (TEBs), which emerge from the rudimentary gland and extend into the surrounding mesenchyme (2). TEBs contain a mixture of proliferating stem-progenitor cells and differentiating cells fated to the secretory luminal-epithelial or contractile basal-myoepithelial lineages. During morphogenesis, TEB cells are dynamically exposed to different microenvironments that inform final organization of the gland (3). Some microenvironmental cues are supportive or instructive to cells [hormones (4), growth factors (5), and basement membrane (6)]. Others are deleterious or lethal [loss of polarity (7), detachment (8), and endoplasmic reticulum stress (9)]. All of these cues are reconfigured aberrantly and heterogeneously during the early stages of breast-mammary cancer (10–12).
Stress and survival signals also juxtapose when breast-mammary epithelial cells are grown in three-dimensional (3D) culture with reconstituted basement membrane extracellular matrix (ECM) (13, 14). Combining the appropriate adhesive and soluble cues yields TEB-like behavior in 3D-cultured multicellular epithelial fragments from the mammary gland (7). For single-cell cultures that reliably organize as 3D structures, clones or progenitors must iteratively proliferate, maintain cell-cell adhesions, and coordinate function to establish a multicellular ecosystem (15, 16). Cell regulatory states diversify within 3D organoids of primary breast-mammary epithelia (17–19) and also in the simplest 3D spheroids of isogenic cell lines (20–23). Identifying such cell regulatory heterogeneities is important, because there are parallels to in situ lesions of the breast, where premalignant cells must survive and proliferate in the duct (24, 25).
Previously, we identified a cluster of transcripts (Fig. 1A, top) that covaries heterogeneously among hormone-negative, basal breast epithelial cells grown as 3D spheroids (20). The cluster contains KRT5 (a PAM50 classifier for basal-like breast cancer) (26) along with multiple stress tolerance genes, including JUND (27), CDKN1A (28), MUS81 (29), and HSPE1 (30). The transcripts in this cluster were among the strongest and most-negative predictors of breast cancer response to chemotherapy and targeted agents in an independent clinical trial (31). We reported that individual genes in the cluster have complex time- and microenvironment-dependent relationships in 3D spheroids, animal models of ductal carcinoma in situ (DCIS), and clinical hormone-negative premalignancies (24). However, the overarching regulation of the cluster was not determined.
(A) Maximum likelihood inference parameterization (bottom) of a two-state distribution of transcript abundances for the gene cluster from microarray profiles (top) of ECM-attached basal-like MCF10A-5E breast epithelial cells, randomly collected as 10-cell pools (n = 16) from 3D-cultured spheroids after 10 days, extracted from (20). Inferred expression frequencies are the maximum likelihood estimate with 90% confidence interval (CI). (B) Venn diagram summarizing the candidate TFs predicted from four different bioinformatics algorithms (data file S1). (C and D) Quantitative immunofluorescence of (C) hyperphosphorylated RB (pRB, an upstream proxy of active E2F1) and (D) NRF2 in 3D culture with ECM (top), 2D culture (middle), and 2D culture with ECM (bottom). Expression frequencies for a two-state lognormal mixture model (preferred over a one-state model by F test; P < 0.05) were calculated by nonlinear least squares of 60 histogram bins collected from n = 1100 to 1600 of cells quantified from 100 to 200 spheroids from two separate 3D cultures. For each subpanel, representative pseudocolored images are shown in the top right inset and merged (magenta) with DAPI nuclear counterstain (blue) in the bottom right inset. Scale bars, 10 μm.
Here, we found that regulatory state heterogeneity emerges from the coordinated action of two stress-responsive transcription factors (TFs)—NFE2L2 (NRF2) (32, 33) and TP53 (p53) (34)—which become stabilized posttranslationally when breast epithelial cells variably experience oxidative stress in 3D culture. Genetic disruption of NRF2 signaling altered the transcriptional cluster, but 3D phenotypes were buffered or redirected by compensatory increases in p53 signaling. Disabling p53 function synergized with NRF2 deficiency, suppressing normal 3D proliferation and promoting irregular hyperproliferation in a transformed-yet-premalignant derivative. These observations were consistent with an integrated systems model of NRF2-p53 signaling that encoded a shared oxidative stress trigger and common pool of antioxidant target genes without any further cross-talk. Among clinical specimens, NRF2-p53 coordination was retained in normal primary breast tissue and hormone-negative DCIS. However, the two pathways were largely uncoupled in triple-negative breast cancers (TNBCs), in which p53 is usually mutated (35). The integrated NRF2-p53 model predicted variable extents of uncoupling among TNBCs lines, and high uncoupling coincided with the most severe 3D growth alterations upon NRF2 knockdown. Past work on NRF2 in breast cancer has focused on its direct interactions with TNBC-associated tumor suppressors (36, 37). Our results suggest a broader systems-level role for NRF2 and p53 in oxidative stress tolerance of normal breast-mammary epithelia and hormone-negative premalignancies.
RESULTS
Statistical bioinformatics links gene cluster regulation to NRF2 and p53
We began by looking within the gene cluster (Fig. 1A, top) for potential regulatory mechanisms. The only TF in the cluster is JUND, and we showed previously that its chronic knockdown in MCF10A-5E cells (20) causes specific morphometric defects during spheroid growth (24). We revisited these results by acutely knocking down the expression of JUND with inducible short hairpin RNA (shRNA) and measuring transcript abundance of cluster genes by quantitative polymerase chain reaction (qPCR) (see Materials and Methods). Unexpectedly, other than JUND itself, no transcripts were reliably altered by its knockdown (see later in this section), supporting a regulatory role for other factors outside of the cluster.
We constrained the search for candidate regulators by using maximum likelihood inference (38) to estimate a frequency of bimodal transcriptional regulation (39) for the gene cluster. Given the 10-cell-averaged fluctuations from the original study (20), the maximum likelihood approach inferred two lognormal regulatory states defined by transcript abundance (Fig. 1A, bottom). The data supported a low-abundance regulatory state predominating in 58% of ECM-attached cells along with a second, high-abundance subpopulation in the remaining 42%. The frequency estimates placed quantitative bounds on the bimodal characteristics of upstream regulatory mechanisms.
Next, we applied a panel of bioinformatics approaches to search for TFs that might impinge upon the gene cluster (see Materials and Methods). The informatic methods adopt different strategies for assessing binding site overrepresentation (40–43). Therefore, we intersected their respective outputs to arrive at predictions that were robust to algorithmic details. The analysis converged upon two TFs: the G1/S regulator E2F1 and the stress response effector NRF2 (Fig. 1B and data file S1). We assessed the relative activation of the NRF2 and E2F1 pathways in single cells by quantitative immunofluorescence for the total stabilized NRF2 protein or phosphorylated RB1 (pRB indicates disinhibited E2F1; see Materials and Methods). In 3D spheroid cultures, pRB immunostaining was bimodal, but high-pRB cells were much rarer than the inferred regulatory frequency of the gene cluster (Fig. 1C, top). In 2D cultures, pRB staining was more than twice as immunoreactive and nearly twice as prevalent in the population (Fig. 1C, middle). The reduced proportion of high-pRB cells in 3D is consistent with the proliferative suppression of late-stage spheroid cultures (23). A 3D-like distribution of pRB was achieved in 2D cultures upon addition of dilute ECM (Fig. 1C, bottom), stemming from soluble proliferation-suppressing factors in the reconstituted basement membrane preparation (44). By contrast, NRF2 stabilization was only distinctly bimodal in 3D spheroids, and the observed frequency of low- and high-NRF2 states almost exactly coincided with that inferred for the gene cluster (Fig. 1D). Stabilization of hypoxia-inducible factor 1α (HIF-1α) was negligible in 3D spheroids overall (fig. S1, A and B), excluding irregular hypoxic stress as a contributor to the two-state distribution of NRF2. These results build a strong statistical argument for NRF2 as a covarying regulator of the gene cluster.
The NRF2-associated gene cluster (Fig. 1A, top) was originally identified by quantitative analysis of transcriptomic fluctuations among 4557 genes profiled by oligonucleotide microarray (20). The same samples were later reprofiled by 10-cell RNA sequencing (10cRNA-seq) (45), creating an opportunity to look more deeply at covariates with the NRF2-associated gene cluster. We used the median ranked fluctuations of the cluster across 10 cell samples (Fig. 1A, top) and surveyed the 10cRNA-seq data for genes that covaried (Spearman ρ > 0.5), identifying 633 candidates (Fig. 2A). When this expanded cluster was assessed for functional enrichments by Gene Ontology (GO) (data file S2) (46), we noted multiple GO terms linked to cell stress (“Response to stress” and “Oxidative stress”) and the TF p53 (“DNA damage response” and “p53 pathway”). p53 is sporadically stabilized in regenerating epithelia such as the intestine and skin, but p53 activation in quiescent tissues is rare (47). Recognizing the residual proliferation observed in 3D cultures (Fig. 1C), we immunostained for p53 and found nonuniform stabilization associated with the abundance of NRF2 in single cells [Fig. 2B, estimated mutual information (MI) = 0.15 (0.12 to 0.18); see Materials and Methods]. The analysis raised the possibility of a coordinated NRF2-p53 regulatory event triggered heterogeneously when breast epithelial cells proliferate and organize in reconstituted ECM.
(A) Transcripts covarying with the median NRF2-associated fluctuation signature (Fig. 1A, top) (20) measured by 10cRNA-seq (45) of ECM-attached MCF10A-5E cells grown as 3D spheroids (n = 18 10-cell pools from GSE120261). Selected GO enrichment analysis (green and purple) is shown for the transcripts with a Spearman correlation (ρ) greater than 0.5. The complete list of enrichments is available in data file S2. (B) Quantitative immunofluorescence of NRF2 and p53 abundance in ECM-attached MCF10A-5E cells grown as 3D spheroids. Representative pseudocolored images for NRF2 (top left) and p53 (middle left) are shown merged with DAPI nuclear counterstain (bottom left). White arrows indicate concurrent NRF2 and p53 stabilization. Median-scaled two-color average fluorescence intensities are quantified (right) along with the log-scaled and background-subtracted mutual information (MI) with 90% CI for n = 1691 cells segmented from 50 to 100 spheroids from two separate 3D cultures. (C) Genetic perturbation of NRF2 by inducible shRNA knockdown (top) and p53 by inducible expression of a FLAG-tagged carboxy terminal (residues 1 to 13 and 302 to 390) dominant-negative p53 (DNp53; bottom). NRF2 knockdown reduced NRF2 protein abundance to 22 ± 4% of control knockdown (fig. S3B). MCF10A-5E cells were treated with doxycycline (1 μg/ml) for 72 (top) or 24 (bottom) hours and immunoblotted for NRF2 or FLAG with vinculin, tubulin, and p38 used as loading controls and p21 used to confirm efficacy of DNp53. The negative control for shNRF2 was an inducible shGFP, and the negative control for DNp53 was an inducible FLAG-tagged LacZ. (D) Abundance changes in the gene cluster after single and combined perturbations of NRF2 and p53. NQO1 was used as a control for efficacy of shNRF2, and CDKN1A shows efficacy of DNp53. MCF10A-5E cells with or without NRF2 knockdown or DNp53 were treated with doxycycline (1 μg/ml) for 48 hours, grown as 3D spheroids for 10 days, and profiled for the indicated genes by qPCR. Data are log2 geometric mean relative to the negative control (shGFP + FLAG-tagged LacZ), with asterisks indicating statistically significant changes (left and middle columns) or interaction effects (right column) by two-way ANOVA of n = 8 independent 3D-cultured samples and an FDR of 5%. The complete set of transcripts in the gene cluster is shown in fig. S2C. (E) Dual inactivation of NRF2 and p53 causes synergistic proliferative suppression in MCF10A-5E 3D spheroids. Black arrows indicate proliferation-suppressed spheroids. Data are mean percentage of proliferation-suppressed spheroids ± SEM of n = 8 independent 3D-cultured samples after 10 days. Statistical interaction between NRF2 and p53 (Pint) was assessed by two-way ANOVA with replication. Scale bars, 20 μm (B) and 100 μm (E).
NRF2 coimmunoprecipitates with p53 in TNBC cells harboring gain-of-function p53 mutations, but this complex is absent in MCF10A cells with wild-type p53 (37). Loss of wild-type p53 function in MCF10A cells yields only minor 3D culture defects, but gain-of-function p53 mutants strongly perturb 3D architecture (48). Suspecting that some of p53’s effects could be explained through NRF2, we inducibly knocked down NRF2 with shRNA and inducibly coexpressed a truncated p53 (49) that acts as a dominant negative (DNp53; Fig. 2C). Compared with the gene cluster response to JUND knockdown or constitutive E2F1 activation through RB inhibition with overexpressed human papillomavirus E7 protein, we observed substantially more alterations upon NRF2 knockdown (66%) or inhibition of p53 (31%; Fig. 2D and fig. S2, A to D). Using public chromatin immunoprecipitation sequencing (ChIP-seq) datasets (50, 51), we found significant enrichment of proximal NRF2 binding among transcripts reduced by NRF2 knockdown and a slight enrichment in p53 binding among those increased by NRF2 knockdown (fig. S2C). Compound perturbation of NRF2 and p53 elicited further nonadditive changes to multiple genes in the cluster, including synergistic reduction in CDKN1A, encoding a cyclin-dependent kinase inhibitor, and KRT5, encoding a basal cytokeratin. Although p53 can antagonize certain NRF2 target genes in reporter assays (52), significant antagonism was detected for only one transcript in the cluster (MRPL33; fig. S2C). Phenotypically, disruption of NRF2 reduced mean 3D growth by 10 to 13% (fig. S3, A to D), but dual perturbation with p53 gave rise to an increase in aborted spheroids unable to grow in the culture (Fig. 2E). The penetrance of the phenotype (37%; range, 34 to 44%) was close to the percentage of cells showing stabilized NRF2 at the same time point in 3D culture (43%; Fig. 1E). For this clonal basal-like breast epithelial line (20), we conclude that 3D culture heterogeneously elicits NRF2- and p53-inducing stresses, which must be withstood for extended proliferation.
NRF2 disruption in basal-like premalignancy causes similar p53 adaptations but different 3D phenotypes
We next asked how the cellular, molecular, and phenotypic relationships between NRF2 and p53 change in basal-like premalignancy by using isogenic MCF10DCIS.com cells (53) as a proxy for DCIS (54). MCF10DCIS.com cells express oncogenic HRAS (55) and hyperproliferate as 3D spheroids (confirmed in fig. S4A), but they retain wild-type p53 function, albeit at reduced levels compared with parental MCF10A cells (fig. S4, B and C). By two-color immunostaining, we found that NRF2-p53 costabilization was even more pronounced in MCF10DCIS.com cells [MI = 0.30 (0.27 to 0.33); Fig. 3A]. To identify common adaptive programs downstream of NRF2 deficiency, we inducibly knocked down NRF2 and profiled 3D spheroids by RNA-seq (see Materials and Methods). Among transcripts consistently increased or decreased in both MCF10A-5E and MCF10DCIS.com spheroids, there was a significant enrichment in gene signatures encompassing p53, including transcriptional programs downstream of BRCA1, ATM, and CHEK2 (Fig. 3B and data file S3). Consistent with these results, NRF2 knockdown in MCF10DCIS.com cells was sufficient to significantly stabilize p53 (fig. S5A). Stabilization of wild-type p53 upon NRF2 knockdown was also observed in premalignant CHEK21100delC SUM102PT cells (56) and became even more pronounced when these cells were reconstituted with inducible wild-type CHEK2 (fig. S5, B and C), as expected, given the feedforward stabilization of p53 by ATM and ATM-activated CHEK2 (57). Thus, NRF2 impairment promotes p53 pathway activity in basal-like breast epithelia without the need for specific oncogenic drivers.
(A) Quantitative immunofluorescence of NRF2 and p53 abundance in ECM-attached MCF10DCIS.com cells grown as 3D spheroids. Median-scaled two-color average fluorescence intensities are quantified along with the log-scaled and background-subtracted MI with 90% CI for n = 1832 cells segmented from 70 to 110 spheroids from two separate 3D cultures. (B) Common changes in transcript abundance identified by RNA-seq of MCF10A-5E (5E) and MCF10DCIS.com (DCIS.com) cells grown as 3D spheroids with or without NRF2 knockdown. The negative control for shNRF2 was an inducible shGFP (5E) or shLacZ (DCIS.com). Data are log2-transformed Z scores for genes detected at >5 transcripts per million from n = 4 biological replicates. Enriched gene sets for the BRCA1, ATM, and CHEK2 networks are indicated, with black denoting multiple enrichments. The complete list of enrichments is available in data file S3. (C) Quantification of rounded spheroids (circularity >0.9) in 3D-cultured MCF10DCIS.com cells with or without NRF2 knockdown. The negative control for shNRF2 was an inducible shLacZ. (D) Quantification of large spheroids (size > e8.5 ≈ 5000 μm2) in 3D-cultured MCF10DCIS.com cells with or without p53 disruption. The negative control for p53 constructs was an inducible FLAG-tagged LacZ. (E) Quantification of size and circularity in 3D-cultured MCF10DCIS.com cells with or without NRF2 knockdown, p53 disruption, or both. For (C) to (E), cells with or without inducible perturbations were treated with doxycycline (1 μg/ml) for 48 hours, grown as 3D spheroids for 10 days, imaged by brightfield microscopy, and segmented. For (C) and (D), data are mean ± 90% bootstrap-estimated CI from n = 8 biological replicates, with P values by rank sum test estimated by bootstrapping. For (E), data are means ± SEM of n = 8 biological replicates. Statistical interaction between NRF2 and p53 perturbations (Pint) was assessed by two-way ANOVA with replication. Scale bars, 100 μm.
Despite many transcriptomic alterations in common with MCF10A-5E cells (Fig. 3B), MCF10DCIS.com cells yielded very different 3D phenotypes when NRF2 or p53 was perturbed. NRF2 knockdown did not detectably alter 3D growth (fig. S6A) but instead gave rise to more rounded, organized MCF10DCIS.com spheroids of high circularity compared with control (Fig. 3C), which reverted upon addback of an RNA interference (RNAi)–resistant (RR) NRF2 mutant (fig. S6B). NRF2 deficiency also increased rounding in 3D cultures of SUM102PT cells with or without CHEK2 reconstitution (fig. S6C). By contrast, p53 disruption in MCF10DCIS.com cells with either DNp53 or a gain-of-function p53R280K mutant increased the prevalence of hyper-enlarged outgrowths (Fig. 3D). Combined NRF2-p53 perturbation elicited a synergistic increase in nonspherical hyper-enlargement (Fig. 3E), starkly contrasting the proliferative suppression observed with the same combination in nontransformed MCF10A-5E cells (Fig. 2E). The data suggested that the coordinate transcriptional adaptations of NRF2 and p53 are conserved in premalignant cells but insufficient to buffer the cellular phenotypes caused by single-gene perturbations in either pathway.
NRF2 and p53 are coordinately stabilized by sporadic oxidative stress
Coordination of the NRF2-p53 pathways could be achieved if they shared the same inducer. We thus considered various potential upstream and intermediate triggers for NRF2 and p53 stabilization in basal-like breast epithelia. Inhibition of KEAP1 with the electrophile sulforaphane (58) stabilized NRF2 but not p53, and pharmacologic inhibition of MDM2 with nutlin-3 (59) stabilized p53 but not NRF2 (fig. S4, B to E), suggesting they act as parallel pathways downstream of a common inducer. An obvious candidate was DNA damage, given CDKN1A and MUS81 in the gene cluster (Fig. 1A, top) and the most recognized function of p53 (60). However, chemotherapy-induced double-strand breaks did not appreciably stabilize NRF2 in cells with wild-type p53 (Fig. 4A and fig. S7, A and B), and genetically driving increased proliferation (61) did not detectably affect regulation of the gene cluster in 3D spheroids (fig. S2, B and D). The lack of NRF2-p53 coinduction by conventional agonists prompted a search for less canonical activators.
(A and B) NRF2 and p53 stabilization by oxidative stress compared with DNA double-strand breaks. MCF10A-5E cells were treated with 5 μM doxorubicin (double-strand breaks) or 200 μM H2O2 (oxidative stress) for the indicated time points, and NRF2 (magenta) or p53 (green) protein abundance was estimated by quantitative immunoblotting. Data are means ± SEM of n = 3 (A) or 4 (B) independent perturbations. n.s., not significant. (C) Endogenous oxidative stress association with NRF2 stabilization in 3D spheroids. MCF10A-5E cells stably expressing HyPer-2 (67) and mRFP1-NRF2 reporter (NRF2rep) were grown as 3D spheroids for 10 days and imaged by laser scanning confocal microscopy. Representative pseudocolored images for HyPer-2 ratio (top left) and mRFP1-NRF2 reporter (bottom left) are shown. HyPer-2 ratios and mRFP1-NRF2 reporter fluorescence are quantified (right) along with the log-scaled MI with 90% CI for n = 605 cells segmented from 10 to 25 spheroids from four separate 3D cultures. (D) Suppression of endogenous NRF2-p53 coordination during 3D culture with the antioxidant Trolox. Representative pseudocolored images for NRF2 (top left) and p53 (middle left) are shown merged with DAPI nuclear counterstain (bottom left). White arrows indicate concurrent NRF2 and p53 stabilization. The log-scaled and background-subtracted MI (right) is shown with 90% CI estimated from n = 1000 bootstrap replicates. (E) Trolox interference with the synergistic proliferative suppression caused by dual inactivation of NRF2 and p53 in MCF10A-5E cells. Data are mean percentage of proliferation-suppressed spheroids ± SEM of n = 8 independent 3D-cultured samples after 10 days. The overall effect of Trolox on spheroid size is shown in fig. S10. Statistical interaction between Trolox and NRF2-p53 (Pint) was assessed by three-way ANOVA with replication. For (A) and (B), change in protein abundance over time was assessed by one-way ANOVA. For (D) and (E), MCF10A-5E cells cultured for 10 days in 3D with or without 50 μM Trolox supplemented every 2 days. Scale bars, 10 μm (C) and 20 μm (D).
One shared inducer of the KEAP1-NRF2 and ATM-CHEK2-p53 pathways is oxidative stress (62, 63). In human breast tissue, increased levels of reactive oxygen species (ROS) are generated and tolerated by basoluminal progenitors (64), which are the cells of origin for basal-like breast cancer (65). We documented local niches of Nrf2 stabilization in the murine mammary gland during puberty (fig. S8, A to F), potentially linking NRF2 and oxidative stress in expanding progenitor(-like) cells, such as MCF10A. When MCF10A-5E cells were exogenously stimulated with H2O2, NRF2 was rapidly stabilized, and p53 also accumulated after several hours (Fig. 4B and fig. S7, A and B). Recognizing oxidative stress heterogeneities in 3D spheroids (21, 22, 66), we used the genetically encoded sensor HyPer-2 (67) together with an engineered mRFP1-NRF2 reporter (NRF2rep) to colocalize intracellular H2O2 with stabilized NRF2 (see Materials and Methods and fig. S9, A to F). We observed a small but nonzero MI between HyPer-2 fluorescence ratios and NRF2rep [MI = 0.05 (0.02 to 0.10); randomized MI = 0.0004 (0.0001 to 0.0007); Fig. 4C], suggesting a complex connection between the two reporters (see next section). Next, we evaluated whether oxidative stress resided upstream of NRF2-p53 coordination by using the cell-permeable, vitamin E analog Trolox to quench overall ROS in the 3D cultures. Trolox treatment halved the MI between stabilized NRF2-p53 and significantly reduced the synergistic proliferative suppression caused by dual perturbation of NRF2 and p53 (Fig. 4, D and E, and fig. S10). Together, the data suggested that the NRF2 and p53 pathway coregulation involves upstream heterogeneities in oxidative stress.
An integrated NRF2-p53 model of oxidative stress reconciles pathway coordination with 3D phenotypes
To connect NRF2 and p53 costabilization with spontaneous heterogeneities in oxidative stress, we assembled an integrated computational systems model. The model expands or condenses isolated modules of NRF2 and p53 signaling from the literature, fusing them through known or reported mechanisms of cross-talk and convergence (Fig. 5A). For the NRF2 pathway, we streamlined the detailed model of Khalil et al. (68) at several points. Instead of relying on ill-defined kinetic parameters for KEAP1-mediated ubiquitination, KEAP1-NRF2 complexes were modeled as separate oxidized or reduced species with distinct half-lives estimated by experiment (see Materials and Methods). We likewise abandoned the elaborate multistep encoding of thioredoxins, peroxiredoxins, and glutathione transferases (68) by substituting a simpler, lumped pool of antioxidant enzymes in the model. The resulting architecture is similar to the general negative-feedback control scheme of stress response gene regulatory networks described by Zhang and Andersen (69). Last, we retained the nucleocytoplasmic trafficking of stabilized NRF2 to account for observations that H2O2 stimulation retains NRF2 in the cytoplasm longer than treatment with the electrophilic stress, sulforaphane (fig. S11, A to C). Oxidative stress feeds directly into the NRF2 module according to a basal production rate of ROS, which was adjusted in the final model to yield steady-state intracellular H2O2 concentrations consistent with the literature (70).
(A) Connecting NRF2 and p53 signaling models (68, 69, 71) through oxidative stress activators and antioxidant target enzymes. Additional cross-talk linking oxidative stress to p53 inhibition (73), p53 to NRF2 through p21 (76), and NRF2 to MDM2 (74, 75) was considered (gray). (B) Simulation strategy for quantifying association between signaling intermediates. The model was challenged with various ROS production rates and randomly sampled at multiple intermediate time points (yellow to blue). Integrated intracellular H2O2 (gray) is used for phenotype predictions related to NRF2 and p53 perturbation. (C) Intracellular H2O2 concentration is associated with a reporter of NRF2 stabilization (NRF2rep) following simulated step increases in ROS production rate as illustrated in (B). (D and E) Coordination of NRF2 and p53 stabilization in the oxidative stress model and in simulations of premalignancy through the computational approach illustrated in (B). (F and G) Modeling NRF2 knockdown by reduced synthesis captures the synergistic oxidative stress profile of cells harboring dual perturbation of the NRF2 and p53 pathways. In (F), transcriptional changes secondary to NRF2 knockdown were added to the model according to the results in Fig. 3B. For (C) to (E), simulated time points are log-scaled and background-subtracted MI with 90% CI for 10 time points from n = 100 random ROS generation rates. For (F) and (G), time-integrated intracellular H2O2 profiles are scaled to the unperturbed simulations and reported as the mean oxidative stress with 90% CI from n = 100 random ROS generation rates. Statistical interaction between shNRF2 and DNp53 perturbations (Pint) was assessed by two-way ANOVA with replication.
For the p53 pathway, we built upon the base model of Batchelor et al. (71), which was originally used to describe oscillations in p53 abundance after ionizing radiation. In this model, the kinases ATM and CHEK2 act as aggregate sensor transducers of the DNA damage response (Fig. 5A). They phosphorylate and stabilize p53 against degradation triggered by the ubiquitin ligase MDM2, which is also directly phosphorylated and inactivated by ATM. Stabilized p53 promotes its own degradation by inducing the expression of MDM2 transcripts and deactivates ATM-CHEK2 by enhancing transcription of the phosphatase-encoding gene PPM1D. For the integrated model, oxidative stress replaced DNA double-strand breaks as the pathway trigger, recognizing that ATM autoactivates in the presence of oxidants (63). Furthermore, in response to oxidative stress, proper induction of many antioxidant enzymes requires p53 (72), which contributes to the overall antioxidant pool along with antioxidant response element (ARE) target genes (Fig. 5A). Oxidative stress has also been reported to inhibit p53 DNA binding (73), but we found that p53 stabilized by H2O2 treatment was as capable at increasing MDM2 abundance as was p53 stabilized by nutlin-3 (fig. S12). Likewise, NRF2 increases MDM2 abundance in some settings (Fig. 5A, gray) (74, 75), but we were unable to detect changes in MDM2 when NRF2 was knocked down with shRNA or stabilized with sulforaphane (fig. S13, A to D). As a final candidate for NRF2-p53 cross-talk that was conditionally incorporated in the model, we considered reports that p21, encoded by the p53 target gene CDKN1A, directly stabilizes NRF2 by interfering with KEAP1-catalyzed turnover (Fig. 5A, gray) (76, 77). Together, the modifications provided an integrated model of NRF2-p53 signaling downstream of oxidative stress with enough molecular detail to enable kinetic and functional predictions.
We revisited the oxidative stress time course (Fig. 4B) to append immunoblot quantification of ATM-CHEK2 phosphorylation and p21 abundance after H2O2 addition (fig. S14A). Exogenous H2O2 was encoded as an extracellular spike-in that decayed rapidly and spontaneously (78) amidst a basal ROS generation rate, yielding a realistic intracellular H2O2 burden at steady state (70). The H2O2 partition coefficient in the model was calibrated to capture the magnitude of NRF2 stabilization (see Materials and Methods). Likewise, the parameters for H2O2-induced autoactivation of ATM-CHEK2 and signal inactivation were defined to align with the time-delayed kinetics and duration of p53 stabilization (fig. S14B). In this model, addition of p53-p21-NRF2 cross-talk (76) caused NRF2 stabilization to peak earlier and deactivate faster than observations (fig. S14C). We were also unable to detect even transient short-range interactions between inducible BirA*-fused versions of p21 or NRF2 and endogenous NRF2 or p21 by proximity ligation (fig. S15, A to C). The results, thus, argued against p53-p21-NRF2 cross-talk during oxidative stress in these cells.
With the provisionally calibrated base model, we sought to test whether the encoded mechanisms of regulation were sufficient to capture prior observations relating NRF2, p53, and oxidative stress. The data obtained by quantitative fluorescence microscopy (Figs. 2B, 3A, and 4C) presumably arose from spontaneous oxidative stress that was occurring transiently and asynchronously during imaging. We mimicked oxidative stress transients by triggering a step increase in the rate of ROS production for 2 hours, followed by relaxation of the system for an additional 10 hours (Fig. 5B). The magnitude of the step was sampled lognormally to elicit intracellular H2O2 concentrations within the range of HyPer-2 ratios observed experimentally (see Materials and Methods). We represented the asynchrony of image acquisition by randomly selecting 10 snapshots of the network for each model iteration. This collection of 1000 snapshots (100 random generation rates × 10 random time points) was used to quantify coordination of species within the model.
For connecting oxidative stress to NRF2 stabilization (Fig. 4C), we expanded the base model to include the mRFP1-NRF2 reporter, which does not bind DNA or interact with MAF proteins and requires ~1 hour to mature fully (see Materials and Methods) (79). By contrast, HyPer-2 becomes fully oxidized within ~1 min of H2O2 addition (67), enabling intracellular H2O2 concentration in the model to be used directly as a surrogate of HyPer-2 fluorescence ratio. We calculated the MI from 1000 simulated snapshots and found that the two reporters were statistically coupled in the model [MI = 0.14 (0.10 to 0.19); randomized MI = 0.001 (2.8 × 10−6 to 0.0005); Fig. 5C]. Associations were stronger than those observed by experiment (Fig. 4C) due to the early time points sampled in the model (yellow-orange times in Fig. 5C), suggesting that peak H2O2 transients may be difficult to observe in practice. We next compared endogenous NRF2-p53 costabilization between MCF10A-5E and MCF10DCIS.com cells. The base MCF10A-5E model was adjusted to reflect (i) proportional differences in species abundance estimated from RNA-seq (see Materials and Methods) and (ii) an increased ROS generation rate estimated from HyPer-2 imaging (fig. S9D). NRF2-p53 MI was much less dependent on signaling transients, and coupling was substantially higher in MCF10DCIS.com cells. The simulations were consistent with our immunofluorescence data (Figs. 3A, 4C, and 5, D and E) and support that NRF2-p53 pathway kinetics were accurately encoded in the base model.
We then investigated whether the base model could also relate to the synergistic phenotypes observed upon dual NRF2-p53 perturbation in MCF10A-5E and MCF10DCIS.com cells (Figs. 2E and 3E). We mimicked shNRF2-mediated knockdown by reducing the NRF2 production rate fivefold in the model (fig. S3B) and encoding secondary transcriptional adaptations in other components by using the associated RNA-seq data (Fig. 3B). For DNp53, the p53 species was rendered unable to induce transcription of MDM2, PPM1D, p21, and its share of the antioxidant enzyme pool. After reestablishing steady state, the perturbed models were challenged with the random step increase in ROS production described above. We used the time-integrated intracellular H2O2 concentration as the overall measure of oxidative stress experienced during simulation with either the MCF10A-5E or the MCF10DCIS.com initial conditions. For both cell lines, the base model predicted synergistic increases in oxidative stress beyond the linear superposition of shNRF2 and DNp53 effects (Fig. 5F). Encouragingly, the same conclusions were reached with models that simply encoded the reduced NRF2 production rate without secondary adaptations (Fig. 5G). Beyond oxidative stress inducers and antioxidant target enzymes, we conclude that the NRF2-p53 network does not require any additional mechanisms to capture signaling coordination or phenotypic interactions.
NRF2-p53 coregulation occurs in normal breast tissue and hormone-negative DCIS but not in invasive TNBC
The regulatory heterogeneities observed in 3D culture often reflect adaptations in hormone-negative premalignancy (24) that become further disrupted in TNBCs (25). We thus sought to quantify NRF2-p53 coordination in TNBC and premalignant DCIS lesions, using adjacent normal tissue as a comparator. The TP53 gene is frequently mutated in TNBC (35) and gives rise to loss of p53 protein or hyperstabilization of a dominant-negative mutant in tumors (80). By contrast, prior immunohistochemistry of NRF2 abundance in breast carcinomas was inconclusive (81), owing to an NRF2 antibody that was later shown to be nonspecific (82). There was an opportunity to revisit NRF2-p53 abundance heterogeneities from the perspective of costabilization, with a focus on TNBC and its precursor lesions.
Using a knockout-verified commercial antibody (83), we immunoblotted with our production lot and confirmed detection of basal and induced NRF2 with only ~35% immunoreactivity attributed to nonspecific bands (fig. S16, A to C). By immunohistochemistry, the antibody-detected endogenous NRF2 stabilized with electrophiles in paraffin sections of cell pellets (fig. S16D). The antibody has also been used independently to track NRF2 abundance in other solid tumors (84). However, when we stained adjacent normal epithelium immunohistochemically, NRF2 was not clearly discernible (Fig. 6A, top). In MCF10A 3D spheroids, stabilized p53 is not detected by immunohistochemistry either (85), and yet we readily visualized it by immunofluorescence (Fig. 2B). Therefore, to improve signal-to-background and facilitate multiplex quantification, we used two-color immunofluorescence after antigen retrieval, segmenting 24,949 normal and transformed epithelial cells in 15 cases of TNBC and hormone-negative DCIS.
(A) Immunohistochemistry (top) and immunofluorescence (bottom) for NRF2 and p53 in tumor-adjacent normal breast lobules. Hematoxylin and eosin (H+E, top right) histology is from a serial paraffin section for p53. Images from a tumor-adjacent normal breast duct are shown in fig. S17. (B and C) Multicolor immunofluorescence for NRF2 and p53 in (B) hormone-negative DCIS and (C) TNBC. (D) Quantification of the association between NRF2 and p53 immunoreactivities represented in (A) to (C). (E and F) Median NRF2 and p53 immunoreactivities for the designated tissue type in each clinical case. n.s., not significant (P > 0.05). For (A) to (C), immunofluorescence is shown as representative pseudocolored images for NRF2 (left) and p53 (middle) are shown merged with DAPI nuclear counterstain (right). White arrows indicate concurrent NRF2 and p53 stabilization, and magenta or green arrows indicate stabilization of NRF2 or p53 separately. Scale bars, 20 μm. For (D) to (F) data are means ± SEM of n = 14 cases with tumor-adjacent normal epithelium (Normal), 8 cases with DCIS, and 7 cases with TNBC. Multigroup comparison was made by Kruskal-Wallis rank sum test with Šidák correction for multiple hypothesis testing.
In adjacent normal epithelium, we observed local niches of stabilized NRF2 in lobules and ducts, which often corresponded with stabilized p53 (Fig. 6A, bottom; and fig. S17). Stabilized NRF2 was frequently detected in the cytoplasm, consistent with the prolonged cytoplasmic localization observed in H2O2-treated cells compared with cells stressed with an electrophile (fig. S11). The results corroborated findings that KEAP1 senses oxidative stress differently than electrophilic stress (62). The patterns of NRF2-p53 coaccumulation were largely preserved in hormone-negative DCIS (Fig. 6B and fig. S18, A to C), even in cases with abundantly stabilized p53 that was likely mutated (see later in this section). Nuclear localization of NRF2 was also more prominent, perhaps reflecting the stronger ROS generation rates of transformed cells (86). NRF2 and p53 were almost completely uncoupled in invasive TNBCs (Fig. 6C and fig. S18D), reflecting a profound shift in single-cell regulation. We quantified NRF2-p53 coordination by MI and found that it was largely eliminated in regions of invasive TNBC, irrespective of whether p53 was chronically stabilized (Fig. 6D). Such alterations were not apparent in regional estimates of protein abundance by cell population–averaged fluorescence, where neither NRF2 nor p53 was reproducibly different among groups (Fig. 6, E and F). We conclude that 3D culture in reconstituted basement membrane costimulates the NRF2-p53 pathways akin to that observed in normal breast tissue and hormone-negative premalignancy. Full-blown TNBC, by contrast, evokes a different set of dependencies.
TNBC adaptations to p53 disruption predict variable NRF2 miscoordination, NRF2-deficient oxidative stress profiles, and 3D growth responses
TP53 is the most frequently mutated gene in TNBC (35), and transcriptomic analyses support it as a prevalent founder mutation in the disease (87). Disrupting p53 would undoubtedly affect transcriptional feedback and the overall cellular response to oxidative stress (Fig. 5A). Conversely, neither NFE2L2 nor KEAP1 is mutated in breast cancer (88), but it is unclear whether wild-type NRF2 might serve as a transient “non-oncogene” (89) that promotes stress tolerance during early tumorigenesis. Compared with in situ lesions, the stromal environment of invasive tumors is stiffer and more mesenchymal (90), which may render NRF2 signaling dispensable at later stages. We wondered whether the fragmentation of the NRF2-p53 network in TNBC cells and its origins could be reconciled with the systems model.
Using RNA-seq data from the National Institutes of Health (NIH) Library of Integrated Network-Based Cellular Signatures (LINCS) consortium (91) on 15 TNBC lines with mutated p53 (six claudin-low subtype and nine basal-like subtype), we adjusted initial conditions from the original MCF10A model and removed all transcriptional processes downstream of p53 (Fig. 7A; see Materials and Methods). The individual TNBC models were run to steady state and then challenged with increased ROS generation rates as in Fig. 5B. The coordination between NRF2 and mutant p53 was calculated by MI, and the integrated H2O2 response was scaled to that of MCF10A-5E cells as a relative measure of ROS tolerance. The goal was to associate the model-derived predictions with NRF2-knockdown phenotype in ROS-generating environments such as 3D culture. To the extent possible, we hoped that 3D growth in reconstituted basement membrane might quantify any vestigial requirements for NRF2 signaling from the in situ stage of the TNBC lines.
(A) Transcripts per million for the indicated TNBC cell lines scaled to MCF10A cells from the NIH LINCS dataset (91). (B) NRF2-p53 MI and ROS tolerance for TNBC cell lines using the simulation strategy in Fig. 5B. (C to G) Quantification of mean spheroid area with or without NRF2 knockdown in 3D-cultured TNBC cells with higher simulated ROS tolerance and NRF2-p53 MI (C and D) and with lower simulated ROS tolerance and NRF2-p53 MI (E to G). (H) Transcripts per million for the TNBC lines in (A) scaled to clinical cases of TNBC in TCGA (35). (I) Simulated NRF2-p53 MI and ROS tolerance for TNBC tumors. Vertical lines indicate cases with high MI in the lower quartile of ROS tolerance. For (A) and (H), the clustered transcripts were used to adjust the initial conditions of the model simulations for each cell line and tumor. For (B) and (I), ROS tolerance was defined as the integrated intracellular H2O2 concentration in each cell line compared with that of MCF10A-5E cells in response to an increased ROS production rate as in Fig. 5B. For (C) to (G), TNBC cells with or without inducible NRF2 knockdown were treated with doxycycline (1 μg/ml) for 72 hours, grown as 3D spheroids, imaged by brightfield microscopy, and segmented. Data are means ± SEM of n = 4 to 8 biological replicates. The difference between means was assessed by Student’s t test with Šidák correction for multiple hypothesis testing, and the specific p53 mutation of each line is shown (bottom left).
For all TNBC lines, the model predicted substantially reduced covariation between mutant p53 and NRF2 compared with MCF10DCIS.com cells with wild-type p53 (MI < 0.25; Fig. 7B). We noted a spectrum of residual NRF2-p53 costabilization from weak (HCC1937 and SUM159PT) to virtually nonexistent (MDA-MB-468 and MDA-MB-231). Despite complete p53 deficiency in the model, this residual NRF2-p53 MI correlated strongly with the simulated relative increase in oxidative stress when ROS generation rate was increased (Fig. 7B). Neither of these predictions mapped directly to specific transcripts in the TNBC-specific RNA-seq data (Fig. 7A), reinforcing that the models were making nonobvious predictions about oxidative stress handling.
To connect the model predictions with a continued role for NRF2 signaling in TNBC behavior, we selected five lines along the spectrum of MI and ROS tolerance. HCC1937 and SUM159PT cells were both predicted to have residual NRF2-p53 coordination and moderate ROS tolerance (Fig. 7B). Accordingly, inducible knockdown of NRF2 in these lines did not lead to any consistent changes in 3D growth (Fig. 7, C and D). By contrast, MDA-MB-231, HCC1806, and MDA-MB-468 cells were predicted to have among the least NRF2-p53 costabilization and ROS tolerance (Fig. 7B). Knockdown of NRF2 in these lines with two different shRNAs caused significant increases or decreases in overall cell growth (Fig. 7, E to G). Thus, model and experiment support that, despite p53 mutation, residual NRF2-p53 coupling indicates the primordial susceptibility of triple-negative cell lines to perturbations in the NRF2 pathway.
Last, we sought to extend model predictions to 122 cases of TNBC sequenced by The Cancer Genome Atlas (TCGA) (35). Compared with the TNBC lines, we noted reduced abundance of MDM2, the NRF2 binding partner MAFK, ATM, and CHEK2 (Fig. 7H), which suggested that TCGA tumors would be a considerable deviation from prior simulations. There was also more variability in the abundance of multiple antioxidant genes (SOD1, TXN, and PRDX1), anticipating a greater breadth in model outcomes. Unexpectedly, when tumor-derived profiles were encoded and simulated (see Materials and Methods), the models predicted ROS tolerances that were largely within the range of TNBC lines analyzed before (Fig. 7I). The associated NRF2-p53 coordination, by contrast, was qualitatively different, with various TCGA cases giving rise to strong coordination despite pervasive TP53 mutation (Fig. 7I, purple). The high-coordination, low-tolerance TNBC cases (Fig. 7I, vertical lines) form a subset that could be especially sensitive to changes in NRF2 activation.
DISCUSSION
This work posits ROS as an endogenous, spatially heterogeneous trigger of dual NRF2-p53 activation in breast-mammary epithelia surrounded by basement membrane ECM. NRF2 and p53 regulate target gene abundance—both cooperatively and independently—to promote stress tolerance and adaptation. NRF2 deficits are buffered by compensatory increases in p53 signaling, and notable ROS-dependent phenotypes arise when both pathways are perturbed. In hormone-negative premalignant lesions, stabilization of NRF2-p53 remains coordinated, even in cases where p53 has likely mutated. At this preinvasive stage, NRF2 should be most important for tumorigenesis. After invasion through basement membrane and progression to TNBC, the stromal microenvironment reduces the overall NRF2 signaling and often uncouples it from (now mutant) p53. Here, the effect of activating or inhibiting NRF2 will depend more on the exact cellular context and, thus, be unpredictable (for example, Fig. 7, E to G). Despite the overall complexity of NRF2- and p53-mediated transcriptional programs (92, 93), the coordinated response to oxidative stress is captured by a relatively simple mathematical encoding. Known core mechanisms of NRF2-p53 regulation are brought together by a shared ROS inducer and a common pool of detoxifying target genes without the need for any further cross-talk. Therefore, oxidative stress handling in normal breast-mammary epithelia is usefully abstracted as two stability-regulated TFs working independently toward a common homeostatic goal.
Although NRF2 is not an oncogene for breast cancer, it has been connected with multiple breast cancer tumor suppressors previously. In mouse mammary epithelial cells, loss of Brca1 (a predisposing event for basal-subtype TNBC) destabilizes Nrf2 and causes an increase in ROS, favoring the future acquisition of p53 mutations (36, 94). In human breast cancer cells, gain-of-function p53 mutants interact directly with NRF2 and may help retain NRF2 in the nucleus (37). If certain p53 mutations were also to promote NRF2 stabilization, then it would provide a two-for-one benefit to cancer progression by relieving tumor suppression and conferring ROS tolerance constitutively. However, we did not note any association between gain-of-function p53 mutants and NRF2 abundance in TNBC lines, suggesting that KEAP1 regulation predominates, as indicated by the TNBC models. Chronic activation of the NRF2 pathway (for example, by activating NFE2L2 mutation or KEAP1 loss) may be disfavored if increased intracellular ROS is not permanent. The models suggest that supraphysiological activation of NRF2 would lead to runaway induction of antioxidant enzymes, causing reductive stress as documented for NRF2 in other tissues (95). Wild-type NRF2 function must be sufficient to buffer cells from the early stresses of premalignancy and p53 disruption, allowing invasive TNBCs to deactivate the pathway when it is no longer needed. There are parallels to FOXO TFs (96), which are reversibly inactivated by mitogenic signals yet provide critical oxidative stress tolerance when the breast cancer tumor suppressor RUNX1 is disrupted (22, 97, 98).
Breast cancer cell lines organize very differently in 3D culture (99), but their response to perturbations is often less disparate. For example, gain-of-function p53 mutations cause luminal filling in MCF10A 3D cultures (48), similar to the delay in mammary gland involution observed with mutant p53 in vivo (100). Reciprocally, knockdown of mutant p53 in MDA-MB-468 cells promotes luminal hollowing (101). Among p53-mutant TNBC lines, the impact of NRF2 knockdown on 3D growth was nonuniform but explainable through the stress profiles inferred from TNBC-specific systems models. The balance of complexity and tractability makes 3D spheroid-organoid cultures a compelling platform for systems-level dissection of cell state heterogeneity and early tumorigenesis.
The 3D behavior of breast-mammary cancer cells is highly dependent on the surrounding ECM (102). Invasive cancers no longer encounter basement membrane ECM but must have bypassed it upon progression to carcinoma. Although multiple TNBC lines will grow as 3D colonies in reconstituted basement membrane, others cannot, suggesting a type of cellular “amnesia” toward that past encounter. For cancers that do grow in 3D, the use of reconstituted basement membrane (as a more normal-like microenvironment) may give rise to cellular changes reminiscent of premalignancy. We exploited these changes to evaluate the relative importance of NRF2 signaling in different TNBC backgrounds. There are likely other opportunities to examine hurdles of premalignancy by using basement membrane 3D cultures. For 3D organoid biobanks (19), however, it is a reminder that such cultures are not propagating the primary breast tumor but rather tumor-derived cells in a more primitive state.
Cancer mutations engage and cooperate with cell signaling in ways that are not captured by DNA sequencing (103). The coupling of the NRF2 and p53 pathways described here provides a robust oxidative stress–handling network for glandular morphogenesis and maintenance. However, this same coupling creates a redundancy upon which p53 mutations can occur and neoplasms can evolve. Our results give pause to the nutraceutical use of sulforaphane as a potent NRF2 stabilizer (104)—in lung cancer, where KEAP1-NRF2 mutations are common and TP53 mutation is secondary, antioxidants accelerate tumor progression (105). The extraordinary complexity of ROS generation and its cellular effects reinforce the value of modeling redox networks at a granularity suited to a given physiology or pathology (106).
MATERIALS AND METHODS
Plasmids
shRNA targeting sequences from the RNAi consortium (107) were cloned into tet-pLKO.1-puro as previously described (38) for shLuc (TRCN0000072250, Addgene #136587), shNRF2 #1 (TRCN0000281950, Addgene #136584), shNRF2 #2 (TRCN0000284998, Addgene #136585), shJUND #1 (TRCN0000416347, Addgene #136581), and shJUND #2 (TRCN0000416920, Addgene #136583).
For the mRFP1-NRF2 reporter (Addgene #136580), the DNA binding domain of NRF2 was mutated (C506S) along with four leucines (L4A) in the leucine zipper region of the bZIP (basic leucine zipper) domain by site-directed mutagenesis of the pBabe mRFP1-NRF2 hygro plasmid (Addgene #136579) originally prepared by subcloning into pBabe mRFP1 hygro. The RR version of NRF2 (Addgene #136522) was prepared by introducing four silent mutations into the sequence targeted by shNRF2 #1 in pEN_TT 3xFLAG-NRF2 (Addgene #136527). Site-directed mutagenesis was performed with the QuikChange II XL kit (Agilent).
pDONR223 CHEK2 was obtained from the human Orfeome V5.1 (108). CHEK2 amplicon was prepared with Xba I and Mfe I restriction sites and cloned into pEN_TTmiRc2 3xFLAG (Addgene #83274) that had been digested with Spe I and Mfe I (Addgene #136526). BirA* was cloned out of pcDNA3.1 mycBioID (Addgene) (109) with Xba I and Spe I restriction sites and cloned into pEN_TTmiRc2 digested with Spe I and Mfe I (Addgene #136521). CDKN1A and NRF2 PCR amplicons were prepared with Spe I and Mfe I restriction sites and cloned into pEN_TTmiRc2 BirA* (Addgene #136521). Luciferase PCR amplicon was prepared with Spe I and Eco RI restriction sites and cloned into pEN_TTmiRc2 3xFLAG digested with Spe I and Mfe I sites (Addgene #136519). p53DD (p53DN) and p53(R280K)-V5 PCR amplicon was prepared with Spe I and Mfe I restriction sites and cloned into pEN_TTmiRc2 (Addgene #25752), digested with Spe I and Mfe I (Addgene #136520 and #136525).
pEN_TT donor vectors were recombined into pSLIK neo (Addgene #25735), pSLIK zeo (Addgene #25736), or pSLIK hygro (Addgene #25737) by LR recombination to obtain pSLIK 3xFLAG-Luciferase zeo (Addgene #136533), pSLIK p53DD zeo (Addgene #136534), pSLIK 3xFLAG-Luciferase hygro (Addgene #136528), pSLIK 3xFLAG-NRF2(RR) hygro (Addgene #136535), pSLIK BirA* hygro (Addgene #136537), pSLIK BirA*-CDKN1A hygro (Addgene #136538), pSLIK BirA*-NRF2 hygro (Addgene #136539), pSLIK p53(R280K)-V5 hygro (Addgene #136540), and pSLIK 3xFLAG-CHEK2 neo (Addgene #136536).
pLXSN HPV16E7 (110) and the ∆DLYC mutant (Addgene #136588) were provided by S. Vande Pol (University of Virginia). pCDH–HyPer-2–puro (66) was provided by J. Brugge (Harvard Medical School).
Cell lines
The MCF10A-5E clone was previously reported and cultured as described for MCF-10A cells (13, 20). MCF10DCIS.com cells were obtained from Wayne State University and cultured in Dulbecco’s modified Eagle’s medium/F-12 medium (Gibco) plus 5% horse serum (Gibco). SUM102PT cells were obtained from Asterand Biosciences and cultured in Ham’s F-12 (Gibco) plus 10 mM Hepes (Gibco), epidermal growth factor (10 ng/ml; PeproTech), 5 mM ethanolamine (Sigma-Aldrich), 50 nM sodium selenite (Sigma-Aldrich), apo-Transferrin (5 μg/ml; Sigma-Aldrich), 10 nM triiodo-l-thyronine (VWR), insulin (5 μg/ml; Sigma-Aldrich), hydrocortisone (1 μg/ml; Sigma-Aldrich), and 5% fatty acid–free bovine serum albumin (VWR). SUM159PT cells were obtained from Asterand Biosciences and cultured in Ham’s F-12 (Gibco) plus 10 mM Hepes (Gibco), insulin (5 μg/ml; Sigma-Aldrich), hydrocortisone (1 μg/ml; Sigma-Aldrich), and 5% fetal bovine serum (Hyclone). All other cell lines were obtained directly from the American Type Culture Collection (ATCC). MDA-MB-231 and MDA-MB-468 cells were cultured in L-15 medium plus 10% fetal bovine serum without supplemental CO2. HCC1806 and HCC1937 cells were cultured in RPMI 1640 medium plus 10% fetal bovine serum. All base media were further supplemented with 1× penicillin and streptomycin (Gibco). All cell lines are female, were grown at 37°C, authenticated by short tandem repeat profiling by ATCC, and confirmed negative for mycoplasma contamination.
Viral transduction and selection
Lentiviruses were prepared in human embryonic kidney 293 T cells (ATCC) by triple transfection of the viral vector with psPAX2 + pMD.2G (Addgene) and transduced into MCF10A-5E, MCF10DCIS.com, HCC1937, SUM159PT, MDA-MB-231, HCC1806, and MDA-MB-468 as previously described (25). Retroviruses were prepared similarly by double transfection of the viral vector with pCL ampho (Addgene) and transduced into MCF10A-5E cells as previously described (22). Transduced cells were selected in growth medium containing puromycin (2 μg/ml), G418 (300 μg/ml), hygromycin (100 μg/ml), or zeocin (25 μg/ml) until control plates had cleared. For RR addback, viral titers were adjusted to match the endogenous protein abundance as closely as possible. For mRFP1-NRF2 fluorescent reporter, we used the minimum viral titer that gave sufficient signal in sulforaphane-treated cells compared with dimethyl sulfoxide (DMSO)–treated cells. Reporter abundance was roughly equal to endogenous NRF2 expression (fig. S9F).
3D culture
3D overlay cultures were performed on top of Matrigel (BD Biosciences) as described previously for MCF-10A cells (111) with culture media previously optimized for each cell line (25). In addition, HCC1806 cells were cultured in MCF10A assay media (111), and SUM159PT cells were cultured in SUM159PT growth media (described above) plus 2% fetal bovine serum. For each culture, 45 μl of Matrigel was spread with a pipette tip on the bottom of an eight-well chamber slide. A suspension of 5000 single cells per well was laid on top of the Matrigel in culture media supplemented with 2% Matrigel. 3D culture medium was replaced every 4 days as originally described (111). For antioxidant supplementation, cells were treated with 50 μM Trolox (Calbiochem) for 2 days before 3D culture, and Trolox was included in media refeeds and supplemented every 2 days between refeeds. For long-term knockdown experiments, cells were treated with doxycycline (1 μg/ml; Sigma-Aldrich) for 3 days before 3D culture, and doxycycline was maintained in the 3D culture medium throughout the experiment. For experiments with long-term knockdown and inducible overexpression, cells were treated with doxycycline (1 μg/ml) for 2 days before 3D culture, and doxycycline was maintained in the 3D culture medium throughout the experiment.
RNA purification
RNA from cultured cells was isolated with the RNeasy Plus Mini Kit (QIAGEN) according to the manufacturer’s protocol. RNA from 3D cultures at day 10 was extracted by lysing individual wells in 500 μl of RNA STAT-60 (Tel-Test) and purified as described previously (25).
RNA-seq and analysis
Total RNA was diluted to 50 ng/μl and prepared using the TruSeq Stranded mRNA Library Preparation Kit (Illumina). Samples were sequenced on a NextSeq 500 instrument with NextSeq 500/550 High Output v2.5 kits (Illumina) to obtain 75–base pair (bp) paired-end reads at an average depth of 15 million reads per sample. Adapters were trimmed using fastq-mcf in the EAutils package (version ea-utils.1.1.2-537) with the following options: -q 10 -t 0.01 -k 0 (quality threshold 10, 0.01% occurrence frequency, and no nucleotide skew causing cycle removal). Quality checks were performed with FastQC (version 0.11.7) and multiqc (version 1.5). Datasets were aligned to the human (GRCh38.86) genome using HISAT2 with the option: --rna-strandness RF (for paired-end reads generated by the TruSeq strand–specific library). Alignments were assembled into transcripts using StringTie (version 1.3.4) with the reference guided option. Transcripts that were expressed at greater than five transcripts per million across all samples were retained for downstream analysis. Differential gene expression analysis was carried out using edgeR (version 3.8) (112) on raw read counts corresponding to transcripts that passed the abundance-filtering step. Trimmed mean of M values normalization using the calcNormFactors function was performed before differential expression analysis using exactTest in edgeR. The 1132 transcripts that were commonly differentially expressed [5% false discovery rate (FDR)] between MCF10A-5E shControl and shNRF2 #1, shControl and shNRF2 #2, and MCF10DCIS.com shControl and shNRF2 #1 are shown in Fig. 3B. Gene set enrichment analysis was done on transcripts that were differentially increased or decreased in shNRF2 compared with shControl using the Molecular Signatures database collections C1, C2, C3, C4, C6, and C7 (113, 114). The full list of enrichments (5% FDR) is provided in data file S2.
Quantitative PCR
cDNA synthesis and qPCR were performed as previously described (25, 115) with the primers listed in table S1. Human samples were normalized to the geometric mean of ACTB, HINT1, PP1A, and TBP (Fig. 2D and fig. S2C); B2M, GAPDH, GUSB, HINT1, and PRDX6 (fig. S2A); or ACTB, B2M, GUSB, PPIA, and PRDX6 (fig. S2B).
Brightfield imaging and quantification of spheroid phenotypes
Brightfield 3D images were acquired on an Olympus CKX41 inverted microscope with a 4× plan objective (four fields per chamber) and a qColor3 camera (Q-Imaging). Images were segmented using OrganoSeg (116) to produce morphometric measures for each segmented spheroid. “Rounded” spheres were classified as having circularity greater than 0.9 (Fig. 3, C and E; and fig. S6, B and C). “Hyper-enlarged” spheres were classified as having an area greater than e8.5 ≈ 5000 µm2 (Fig. 3, D and E). “Proliferation suppressed” spheres were classified as having an area less than 1600 μm2 for MCF10A-5E cells after 10 days of 3D culture (Figs. 2E and 4E).
Clinical samples
Cases were identified from the pathology archives at the University of Virginia and build upon a cohort of samples previously described (24, 25). Hormone-negative DCIS lesions were deemed negative (less than 10% expression frequency) for estrogen receptor and progesterone receptor by clinical immunohistochemistry, and TNBC cases were additionally scored negative for HER2 amplification by clinical DNA chromogenic in situ hybridization. All clinical work was done according to the Institutional Review Board for Health Sciences Research approval #14176 and Protocol Review Committee approval #1363 (502-09).
Immunofluorescence
MCF10A-5E and MCF10DCIS.com 3D cultures were embedded at day 10 of morphogenesis, and 5-μm sections were cut and mounted on Superfrost Plus slides (Fisher Scientific). For clinical samples, paraffin tissue sections were dewaxed and antigens were retrieved on a PT Link (Dako) with low-pH EnVision FLEX Target Retrieval Solution (Dako) for 20 min at 97°C. Immunofluorescence on cryosections and antigen-retrieved slides was performed as previously described (20) with the following primary antibodies: NRF2 (1:100; Santa Cruz Biotechnology, #sc-13032), phospho-Rb (1:1600; Cell Signaling, #8516), HIF-1α (1:800; Cell Signaling, #79233), and p53 (1:200; Santa Cruz Biotechnology, #sc-126). Slides were incubated the next day for 1 hour in the following secondary antibodies: Alexa Fluor 555–conjugated goat anti-rabbit (1:200; Invitrogen) and Alexa Fluor 647–conjugated goat anti-mouse (1:200; Invitrogen).
Image acquisition analysis and MI calculation
Fluorescence images were collected unblinded on an Olympus BX51 fluorescence microscope with a 40× 1.3 numerical aperture (NA) UPlanFL oil immersion objective and an Orca R2 charge-coupled device (CCD) camera (Hamamatsu) with no binning. Images were segmented in CellProfiler (117) using 4′,6-diamidino-2-phenylindole (DAPI) to identify nuclei. Nuclear objects were dilated to a median diameter of 15 μm to capture about one whole cell. NRF2 staining was quantified in the nucleus, the whole cell, and the cytoplasm (whole cell area − nuclear area). p53 staining was quantified in the whole cell. Immunoreactivity was quantified as the median fluorescence intensity of the whole cell unless otherwise noted.
For pRB and NRF2 immunofluorescence (Fig. 1, D and E), log-transformed distributions were analyzed with the MClust function in R using the unequal variance model with either one or two mixture components specified. Model fit was evaluated by F test.
MCF10A-5E cells stably expressing pCDH–HyPer-2–puro were imaged at 37°C in Hanks’ balanced salt solution (Gibco) with a 40× 1.3 NA EC plan Neofluar oil immersion objective on a Zeiss LSM 700 laser scanning confocal microscope. Lasers (405 and 488 nm) were used to sequentially excite two excitation peaks of HyPer-2 and collect fluorescence emission from 500 to 550 nm. To calculate HyPer-2 ratios on a pixel-by-pixel basis, 488-nm images were divided by 405-nm images and thresholded in ImageJ to remove background pixel values (~10%). For quantification of cells cultured in 2D (fig. S9, B to D), the mean HyPer-2 ratio per image was used for analysis. For quantification of cells cultured as spheroids (Fig. 4C), cells were manually segmented to calculate the median HyPer-2 ratio per cell.
Clinical samples were imaged on an Olympus BX51 fluorescence microscope with a 40× 1.3 NA UPlanFL oil immersion objective and an Orca R2 CCD camera (Hamamatsu) with 2 × 2 binning and fixed exposure times for NRF2 (150 ms) and p53 (50 ms). Images were autoexposed in the DAPI channel for nuclear segmentation and in the unlabeled fluorescein isothiocyanate (FITC) channel for autofluorescence estimation. Image fields were classified as follows: normal—bilayered epithelium, intact basement membrane (visualized by FITC autofluorescence), and normal cytoarchitecture; DCIS—multilayered and disorganized epithelium (with partial or complete luminal filling), intact basement membrane, and cytologic atypia; and TNBC—invasive carcinoma cells with cytologic atypia and no discernable basement membrane. All images were segmented in CellProfiler as described above. After nuclear identification, nuclei outside of the ductal epithelium (fibroblasts, endothelial cells, and immune cells) were manually removed using the IdentifyObjectManually module. Because paraffin fixation of tissue increases autofluorescence (118), the analysis excluded images that were dominated by autofluorescent bleedthrough into the Alexa 555 channel localizing NRF2. Spearman correlation was calculated between cellular FITC-555 channels and FITC-DAPI channels on a pixel-by-pixel basis for each image. Images with a FITC-555 correlation coefficient above the 95th percentile for FITC-DAPI correlation (in which autofluorescent artifacts were negligible due to the low exposure time) were excluded from further analysis.
For NRF2 quantification in neighboring cells (fig. S8), spheroid and mouse mammary gland images were loaded into CellProfiler, and the IdentifyObjectManually module was used to manually identify regions of ductal epithelium. The images were cropped manually, and cell nuclei within the cropped area were identified by DAPI staining. Nuclear area was dilated to a median diameter of ~15 μm to define a cell. Position, area, and median NRF2 staining intensity were measured for each cell. Measurements were loaded into MATLAB, and single-cell NRF2 intensities were normalized to the median intensity of all exposure-matched cells. Neighboring cells were defined as cells located within a radius of 1.5 times the median cell diameter. For more distant neighbors, annular areas of 3 to 5 and 5 to 10 times the median cell diameter were used. Cells that fell within the applied search area were used to calculate the median neighbor NRF2 intensity. The original cell at the center of the search area was not included in the intensity calculations, and cells with NRF2 intensity values equal to 0 or lacking neighboring cells within the defined search area were excluded from calculations.
To quantify the association between fluorescence channels, we used MI in lieu of standard correlation measures (Pearson and Spearman). After appropriate transformation and binning into discrete high-low states, MI provides greater flexibility to capture nonlinear relationships (119) and more stringency to detect compressions in dynamic range (120). Median fluorescence intensity distributions were transformed by their respective cumulative distribution functions (probability integral transform) to produce uniformly distributed random variables (121). The uniform distributions were split into low and high states at the 67th percentile, and the joint marginal state probabilities estimated for the two fluorescence channels (R and G) were used to calculate the MI as follows
MI confidence intervals were estimated by bootstrapping the segmented cell population 1000 times. To create a randomized (null) dataset, the values of one fluorescence channel were randomly shuffled before analysis.
Clinical samples often had fewer areas of classified cells for imaging, which require an added analysis step in the MI calculation. For a classification (normal, DCIS, and TNBC) composed of two images from one case, we evaluated batch effects by hypergeometric test to determine if the two images separated by high versus low staining intensity. If so, the case for that classification was excluded.
Quantitative immunoblotting
Quantitative immunoblotting was performed as previously described (122). Primary antibodies recognizing the following proteins or epitopes were used: NRF2 (1:1000; Santa Cruz Biotechnology, #sc-13032), p53 (1:1000; Santa Cruz Biotechnology, #sc-126), p21 (1:1000; ProteinTech, #10355-1-AP), total Chk2 (1:1000; Cell Signaling, #2662), phospho-Chk2 (Thr68; 1:1000; Cell Signaling, #2197), phospho-ATM (Ser1981; 1:1000; Abcam, #ab81292), KEAP1 (1:1000; Santa Cruz Biotechnology, #sc-15246), CDK4 (1:1000; Cell Signaling, #12790), CDK2 (1:200; Santa Cruz Biotechnology, #sc-6248), vinculin (1:10,000; Millipore, #05-386), glyceraldehyde-3-phosphate dehydrogenase (GAPDH; 1:20,000; Ambion, #AM4300), tubulin (1:20,000; Abcam, #ab89984), p38 (1:5000; Santa Cruz Biotechnology, #sc-535), and Hsp90 (1:5000; Santa Cruz Biotechnology, #sc-7947).
Proximity ligation using BirA* fusions of p21 and NRF2
MCF10A-5E cells inducibly expressing the promiscuous biotin ligase BirA* (123), BirA*-NRF2, or BirA*-CDKN1A were plated on 10-cm plates and induced with doxycycline (1 μg/ml) at 50% confluency. After 24 hours, medium was refed with doxycycline (1 μg/ml), 10 μM sulforaphane (Sigma-Aldrich), 10 μM nutlin-3 (Calbiochem), and 1 mM biotin (Sigma-Aldrich). After 24 hours, cells were lysed in 200 μl of radioimmunoprecipitation assay (RIPA) buffer [50 mM tris (pH 8.0), 150 mM NaCl, 5 mM EDTA, 1% Triton X-100, 0.1% SDS, and 0.5% sodium deoxycholate]. Anti-biotin antibody enrichment of biotinylated peptides was performed as previously described (124). Briefly, biotin antibody bound agarose beads (ImmuneChem Pharmaceuticals Inc., #ICP0615) were washed three times in immunoaffinity purification (IAP) buffer [50 mM MOPS (pH 7.2), 10 mM sodium phosphate, and 50 mM NaCl]. Antibody (500 μg; 50 μl) was added to each RIPA lysate on ice. Ice-cold IAP buffer was added up to 1 ml, and samples were incubated on a nutator overnight at 4°C. The next day, beads were washed four times with ice-cold IAP buffer, boiled in dithiothreitol-containing 2× Laemmli sample buffer, and used for immunoblotting against the indicated targets.
Promoter bioinformatics
The 36 transcripts of the Fig. 1A gene cluster (20, 24) were assessed with four promoter analysis algorithms to identify recurrent TF candidates (125). First, distant regulatory elements (DiRE) analysis was conducted using the DiRE website (https://dire.dcode.org) (40) searching evolutionarily conserved 5′ untranslated regions (5′ UTR ECRs) and evolutionarily conserved promoter regions (promoter ECRs) for genes on the human genome (hg18). A random set of 7500 genes was selected as background control genes. Second, Expression2Kinases (X2K) software was used to identify upstream TFs for the Fig. 1A gene cluster (41). The potential TFs were selected from ChIP-X Enrichment Analysis (ChEA) database using “mouse + human” as the background organisms (126). The P value from the Fisher’s exact test and Z score were used for sorting and ranking. Third, from the National Center for Biotechnology Information (NCBI), we collected the proximal promoter of each transcript—defined as 1416 bp upstream and 250 bp downstream of the transcription start site to remain within the 60-kb sequence limit—for use as an input set for MEME (127, 128). Using MEME-defined motifs from classic discovery mode, the top three enriched motifs were searched against the JASPAR CORE (2018) database (containing 1404 defined TF binding sites for eukaryotes) (129) or HOCOMOCO (Homo Sapiens Comprehensive Model Collection) Human (v11) database (containing 769 TF binding motifs) (130) using Tomtom (131) to identify TF recognition sequences. A recognition sequence for E2F6 was considered within the E2F group, and a recognition sequence for the NRF2 binding partner MAFK was considered within the NRF2 group. Last, oPOSSUM (43) was used to identify potential TFs targeting transcripts in the cluster. We selected Single Site Analysis–Human mode and used all 24,752 genes in the oPOSSUM database as a background. All vertebrate profiles with a minimum specificity of eight bits in the JASPAR CORE Profiles were selected as TF binding site sources. oPOSSUM was run with the following parameters: conservation cutoff of 0.4, matrix score threshold of 85%, amount of upstream/downstream sequence: 2000/0, and sort results by Fisher score. Outputs of the promoter bioinformatics are available in data file S1.
ChIP-seq bioinformatics
NRF2 ChIP-seq raw data files were downloaded from ENCODE (ENCAB800OND) (50), consisting of fastq files from three cell lines (K562, A549, and HepG2) with two biological replicates each. Quality of the sequenced reads was analyzed using FastQC. Reads were aligned to the human genome (hg19) using BWA with the -M option. Peaks were identified using MACS2 (version 2.1.0) with an FDR cutoff of 0.01 to reduce the number of spurious peaks. Irreproducibility discovery rate analysis was performed on biological replicates, and a cutoff of 0.05 was used to generate a list of high-confidence peaks for each cell line. Peaks were annotated using the Homer annotatePeaks program. p53 ChIP-seq binding sites were used from a ChIP-seq dataset (GSE86164) (51). Briefly, ChIP-seq was performed on three cell lines (HCT116, MCF7, and SJSA) treated with or without nutlin-3. Reads were mapped to hg19 using Bowtie2, and peaks were identified and annotated using the Homer suite. NRF2-p53 binding sites were indicated (fig. S2C) if a peak was present for a gene in at least one cell line analyzed for each ChIP-seq dataset.
Computational modeling
The NRF2 pathway was encoded as first- and second-order rate equations for KEAP1 oxidation and NRF2 stabilization; NRF2-mediated transcription of antioxidant enzymes was modeled as a Hill function (68, 69). The p53 pathway was reconstructed from a delay differential equation model of p53 signaling in response to DNA damage (71). Abundances in the original p53 model were unitless, but abundances were cast as concentrations in the earlier NRF2 models. Consequently, the integrated model adopted unitless abundances in its initial conditions and second-order parameters (table S2). To adapt the p53 DNA damage model to respond to oxidative stress, we changed the “Signal” activation (representing activation of upstream kinases p-ATM and p-CHEK2) from a Heaviside step function to a first-order oxidation reaction of ATM/CHEK2 by intracellular H2O2 (63). A basal ROS generation rate was added yielding a realistic intracellular H2O2 burden at steady state (70). Transcription of antioxidant enzymes by p53 (72) was modeled using the same model parameters describing the p53-mediated induction of MDM2 (71). p53- and NRF2-mediated antioxidant gene transcription contributes to a shared pool of antioxidant enzymes, which catalytically reverse the oxidation states of KEAP1 and p-ATM/CHEK2. Transcription of CDKN1A by p53 (132) was included for model calibration (fig. S14) and for testing the relevance of p53-p21-NRF2 cross-talk (see next paragraph). The integrated base model of NRF2-p53 oxidative stress signaling contains 42 reactions and 22 ordinary differential equations (ODEs). The model was simulated with dde23 in MATLAB to reach steady state before the addition of oxidative stress.
The integrated model was calibrated to capture the dynamics of MCF10A-5E cells stimulated with 200 μM H2O2 (fig. S14). Bolus addition of H2O2 was simulated as an impulse of intracellular H2O2. We used an H2O2 partition coefficient that gave rise to NRF2 stabilization levels comparable to immunoblot quantification (extracellular/intracellular partition = 3). We approximated p-ATM/CHEK2 in the integrated model as the maximum normalized increase in p-ATM or p-CHEK2 over baseline at each experimental time point. Robustness of the system output to initial conditions was evaluated by randomly varying the concentration of model species with a log coefficient of variation of 10%, taking the base model as the geometric mean.
For simulations involving the mRFP1-NRF2 reporter (NRF2rep; Fig. 5C), NRF2rep and mature fluorescent species (Nrf2repmat) were added to the MCF10A-5E base model. Both reporter species were allowed to react with KEAP1, but neither could bind MAF proteins or antioxidant response elements in the model (fig. S9E). We used an mRFP1 maturation time of 1 hour (79) to model the conversion of NRF2rep to NRF2repmat. The modifications added 19 additional reactions and eight additional ODEs to the MCF10A-5E base model.
For simulations involving p53-p21-NRF2 cross-talk (fig. S14C), we added reactions involving p21 binding to NRF2 to the MCF10A-5E base model. p21 was assumed to interact with NRF2 like KEAP1 and compete with KEAP1 for binding NRF2 through its DLG and ETGE domains (76). The p21-NRF2 complex was assumed to degrade at the same reduced rate as when NRF2 is bound to oxidized KEAP1 (k_nrf2degox). These modifications added eight additional reactions and two additional ODEs to the MCF10A-5E base model.
For simulations involving MCF10DCIS.com cells (Fig. 5, E to G), RNA-seq data (Fig. 3B) were used to estimate proportional differences in model species abundance between MCF10DCIS.com and MCF10A-5E cells. Average gene expression in transcripts per million from the four biological replicates of MCF10DCIS.com and MCF10A-5E control cell lines was calculated for each gene. Fold changes in model species of MCF10DCIS.com relative to MCF10A-5E were used to adjust each initial condition in the model. For the “MAF” species, we used the median fold change in NRF2-binding small MAFs MAFF, MAFG, and MAFK (133). For the antioxidant species, we used the median fold change in TXN, SOD1, PRDX1, and HMOX1 to include antioxidants that react with both free radicals and oxidized proteins (134). In addition, the MCF10DCIS.com model included a 1.4-fold increase in the basal ROS generation rate, informed by the increased median HyPer-2 ratio in MCF10DCIS.com cells compared with MCF10A-5E cells (fig. S9D). The increased ROS generation rate was paired with an increased basal turnover of the antioxidant pool to arrive at steady-state antioxidant gene expression levels consistent with MCF10DCIS.com RNA-seq data.
For simulations involving bursts of oxidative stress, an increased ROS production rate was added for 2 hours to match the duration of transient stabilizations of JUND (a gene in the NRF2-associated gene cluster) in 3D (24). We selected the minimum increase in ROS generation (20-fold; log coefficient of variation = 20%) that gave rise to a detectable stabilization of both the NRF2 and p53 pathways in the MCF10A-5E base model. Under these conditions, overall oxidative stress burden was within 4 to 16% of that observed with 200 μM H2O2. For MCF10DCIS.com and TNBC models, the mean ROS generation rate was scaled 1.4-fold to reflect the increased basal ROS generation rate described above. NRF2 knockdown was encoded by decreasing the net synthesis rate of NRF2 fivefold to mimic the fivefold decrease in NRF2 protein resulting from short-hairpin knockdown (fig. S3B). To account for secondary transcriptional adaptations (Fig. 5F), initial conditions were also adjusted by RNA-seq–based fold changes in model species for shNRF2 cells relative to negative control cells (Fig. 3B). DNp53 was encoded by removing all reactions downstream of p53 (transcriptional activation of MDM2, PPM1D, CDKN1A, and the p53 share of the antioxidant enzyme pool).
For the control case and all genetic perturbations (shNRF2, DNp53, and shNRF2 + DNp53), 100 simulations were run with random ROS generation rates varied with a log coefficient of variation of 25% to capture the variability of HyPer-2 ratios observed experimentally (Fig. 4C). Each simulation was run for 2 hours with increased ROS production rate and then an additional 10 hours to allow relaxation back to steady state. For assessment of species coordination (Figs. 5, C to E, and 7, B and I), species abundances were captured at 10 random time points from each simulation, and MI was calculated as it was for quantitative immunofluorescence datasets. For oxidative stress analysis (Figs. 5, F and G, and 7, B and I), the time-integrated intracellular H2O2 concentration was used as an overall measure of oxidative stress.
For simulations involving TNBC cells (Fig. 7, A and B), RNA-seq data from the NIH LINCS consortium (91) (Harvard Medical School dataset ID: 20348) was used to estimate proportional differences in model species abundance between 15 TNBC cell lines and MCF10A cells. Reads per kilobase per million mapped reads values were normalized as transcripts per million before fold change calculation. MAF and antioxidant species were estimated as described above. TNBC models used the same increased basal ROS generation rate as in the MCF10DCIS.com model (135). To simulate p53 mutation in the 15 p53-mutant TNBC cell lines, all reactions downstream of p53 were removed.
For simulations involving TNBC tumors from TCGA (Fig. 7, H and I), breast cancer sequencing data and associated clinical information were downloaded from TCGA Data Portal (https://portal.gdc.cancer.gov/). We identified 122 TNBC cases as tumors that were scored negative for estrogen receptor expression, progesterone receptor expression, and HER2 amplification in the clinical record. Fragments per kilobase per million mapped reads were normalized as transcripts per million before fold change calculation to model species abundance relative to MCF10A cells. Simulations were carried out exactly as described for TNBC cell lines. All computational models and associated results are available in data file S4.
Statistical analysis
For analysis of the 10cRNA-seq dataset (Fig. 2A), Spearman correlation between transcripts and the median expression of the NRF2-associated gene cluster was calculated at an FDR of 10%. Transcripts with a Spearman correlation coefficient above 0.5 were examined by GO analysis. Statistical enrichment of GO terms was assessed by Fisher’s exact test with FDR-corrected P values. Statistical enrichments in ChIP-seq binding were determined by hypergeometric test (fig. S2C). For qPCR data, differences in geometric means were assessed by Welch’s t test after log transformation (fig. S2, A and B). Statistical interaction between shNRF2 and DNp53 and differences between immunoblotting time courses were assessed by one-way analysis of variance (ANOVA) (Figs. 2, D and E; 3E; 4, A and B; and 5, F and G; and fig. S2C). Statistical interaction between shNRF2, DNp53, and Trolox was assessed by three-way ANOVA (Fig. 4E). For fig. S15B, two-way ANOVA without replication was used. For unpaired clinical data, multigroup comparison was made by Kruskal-Wallis rank sum test (Fig. 6, D to F). For 3D spheroid growth, mean differences in area were assessed by Kruskal-Wallis test with Dunn’s post hoc test (Fig. 7, C to G). Distributions were compared by Kolmogorov-Smirnov test (figs. S2D, S9, and S11). All other two-sample comparisons were performed by Student’s t test.
SUPPLEMENTARY MATERIALS
stke.sciencemag.org/cgi/content/full/13/627/eaba4200/DC1
Fig. S1. HIF-1α is not appreciably stabilized in 3D culture.
Fig. S2. Abundance of the heterogeneously regulated gene cluster is perturbed by NRF2 knockdown or p53 disruption, but not by JUND knockdown or human papillomavirus E7–induced inhibition of RB.
Fig. S3. NRF2 knockdown and 3D phenotype quantification in MCF10A-5E cells.
Fig. S4. Proliferation differences and signaling similarities between MCF10A-5E and MCF10DCIS.com cells.
Fig. S5. NRF2 knockdown causes p53 stabilization in premalignant breast epithelial cell lines.
Fig. S6. Premalignant breast epithelial cell lines have similar adaptations to NRF2 knockdown in spheroid culture.
Fig. S7. Representative immunoblot images for the double-strand break and oxidative stress time courses in MCF10A-5E cells.
Fig. S8. Local niches of NRF2 stabilization in MCF10A-5E 3D spheroids and pubertal murine mammary glands.
Fig. S9. Description and validation of the HyPer-2 probe for H2O2 and the mRFP1-NRF2 reporter.
Fig. S10. Antioxidant treatment causes an overall increase in MCF10A-5E spheroid size.
Fig. S11. Oxidative stress stabilizes NRF2 in the cytoplasm more so than electrophilic stress.
Fig. S12. Oxidative stress does not measurably inhibit MDM2 induction by stabilized p53.
Fig. S13. NRF2 perturbations do not detectably alter MDM2 abundance.
Fig. S14. Calibration of an integrated NRF2-p53 systems model for oxidative stress.
Fig. S15. Endogenous NRF2 and p21 are not proximity labeled by BirA* fusions of each other.
Fig. S16. Anti-NRF2 antibody validation for immunohistochemistry.
Fig. S17. NRF2 and p53 are costabilized in breast epithelial ducts.
Fig. S18. Low-magnification hematoxylin-eosin images of the tissues and tumors in the work.
Table S1. qPCR primer sequences.
Table S2. Parameter summary for the integrated NRF2-p53 computational model.
Data file S1. Promoter analysis results underlying the summary Venn diagram in Fig. 1B.
Data file S2. GO enrichment analysis.
Data file S3. Gene set enrichment analysis of differentially abundant transcripts in MCF10A-5E and MCF10DCIS.com cells upon NRF2 knockdown compared with control.
Data file S4. NRF2-p53 computational model and associated files.
This is an article distributed under the terms of the Science Journals Default License.
