Tamoxifen Response at Single-Cell Resolution in Estrogen Receptor–Positive Primary Human Breast Tumors

Breast cancer is the most common cancer in women, and around 600,000 women die from breast cancer worldwide each year (1). Approximately 75% of breast cancers are clinically positive for the estrogen receptor–positive (ER⁺) and categorized in the luminal subtypes. In general, ER⁺/luminal breast cancers respond to endocrine therapy targeting ER (2, 3). However, up to a third of early ER⁺ breast cancers are inherently unresponsive or develop resistance to endocrine therapy, and resistance eventually develops in all patients with metastatic disease (4, 5).

Resistance to endocrine therapy is complex and multifactorial (6). Known mechanisms of resistance include autologous activation of ER signaling (7–9), alterations of transcription factors (10, 11), and upregulation of growth factor signaling pathways (12–15). However, these account for only a minority of cases, and mechanisms have not been described for the majority of endocrine-resistant tumors.

Signatures of estrogen-response genes defined from the ER⁺ cell line MCF-7 are highly prognostic in ER⁺ breast cancer patients (16), demonstrating the importance of ER transcriptional activity on tumor biology and outcomes. Similarly, in patients treated with endocrine therapy, gene-expression changes, including reduced proliferation, correlate strongly with outcome (17–19). However, bulk profiling averages the signal obtained from multiple cell types including stromal, immune, and normal epithelial cells, which are relatively abundant in ER⁺/HER2⁻ tumors (20). Tumor heterogeneity within human breast cancers is associated with metastasis and worse response to treatment, including endocrine therapy (21–24). Distinct response to treatments has been described within breast tumors in cell populations with unique gene expression profiles (25, 26). Therefore, bulk analysis that includes cell types with distinct responses could mask relevant tumor-specific response elements and low-abundance subpopulations that may contribute to resistance. In addition, nontumor cell types with distinct responses could limit precision in defining tumor-specific transcriptional changes. Improved understanding of how distinct populations within human tumors respond to endocrine therapy, as well as the heterogeneity in response within tumors, is needed to reveal underlying mechanisms of resistance.

In this study, we sought to determine how subpopulations of cells within normal breast epithelium and ER⁺/HER2⁻ breast tumors respond to tamoxifen treatment using an ex vivo treatment model coupled to single-cell RNA sequencing (scRNA-seq). We hypothesize that subpopulations within ER⁺/HER2⁻ human breast tumors have distinct transcriptional responses to tamoxifen and that discerning heterogeneity in cellular response will identify patients at risk for treatment failure.

Patients were identified with ER⁺/HER2⁻ breast cancer or without breast cancer undergoing reduction mammaplasty. Patients received no preoperative therapy. Written informed consent was obtained from all patients under a protocol that was approved by the Institutional Review Board at the University of North Carolina at Chapel Hill in accordance with the U.S. Common Rule. Consent included access to deidentified patient data, which was obtained through an honest broker. Participants were not compensated for participation.

Primary tissue specimens were obtained in the operating room and placed in DMEM/F12 (Gibco) media with 1% Penicillin–Streptomycin (Life Technologies). Specimens were transferred immediately to the laboratory and processed into single-cell suspensions by sharply mincing into 2- to 4-mm fragments. Enzymatic dissociation was then performed using Gentle Collagenase/Hyaluronidase (Stemcell Technologies Inc. 07919) in DMEM/F12 supplemented with 5% BSA (Sigma Aldrich), Hydrocortisone (Stemcell Technologies Inc.), HEPES (Corning), and Glutamax (Gibco) for 16 hours at 37°C with cell agitation. The cells were gently centrifuged and washed twice with PBS (Gibco) supplemented with FBS and HEPES buffer. Cells were resuspended in cold ammonium chloride solution (Stemcell Technologies Inc. 07800) and incubated at room temperature to remove red blood cells. Cells were centrifuged and briefly trypsinized in warm 0.05% Trypsin-EDTA (Gibco) and DNase I (Stemcell Technologies Inc. 07900). Cells were further dissociated with dispase (Stemcell Technologies Inc. 07923) and DNase I. Cells were centrifuged and washed then resuspended in DMEM/F12 with 10% FBS and 1% Penicillin–Streptomycin. Cells were then strained using multiple rounds of sequential straining with 40-μm mini strainers (Pluriselect USA Inc.) to remove cell debris.

We obtained the ER⁺/HER2⁻ cell line T47D from the ATCC. T47D cultures were maintained at 37°C with 5% CO₂ in DMEM with 10% FBS and 1% Penicillin–Streptomycin.

