Mapping the tumor stress network reveals dynamic shifts in the stromal oxidative stress response

Abstract

The tumor microenvironment (TME) presents cells with challenges such as variable pH, hypoxia, and free radicals, triggering stress responses that affect cancer progression. In this study, we examine the stress response landscape in four carcinomas-breast, pancreas, ovary, and prostate-across five pathways: heat shock, oxidative stress, hypoxia, DNA damage, and unfolded protein stress. Using a combination of experimental and computational methods, we create an atlas of stress responses across various types of carcinomas. We find that stress responses vary within the TME and are especially active near cancer cells. Focusing on the non-immune stroma we find, across tumor types, that NRF2 and the oxidative stress response are distinctly activated in immune-regulatory cancer-associated fibroblasts and in a unique subset of cancer-associated pericytes. Our study thus provides an interactome of stress responses in cancer, offering ways to intersect survival pathways within the tumor, and advance cancer therapy.

Keywords: CP: Cancer; NRF2; cancer; cancer-associated fibroblasts; fibroblasts; oxidative stress; pericytes; scRNA-seq; stress responses; stroma; tumor microenvrionemnt.

Graphical abstract

Figure 1

Stromal stress response activation increases with proximity to cancer cells Formalin-fixed paraffin-embedded (FFPE) tumor microarrays (TMAs) of breast (n = 57), pancreas (n = 71), ovary (n = 102), and prostate (n = 43) cancer patients were stained by MxIF for the indicated proteins. DAPI was used to stain nuclei. (A) Representative images are shown. Scale bars = 50 μm. White arrows mark stressed stromal cells; yellow arrows mark stromal cells that are negative for stress markers. (B) Images were analyzed using QuPath software. CD45⁻ CK⁻ cells were defined as non-immune stromal cells and stratified into stressed and unstressed cells based on staining for either NRF2, HSF1, or ATF4. The distance of each non-immune stromal cell to its nearest cancer cell was calculated and averaged. (C) For each patient, the Pearson correlation coefficient between the intensity of the indicated protein of non-immune stromal cells and the distance to the nearest cancer cell was calculated. Results are shown as mean ± SD. The p values were calculated using Student’s t test (B, paired t test; C, one-sample t test; μ = 0).

Figure 2

scRNA-seq data analysis uncovers shared and unique stress response patterns across organs scRNA-seq datasets from human breast, pancreas, ovary, and prostate tumors were reanalyzed using the Seurat R toolkit. (A–D) UMAP plots of 265,034 cells from 51 breast cancer patients^, (A), 199,938 cells from 59 pancreatic cancer patients^,, (B), 84,369 cells from 20 ovarian cancer patients^,, (C), and 32,823 cells from 13 prostate cancer patients (D). UMAPs are colored by cell type, defined by differential gene expression and canonical cell type markers. (E–I) Quantification of stress scores per cell type for each tumor across patients: (E) OSR, (F) HySR, (G) HSR, (H) UPR, and (I) DDR. Different letters denote significant differences in stress scores as determined by ANOVA followed by Tukey's HSD (honestly significant difference) test. Groups with the same letter are not significantly different from each other.

Figure 3

Correlation analysis of stress signatures reveals coordinated activation of the HSR in prostate tumors, the DDR in pancreatic tumors, and the UPR in ovarian tumors (A–D) Correlation matrix of stress scores of different cell types across patients calculated from the scRNA-seq datasets listed in Figure 2. Per-patient average scores were quantified, and Pearson coefficients of all possible pairs were calculated. Outlier patients were removed to avoid bias. Color bars indicate the stress or cell type.

Figure 4

The hypoxia and OSR are differentially activated across subpopulations of the non-immune tumor stroma (A–D) UMAP plots of 20,754 fibroblasts and pericytes from 51 breast cancer patients (A), 14,516 fibroblasts and pericytes from 57 pancreas cancer patients (B), 10,762 fibroblasts and pericytes from 20 ovarian cancer patients (C), and 1,697 fibroblasts and pericytes from 12 prostate cancer patients (D). UMAPs are colored by cell type, defined by the gene signatures we defined (Table S1). The scRNA-seq data originate from the same datasets highlighted in Figures 2E–2I. Shown is projection of the five stress signatures scores: (E) OSR, (F) HySR, (G) HSR, (H) UPR, and (I) DDR.

Figure 5

NRF2 and the OSR are associated with the iCAF signature (A and B) Pearson correlation coefficients between stress and cell type scores calculated from the scRNA-seq data described in Figure 4: (A) iCAFs, (B) musclePer. Results are shown as mean ± SD. The p values were calculated using one-sample t test (μ = 0). (C and D) FFPE TMAs of breast (n = 114), pancreas (n = 125), ovary (n = 93), and prostate (n = 104) cancer patients were stained by MxIF and analyzed using QuPath software. CD45⁻ CK⁻ cells were defined as non-immune stromal cells. Representative images are shown (C). Scale bars, 40 μM. For each patient, Pearson correlation coefficients between the staining intensities of NRF2 and the different CAF markers CLU, ɑSMA, and MHC class II were calculated (D).Results are shown as mean ± SD. The p values were calculated using one-way ANOVA, followed by Tukey’s multiple comparisons test.

Figure 6

NRF2 and the OSR contribute to the transition of normal fibroblasts to iCAFs (A–E) Patient data from the TCGA breast cancer dataset were analyzed for cellular composition using the CIBERSORTx algorithm, and the results were ordered by the stress score. n = 1,187. (F–O) Trajectory analyses. (F and K) UMAP plots of scRNA-seq data of breast (E) or pancreas (J) fibroblasts with data from normal samples, re-analyzed from publicly available datasets.^, (G and L) Clusters were annotated based on the stromal gene signatures we described in Figure S3. (H and M) Trajectory analyses of breast and pancreas tumors and normal samples using the Monocle3 R toolkit. (I and N) Projections of the OSR score we defined previously. (J and O) Patient-level quantification of OSR scores. The p values were calculated using one-way ANOVA, followed by Tukey’s multiple comparisons test. (P and Q) Kaplan-Meier analysis of overall survival for low-grade breast cancer patients from the METABRIC cohort. Patients were stratified based on their OSR signature and iCAF/myCAF (P) or musclePer/matriPer (Q) ratio, calculated by CIBERSORTx (median was used as cutoff). The p values were calculated from a log rank test, and paired comparisons were calculated using the Survdiff function in R with false discovery rate (FDR) correction. (R) Heatmap of the differentially expressed oxidative stress-related genes from bulk RNA-seq between cell culture models of iCAFs, myCAFs, and quiescent PSCs, re-analyzed from Öhlund et al. Differentially expressed genes were filtered by a logFC (log fold change) threshold of 1.5 and adjusted p value of 0.05. (S) Immortalized PSCs were seeded in Matrigel for 3 days with either DMEM or KPC organoid conditioned medium (orgCM), and the OSR genes Ho1 and Sod3 were measured using qPCR. (T) Immortalized PSCs were seeded in 2D culture and depleted of Nrf2 using small interfering RNA (siRNA). Cells were then seeded in Matrigel for 3 days with orgCM, and known iCAF markers in this system, Clu and C3, were measured using qPCR. Results are shown as mean ± SD. The p values were calculated using two samples t test. Results are from 2 to 4 independent experiments.