Cells were counted using trypan blue (Invitrogen) and the Countess II cell counter (Life Technologies). Dissociated single cells (n = 100,000) were placed in suspension on an Ultra-Low attachment microplate (Corning 3474) in DMEM/F12 media supplemented with 10 nmol/L estradiol (Sigma Aldrich). Paired suspensions were treated with control media or media containing 10 μmol/L 4-OH tamoxifen, which is a commonly used in vitro dose but up to 300X higher than circulating levels in patients treated with tamoxifen. Suspensions were incubated at 37°C with 5% CO₂ for 12 hours. Graphical representation of workflow from tumor resection to data analysis is shown in Fig. 1A.

scRNA-seq was done with the 10x Genomics Chromium Next GEM Single Cell 3′ v3.1 reagent kits according to the manufacturer's revision D user guide. After 12-hour suspension, viability was assessed using trypan blue, quantified using the Countess, and 7,000 to 10,000 viable cells were targeted per library. Libraries were sequenced on the NextSeq500 and NextSeq2000 (Illumina) platforms with pair-end sequencing and single indexing according to the recommended manufacturer's protocols.

Formalin-fixed, paraffin-embedded (FFPE) blocks were reviewed for high tumor concentrations and 10-μm unstained slides cut. RNA was isolated from the FFPE tissue by use of the High Pure FFPE RNA Micro Kit (Roche 04823125001). Slides were incubated in Hemo-De (Scientific Safety Solvents HD-150) and then washed in 100% and 70% ethanol. The tissue was then air-dried at 55°C. Tissue Lysis Buffer and Proteinase K were then added to the tissue pellet and incubated at 55°C. Binding Buffer and ethanol were then added to the tissue pellet, and the supernatant was filtered out via filter tube. DNase Solution and DNase Incubation Buffer were then added to the filter tube and incubated at room temperature. A series of washes were then performed using Wash Buffer I and Wash Buffer II. Elution Buffer was then added to the filter tube, incubated at room temperature, and RNA was then eluted.

RNA-seq reads were aligned to the human genome hg38 (gencode_v36) from Genomic Data Commons using STAR v2.7.6a. Transcripts were quantified using Salmon v1.4.0. The matrices were upper quartile fixed normalized (UQN) where quantile (0.75) = 1,000. We used the HER2/ER subgroup-specific gene-centering method described by Fernandez-Martinez and colleagues (27) using log₂-transformed UQN data and IHC status. The expression values of the PAM50 genes were then HER2/ER-subgroup-specific gene-centered, and the PAM50 predictor (28) was applied to assign subtype calls using centroid correlation values.

Slides from clinical FFPE blocks were obtained for all tumors. IHC for Ki67 antigen was performed on FFPE tissue sectioned at 4 μm. Staining was performed using the Leica Bond III Autostainer system. Slides were dewaxed in Bond Dewax solution (AR9222) and hydrated in Bond Wash solution (AR9590). Heat-induced antigen retrieval was performed at 100°C in Bond-Epitope Retrieval solution 1 pH-6.0 (AR9961). After pretreatment, slides were incubated with Ki67 Antibody (MIB-1, Dako) at 1:100 for 30 minutes followed by Novolink Polymer (RE7260-K) secondary. Antibody detection with 3,3′-diaminobenzidine was performed using the Bond Intense R detection system (DS9263). Ki67 was scored by a breast pathologist (B.C. Calhoun) according to the Ki67 IHC MIB-1 pharmDx (Dako Omnis) Interpretation Manual for Breast Carcinoma. The Ki67 pharmDx score (%) was calculated as number of Ki67 staining viable invasive tumor cells divided by the total number of viable invasive tumor cells (minimum 2,000/slide), multiplied by 100.

The filtered feature barcode matrices were generated by Cell Ranger from 10x Genomics. Data processing and subsequent downstream analyses were divided into four analysis blocks (ABlock1–4). In ABlock1, two paired normal samples were analyzed as tamoxifen/control pairs (four samples in total). ABlock2 involved T47D cells, ABlock3 included 10 pairs of primary tumor samples, and ABlock4 focused on targeted analysis of four pairs of ER⁺ breast cancer tumors that were selected on the basis of the depletion of the scTAM-response-T47D signature. For each analysis block, Seurat objects were built from the filtered feature barcode matrix for samples assigned to the analysis block using the Seurat R package (29, 30).

Outlier cells were defined by three metrics: a low number of UMI counts, a low number of genes expressed, and a high percent mitochondrial read count (>25%). For ABlock1 and ABlock4, we used the following threshold values to identify outliers: the number of UMI counts <2,000 and the number of genes expressed <1,000. For ABlock2 and ABlock3, we used more stringent threshold values using UMI counts <5,000 and the number of genes expressed <2,000 to select cells with higher quality. We used higher thresholding to select higher quality cells to exclude specimens with a lower number of epithelial cells with higher quality. The outlier cells according to these criteria were removed before downstream analyses.

After removing doublets by DoubletFinder (31) to make the final Seurat object for each sample, we applied Seurat's merge() function (32) to combine multiple Seurat objects for each analysis block. The log normalized procedure was applied to the gene expression matrix. Subsequently, a scaling operation was performed on the expression values of the 2,000 most variably expressed genes to facilitate subsequent principal component analysis (PCA). The confounding effect of mitochondrial genes (32) was regressed out. A reduction in dimensionality was performed through the summarization of the 2,000 most variably expressed genes into 50 principal components (PC) via PCA. The cells were then projected onto a uniform manifold approximation and projection (UMAP) embedding via Seurat (30) and ggplot2 (33) R packages. For further clustering analysis, a shared nearest-neighbor graph was constructed using 30 PCs, and Louvain clustering was performed with a resolution of 0.8.

Cell type annotations were performed using the SingleR (34) R package. The normalized expression values obtained from human bulk RNA-seq data sourced from Blueprint and ENCODE were used as reference datasets, available via celldex, with SingleR (34).

To facilitate analysis for each block, a merged dataset was created by integrating the raw count matrices of samples assigned to each analysis block. Subsequently, a canonical correlation analysis was performed using Seurat (32) to correct for batch effects, and cell clustering was performed using the Louvain clustering method. The representative cell type for each cluster in the merged dataset was annotated with a cell type label based on the predominant cell type within each cluster.

To differentiate malignant cells from normal epithelial cells, copy-number events for each cell cluster were estimated using the InferCNV R package (35, 36) using a normal background composed of immune cells and endothelial cells. The cell clusters were specified in the InferCNV annotations file, enabling the estimation of copy-number variations (CNV) at the level of these clusters. The epithelial cells were subsequently stratified into epithelial tumor and epithelial normal cells by plotting scatter plots of the CNV values and the correlation values with the top 5% of cells with high CNV values (37, 38). The tumor epithelial cells were used for downstream analyses for characterizing the molecular subtypes of tumor cells and identifying molecular signatures of tamoxifen resistance.

Molecular subtype estimation was performed by the nearest centroid method using four centroids of LumA, LumB, HER2⁻ enriched, and basal, obtained from scRNA-seq profiles and the molecular subtype classifications of single cells by SCSubtype (38) of human breast cancers. The proliferation score was computed by the averaged normalized expression of 11 genes (38).

The differential gene expression between groups was determined using Seurat's FindMarkers() function (32), with log₂fold change (log₂FC) threshold set at 0.25 and minimum percentage of expressed cells set at 0.25. The markers were then further filtered on the basis of an adjusted P value threshold of less than 0.01. Heat maps were used to visualize the data. The enrichment of cancer hallmark gene sets (39) and a curated list of breast cancer relevant gene signatures (40) were assessed using hypergeometric tests with enricher() function from the clusterProfiler package (41) with q-value threshold of less than or equal to 0.01. Cell-cycle phase for cells within resistant clusters was determined using the Seurat CellCycleScoring() function using cell-cycle marker genes (36).

The clinical significances of gene signatures were estimated by Kaplan–Meier curves. The signature score was defined by the median expression of signature genes. Patients with ER⁺/HER2⁻ tumors that received endocrine therapy in METABRIC (refs. 42, 43; EGAS00000000083) or SCAN-B (refs. 44, 45; GSE202203) were stratified into high/low (by 50th percentile) and high/med/low (by tertiles) groups by signature scores. Because high expression of genes that increase in expression with tamoxifen treatment may be associated with poor response, we generated a resistance signature of upregulated genes in response to tamoxifen in malignant cells (scTAM-resistance-M) from the top 300 upregulated genes ranked by FC after selecting genes with log₂FC > 0.25 and adjusted P < 1e-6, in ABlock4 using our four pairs of tumor samples (Supplementary Table S5). After comparing transcriptomic profiles in tamoxifen-resistant clusters (clusters 3, 12, and 19) with tamoxifen response cluster (cluster 2), we generated three gene signatures of tamoxifen resistance from genes upregulated in tamoxifen-resistant clusters (Supplementary Table S7). Kaplan–Meier curves were plotted using the survival R package. Statistical significance was assessed using the log-rank test and defined as P < 0.05.

To assess if enrichment of resistant cluster signatures were associated with presurgical response to endocrine therapy, we obtained microarray gene expression data from the ACOSOG Z1031 (18, 46) trial (GSE136644), which randomized patients with stage 2 to 3 ER⁺/HER2⁻ breast cancer to preoperative aromatase inhibitor therapy with anastrozole, letrozole, or exemestane. No difference was seen in clinical efficacy between treatment arms and so treatment arms were grouped. Response was assessed using the authors’ validated endocrine response score (PEPI), in which a PEPI score of 0 indicates a favorable response. Signature scores were compared between samples with a PEPI score of 0 and those with a PEPI score greater than 0, and statistical significance was assessed using the Wilcoxon rank-sum test.

To determine enrichment of resistant cluster signatures in endocrine therapy–resistant ER⁺ tumors, we obtained bulk RNA-seq data from patients profiled on endocrine therapy (47). Serial tissue biopsies were taken at diagnosis (pretreatment), after 2 to 4 weeks on treatment, and on treatment as feasible (median biopsies, 3; range, 2–8). Clinical response was determined by RECIST 1.1 criteria and follow-up to annotate each specific timepoint as “sensitive” or “resistant.” Centroid values (median expression of signature genes) for each of the three resistant signatures were calculated for each sample. Signature scores were compared between sensitive and resistant samples using violin plots and statistical significance was determined using the Wilcoxon rank-sum test.

Raw data (10x FASTQs) and processed data for bulk and scRNA-seq data have been deposited at the NIH Database of Genotypes and Phenotypes (https://www.ncbi.nlm.nih.gov/gap/) and are available under the accession number phs003186.v1.p1.

All original code has been deposited on the Github and is publicly available at the Github repository BC_tamoxifen_response (https://github.com/hyunsoo77/BC_tamoxifen_response). Additional information regarding data analysis is available upon request ([email protected]).

Studying the responsiveness of ER⁺/luminal breast cancer cells to tamoxifen in vivo is challenging due to features of stromal cells and intratumoral heterogeneity. To address these challenges, we developed a novel ex vivo system for short-term culturing coupled with scRNA-seq. To demonstrate the platforms, two normal breast specimens were obtained from patients undergoing reduction mammaplasty and 10 patients with clinically ER⁺/HER2⁻ treatment-naïve primary breast tumors. Patient demographics and tumor characteristics are shown in Table 1 and are representative of our patient population (Supplementary Table S1). All tumors were strongly positive for ER and negative for HER2. Ki67 positivity ranged from 7% to 35% and eight tumors were luminal by bulk PAM50 subtype.

Two normal breast samples were obtained and processed (Fig. 1A). In total, 13,175 single cells were identified after quality control and annotated with canonical cell type markers (ref. 38; Fig. 1B and C). When comparing the treatment conditions, there was strong cluster overlap of control and tamoxifen-treated cells except in luminal progenitor cells, indicating a more robust change in gene expression with tamoxifen in that subpopulation.

Basal (KRT5/14) and luminal (KRT8/18) cytokeratins were enriched in the expected compartments (Fig. 1D) demonstrating robust classification of epithelial cell types. ESR1 transcripts encoding ER were not detected in most cells including luminal epithelial cells. FOXA1 was localized to mature and progenitor luminal epithelial cells, consistent with existing literature. Cell type abundance varied between the two specimens but not with tamoxifen treatment (Fig. 1E).

To evaluate the effect of time in suspension, we compared an immediate library from Normal_01 to the control suspension sample. UMAP plots were used to visualize the data by cell type (Supplementary Fig. S1A). Cell type was the primary driver of clustering rather than time in suspension (Supplementary Fig. S1B and S1C). Fewer cells were sequenced after suspension, but no enrichment or depletion of cell populations was observed (Supplementary Fig. S1D). Differentially regulated genes showed enrichment of inflammatory signaling pathways, epithelial–mesenchymal transition (EMT) genes, and genes overlapping with estrogen response (Supplementary Fig. S1E; Supplementary Table S2).

To determine tamoxifen-induced gene-expression changes in the distinct cell subpopulations that comprise normal human breast tissue, we performed differential gene expression analysis by cell type. Upregulated and downregulated genes with tamoxifen treatment were identified separately for basal epithelial cells (BEp), luminal progenitor cells (LEp_prog), and mature luminal cells (LEp) (Supplementary Table S3). Overlap analysis showed the majority of tamoxifen-regulated genes were cell type specific (Supplementary Fig. S2). Common downregulated genes included cytokeratins KRT8, KRT15, and KRT19, and hallmark estrogen-response genes including CCND1.

To determine differentially regulated pathways with tamoxifen treatment by cell type, we performed enrichment analyses (Fig. 1F; Supplementary Table S3). Luminal epithelial cells demonstrated strong depletion of E2-induced and enrichment of E2-repressed signatures, indicating capture of on-target effects of tamoxifen. Basal epithelial cells demonstrated downregulation of E2-induced genes, although less pronounced compared with luminal epithelial cells and unchanged expression of E2-repressed genes. Fibroblasts did not demonstrate signature enrichment with tamoxifen treatment. Endothelial cells had enrichment of TNFα signaling and induction of estrogen-repressed genes with tamoxifen, indicating a distinct response relative to other cell types. Cumulatively, these findings show that our novel system captures biological effects of tamoxifen treatment in epithelial cells and detects cell populations (i.e., fibroblasts) with differential cellular tamoxifen response. This demonstrates the ability to distinguish heterogeneity in drug response within diverse cell populations in primary human breast tissue.

Heterogeneity in gene expression and response to treatment occur in distinct subpopulations, even within cell line cultures (48). To elucidate heterogeneity in response to tamoxifen using our suspension-to-single-cell sequencing platform, we treated the ER⁺/HER2⁻ breast cancer cell line T47D to recapitulate our human tissue protocol (Fig. 2A). In bulk, T47D has a luminal B gene expression pattern (49). Cells grouped primarily into two subpopulations: groupA (luminal A–like) and a more proliferative groupB (luminal B–like) (Fig. 2B). Cell clustering occurred primarily by cell type rather than treatment condition (Fig. 2C).

Differentially expressed gene analysis identified 110 downregulated genes and 90 upregulated genes with tamoxifen treatment (Supplementary Table S4). The downregulated genes were enriched with canonical estrogen-induced genes and signatures of cell proliferation (Fig. 2D). The upregulated genes were enriched for gene signatures of androgen response and activation of the PI3K/mTOR signaling pathway. We used the 110 downregulated genes (adjusted P < 0.01), to establish a single-cell signature of tamoxifen response (scTAM-response-T47D). As expected, canonical luminal epithelial markers (KRT8, FOXA1, ESR1) were uniformly expressed, and the proliferation score (38) was enriched in the groupB subpopulation (Fig. 2E).

Ten tamoxifen/control pairs from primary breast tumors were processed using our novel experimental platform. In total, 40,428 cells passed quality control and were annotated using canonical lineage markers (Fig. 3A). Epithelial cells were further classified into malignant (Epi. Tumor) and nonmalignant (Epi. Nontumor) cells by estimating the CNV profiles using InferCNV (36).

Cells clustered by cell type identity and not by specimen, indicating high-quality sequencing. There was significant overlap of fibroblasts, endothelial, and immune cells across samples, but not epithelial cells (Fig. 3B). Malignant cells are visualized alone in Supplementary Fig. S3A. Despite strong clinical ER positivity for these tumors, ESR1 transcripts were zero in most tumor cells. This could be due to discordance between transcript abundance and protein levels, or more likely due to the depth of sequencing. The ER-associated transcription factor FOXA1 was enriched in tumor cells compared with nontumor epithelial and nonepithelial cell types (Fig. 3C), and the significant majority of tumor cells were assigned to the luminal PAM50 subtypes (Supplementary Fig. S3B), supporting true ER⁺ identity for most cells.

Cell type abundance and the number of sequenced cells was highly variable across samples (Supplementary Fig. S3B). The variability in cell type abundance between samples supports the use of single-cell assays to account for distinct cell subpopulations and mitigate effects on downstream analyses. In most specimens, cell type abundance was conserved between control and tamoxifen treatment within samples, with three specimens (Tumor_05, Tumor_06, and Tumor_09) showing increased ratio of luminal B/luminal A cells with tamoxifen treatment.

We leveraged the ability to computationally isolate tumor cells using InferCNV to assess changes in three canonical early estrogen response genes (CCND1, PGR, and GREB1; Fig. 3D) in only malignant cells within tumors. Of note, the tamoxifen dose of 10 μmol/L is higher than circulating levels in patients treated with tamoxifen, and so treatment conditions may be driving supraphysiologic cellular responses. Of the 10 tumor pairs, six showed a significant reduction in CCND1 expression with tamoxifen treatment, although one (Tumor_08) showed increased CCND1 expression. Five tumors had a significant reduction in PGR (encoding the progesterone receptor), and four tumors had reduced GREB1 expression. Overall, this shows reduced expression of estrogen-response genes after tamoxifen treatment in most of these 10 tumors. We did not observe a reduction in Ki67 expression in any tumor, which may be due to insufficient treatment time or low abundance of Ki67 transcripts.

To determine tumor-specific tamoxifen response genes and to assess response genes between tumors, we compared differentially expressed genes in malignant cells (i.e., InferCNV+; Supplementary Fig. S3C). Upregulated and downregulated genes were highly variable between samples. In general, samples with more downregulated genes had more upregulated genes, which could be due to intrinsic responsiveness or the number of malignant cells sequenced and the power to detect gene expression differences.

To assess variability in tumor-specific biological processes regulated by tamoxifen between patients, we assessed enrichment of gene sets. With tamoxifen treatment, estrogen-induced genes (16) were depleted in three tumors (Tumor_03, Tumor_05, and Tumor_09) but enriched in four tumors (Tumor_01, Tumor_02, Tumor_06, and Tumor_08; Fig. 3E). The upregulated genes in response to tamoxifen were enriched in gene sets related to the androgen response, p53, mTORC1 signaling, hypoxia, and apoptosis. MYC target genes and GATA3-induced genes were depleted in several tumors with tamoxifen treatment. The gene set for TNFα signaling was depleted in Tumor_05 and enriched in Tumor_06, highlighting the variability in tamoxifen-mediated inflammatory signaling. To account for variability in sequencing between bulk platforms and our single-cell system, we assessed the previously defined scTAM-response-T47D signature, which was depleted in four tumors. scTAM-response-T47D genes were enriched in Tumor_09 which was classified as a HER2-enriched bulk PAM50 subtype (Fig. 3E).

To determine heterogeneity in cellular tamoxifen response within specimens, we focused on the four tumors with depleted scTAM-response-T47D signature. Cells were clustered to visualize cell populations (Fig. 4A and B). In InferCNV+ malignant cells, top tamoxifen-regulated genes were plotted (Fig. 4C; Supplementary Table S5). Tamoxifen treatment was associated with depletion of proliferation signatures and estrogen- and GATA3-induced signatures (Fig. 4D), including ER-regulated genes (50) KRT7, CCND1, and KRT8. Downregulation of MYC targets, including FKBP4 and RANBP1 which alter ER-regulated genes, was observed with tamoxifen treatment (48, 51, 52).

Enrichment of several notable signatures was also observed among malignant cells within these four tumors with tamoxifen treatment (Fig. 4D; Supplementary Table S5). Signatures of EGFR and RAS (53) signaling were enriched including MAPK signaling genes (MAPK6, VEGFA, DUSP1, and DUSP10). The RAS–RAF–MAPK signaling pathway is known to promote resistance to endocrine therapy (8, 47, 54, 55). PI3K/mTORC1 signaling (53), a targetable pathway in patients with endocrine-resistant tumors (56), was also enriched in tamoxifen-treated malignant cells. Genes associated with EMT (57) and FOS/JUN signaling were also enriched (Supplementary Table S5) and are associated with worse outcomes and poor response to endocrine therapy in breast cancer. Cumulatively, enriched pathways show rapid adaptive upregulation of genes that may be driving a complex early transcriptional program to promote resistance to endocrine therapy.

To determine variability in tamoxifen response by cell type across specimens, tamoxifen-regulated genes were identified in nonmalignant cell populations (Supplementary Table S6). In contrast to malignant cells, fewer differentially regulated genes were seen in fibroblasts and macrophages. Similar to malignant cells, RAS-associated genes were enriched, but this was primarily due to mitochondrial genes. Tamoxifen downregulated gene sets in malignant cells and fibroblasts did not overlap. Macrophages demonstrated downregulation of hallmark inflammatory response with tamoxifen. Overlap analysis showed that the significant majority of both up- and downregulated genes were unique to specific cell types (Supplementary Fig. S4A and S4B).

These results demonstrate that distinct cell types within primary ER⁺ breast tumors have a unique response to tamoxifen. To test if a single-cell signature of tamoxifen resistance (upregulated genes with tamoxifen treatment) derived from only malignant cells (scTAM-resistance-M; Supplementary Table S5) was prognostic in ER⁺/HER2⁻ breast cancer patients treated with endocrine therapy, survival annotated transcriptional data were obtained from METABRIC (42) and SCAN-B (44). Patients were ranked by median expression of the signature genes and stratified into “high” or “low” groups and assessed with Kaplan–Meier curves. In both datasets, overall survival was significantly worse with high expression of the malignant cell-specific tamoxifen response signature (METABRIC: HR, 1.63; P = 0.023; SCAN-B: HR, 2.94; P = 0.002; Fig. 4E).

To measure heterogeneity in tamoxifen response in individual cells, we developed a response score for individual cells using the single-cell platform-derived scTAM-response-T47D signature. Each cell from the four targeted analysis tumors was assigned a score as the difference between the treated tumor cell's centroid signature value and the median centroid value for the matched, untreated tumor cells from the same patient. This metric measures the deviation from the starting state with tamoxifen treatment and can be quantified for each cell. Differential response to tamoxifen was thus identified between clusters (Fig. 5A). Cluster 2, composed primarily of LumA-like cells from Tumor_06, displayed a significant reduction in scTAM-response-T47D score, indicating cellular transcriptional response to treatment. In contrast, clusters 3, 12, and 19 demonstrated increased scTAM-response-T47D scores, indicating resistance to tamoxifen treatment at the single-cell transcriptional level. Cluster 19 was the only cluster composed of luminal B cells. In total, 13% of Tumor_03 was composed of resistant cells, 6% of Tumor_05, 2% of Tumor_06, and Tumor_08 was the only tumor in which the majority of cells (74%) were from resistant subpopulations (Fig. 5B).

To determine features of these transcriptionally unresponsive subpopulations, we identified differentially expressed genes for each of the three clusters compared with the responsive cluster 2 (Supplementary Fig. S5A; Supplementary Table S7). Clusters 3 and 12 were both characterized by enrichment of RAS signaling pathways and EMT, which can reduce estrogen responsiveness. Cluster 19 was characterized by high Ki67 and enrichment of proliferative signatures (refs. 40, 58, 59; Supplementary Fig. S5B). Cell-cycle phase analysis demonstrated that cluster 2 (sensitive) had the highest percentage of nonproliferative (GO/G1) cells at 44.5% compared with cluster 3 (35%), cluster 12 (30%), and cluster 19 (0%). Highly proliferative cells could indicate a population that was dividing or committed to dividing during the short-term treatment, thus completing the cell cycle under treatment rather than being truly “unresponsive.” Several reports have demonstrated that high Ki67 is associated with worse outcome in endocrine therapy-treated patients and that a reduction in Ki67 is a key clinical metric of response to endocrine therapy (18, 19, 60). Interestingly, all three clusters had increased expression of TACSTD2 encoding TROP2, the target of the antibody–drug conjugate sacituzumab govitecan, which has proven efficacy in metastatic ER⁺/HER2⁻ breast cancer (61). Therefore, all three clusters identified as resistant at the cellular level had gene expression patterns of known mechanisms of resistance to endocrine therapy and have increased expression of a druggable target that could be used to deliver precision, subpopulation specific therapy. These findings support the validity of this model to identify differential response within tumor cell subpopulations and identify potential treatment strategies.

To assess whether features of transcriptionally unresponsive subpopulations are associated with poor outcomes, we analyzed ER⁺/HER2⁻ patients treated with endocrine therapy in METABRIC and SCAN-B. Resistant cluster signatures were generated for each of the three clusters (Supplementary Table S7). Tumors were ranked by median and tertile cutoffs and survival assessed by Kaplan–Meier curves. The scTAM-resistance-C19 signature demonstrated the most robust stratification for overall survival (Fig. 5D), with high score corresponding to reduced survival in METABRIC (HR, 2.2; P < 0.001) and in SCAN-B (HR, 1.9; P = 0.002). High expression of the scTAM-resistance-C3 signature was significantly associated with reduced overall survival in SCAN-B (HR, 2.2; P = 0.03) but not METABRIC. Similarly, high expression of the scTAM-resistance-C12 signature was associated with reduced survival only in SCAN-B (HR, 1.9; P = 0.003; Supplementary Fig. S6A).

To determine whether the enrichment of tamoxifen-resistant cluster signatures is associated with clinical response in patients treated with endocrine therapy, we obtained two presurgical or unresectable datasets. The ACOSOG Z1031 (18, 46) trial treated patients with stage 2 to 3 ER⁺/HER2⁻ breast cancer with preoperative aromatase inhibitor therapy and evaluated response by their validated PEPI score. Patients with a PEPI score of 0 (good response) had lower pretreatment expression of scTAM-resistance-C19 signature compared with a PEPI score greater than 0 (−0.15 vs. 0.01; P = 0.043). No enrichment of signatures for clusters 3 and 12 based on PEPI score (Supplementary Fig. S6B). Finally, we obtained bulk RNA-seq data from patients with ER⁺/HER2⁻ breast cancer profiled on treatment with endocrine therapy (47). Clinically resistant tumors had significantly higher expression for all three resistant signatures relative to sensitive tumors (−0.11 vs. 0.16 for cluster 3; P = 0.00037; −0.12 vs. 0.21 for cluster 12; P = 0.00081; −0.12 vs. 0.26 for cluster 19; P = 0.00031; Fig. 5E).

Assessing heterogeneity in response to therapy within primary human solid organ tumors represents a significant challenge. For the first time, we report the use of single-cell transcriptional profiling coupled to a primary human tumor ex vivo experimental platform to dissect drug response within and across tumors. Using this platform, this study identified and characterized distinct responses to tamoxifen within normal human breast tissue and primary ER⁺/HER2⁻ human breast tumors. The ability to successfully distinguish differential response has several important implications. First, it allows for computational analysis of distinct cell populations, such as malignant cells, which can be used to derive more precise signatures of transcriptional response. This could be especially important in lobular breast cancers and other solid organ tumors with variable malignant and stromal compositions. Herein, we generated signatures from computationally identified tumor cells and showed that these were associated with prognosis and response to therapy across multiple large datasets.

We have additionally shown the ability to assess differential sensitivity to therapy by measuring gene-expression changes within malignant cells of the same tumor. In some tumors, resistance to therapy may be driven by low-abundance cell populations that are masked in bulk sequencing. Importantly, we identified three low-abundance cell populations that are unresponsive to tamoxifen. These resistant cells were identified within tumors that would have been classified as responsive (depleted E2-induced genes) at the tumor level. Thus, the ability to distinguish these cell populations may unmask hidden resistant cell populations. Supporting the accuracy of this characterization, all three populations had features of described mechanisms of resistance to endocrine therapy and were enriched in endocrine resistant tumors. Interestingly, these subpopulations all had increased expression of TROP2, the target of the antibody–drug conjugate sacituzumab govitecan. The presence of these low-abundance tamoxifen-resistant cells could identify patients in whom a precision medicine dual treatment strategy targeting specific tumor subpopulations is more efficacious than endocrine therapy alone. Further study is needed to develop clinically actionable strategies based on precise characterization of tumor cell subpopulation vulnerabilities.

Studying drug response on pre- and posttreatment specimens at the single-cell level from patients may be more directly relevant as in vivo cellular response may not mirror ex vivo response. Processing tumor cells may alter response, and in vitro doses may not mirror circulating and intratumor levels. An ideal validation experiment would obtain a pretreatment biopsy and assay using our platform at a range of tamoxifen doses. The patient would also undergo treatment with tamoxifen and single-cell libraries generated from pre-/posttreatment specimens. This would allow direct comparison of patient response, ex vivo response, and optimization of ex vivo doses to reproduce a representative response. There are several advantages to an ex vivo model. First, experiments in patients are more expensive and take longer to generate results. Next, our ex vivo primary tumor model allows for experiments that may be impractical or unethical in humans such as genetic perturbations and investigational drugs. Finally, this model could allow rapid parallel testing of multiple therapeutics to assess differential sensitivities and empirically optimize treatment strategies.

We recognize that our study has several limitations. The resistance analysis was done on cells from four patients, which is a small sample size with limited power for subgroup analysis based on menopausal status or ductal/lobular histology. Further, in vitro doses of 4-OH tamoxifen and estradiol used in this study are higher than circulating levels in patients, which could alter our measurement of response and could mask subtle transcriptional differences within tumor cell subpopulations. These doses may be especially relevant in postmenopausal women from whom most of our tumor samples were obtained. In contrast, normal samples were obtained from women in their 20s, which may limit the ability to compare tumor and normal samples in our study. In addition, the analysis of resistant cell populations was based on cells from tumors that were determined to be sensitive to tamoxifen as measured by depletion of a platform-matched T47D derived signature, which may limit the ability to characterize resistant populations. Signatures of resistant cells were assessed using existing clinically annotated bulk sequencing data, which indicates relevance in patients but does not directly validate these findings. Further study is needed to determine precisely how low-abundance resistant subpopulations are associated with risk of recurrence, such as longitudinal follow-up to determine if the presence of resistant subpopulations corresponds with the risk of recurrence. Ideally, this would be done with analysis of recurrent tissue samples to directly assess if recurrences originate from resistant subpopulations. Nevertheless, this work represents the first model system to annotate distinct treatment responses within individual cells of primary human breast tumors and offer a treatment target for resistant cells.

In conclusion, we present a novel ex vivo drug treatment platform coupled to scRNA-seq for human breast tumors. Through this platform, we identified distinctions in cellular response of subpopulations in normal human breast tissue and primary human tumors. Several mechanisms of resistance to tamoxifen were identified in resistant malignant cell subpopulations within human tumors, and signatures of those subpopulations predicted poor outcomes in patients with ER⁺/HER2⁻ breast tumors. Further studies are needed to determine clinically actionable strategies to identify patients at risk for therapeutic failure based on the presence of low-abundance drug-resistant subpopulations and determine cotreatment strategies to overcome resistance.

B.C. Calhoun reports grants from NCI during the conduct of the study and personal fees from PathAI outside the submitted work. C.M. Perou reports personal fees from Bioclassifier LLC outside the submitted work and has a patent for U.S. Patent No. 12,995,459 licensed and with royalties paid from Bioclassifier. P.M. Spanheimer reports grants from NIH/NCI and grants from Society of University Surgeons during the conduct of the study. No disclosures were reported by the other authors.

H. Kim: Data curation, software, formal analysis, visualization, methodology, writing–original draft. A.A. Whitman: Investigation, methodology, writing–original draft, writing–review and editing. K. Wisniewska: Investigation, methodology, writing–review and editing. R.T. Kakati: Formal analysis, investigation, writing–review and editing. S. Garcia-Recio: Investigation, writing–review and editing. B.C. Calhoun: Data curation, writing–review and editing. H.L. Franco: Resources, project administration, writing–review and editing. C.M. Perou: Conceptualization, resources, writing–review and editing. P.M. Spanheimer: Conceptualization, resources, supervision, funding acquisition, writing–original draft, writing–review and editing.

Dr. P.M. Spanheimer is supported by the NIH grants K08CA280388 and P50CA058223 and a Junior Faculty Award from the Society of University Surgeons Foundation. This work was supported in part by P30CA016086 UNC Lineberger Comprehensive Cancer Center Core Support Grant and a grant from the NIH/NCI (R01CA273444) to H.L. Franco.

We would like to thank the patients who donated tissue to make this work possible and Dr. Lisa Carey and Adrianna Warner for facilitating sample acquisition. This work was presented, in part, as a Spotlight Poster Discussion at the 2022 San Antonio Breast Cancer Symposium. We thank Albert Wielgus in the Pathology Services Core for expert technical assistance with histopathology/digital pathology including tissue sectioning, IHC staining, and imaging. The PSC is supported in part by NCI Center Core Support grant P30CA016086.

The publication costs of this article were defrayed in part by the payment of publication fees. Therefore, and solely to indicate this fact, this article is hereby marked “advertisement” in accordance with 18 USC section 1734.