Serial single-cell genomics reveals convergent subclonal evolution of resistance as early-stage breast cancer patients progress on endocrine plus CDK4/6 therapy

Abstract

Combining cyclin-dependent kinase (CDK) inhibitors with endocrine therapy improves outcomes for metastatic estrogen receptor positive (ER+) breast cancer patients but its value in earlier stage patients is unclear. We examined evolutionary trajectories of early-stage breast cancer tumors, using single cell RNA sequencing (scRNAseq) of serial biopsies from the FELINE clinical trial (#NCT02712723) of endocrine therapy (letrozole) alone or combined with the CDK inhibitor ribociclib. Despite differences in subclonal diversity evolution across patients and treatments, common resistance phenotypes emerged. Resistant tumors treated with combination therapy showed accelerated loss of estrogen signaling with convergent up-regulation of JNK signaling through growth factor receptors. In contrast, cancer cells maintaining estrogen signaling during mono- or combination therapy showed potentiation of CDK4/6 activation and ERK upregulation through ERBB4 signaling. These results indicate that combination therapy in early-stage ER+ breast cancer leads to emergence of resistance through a shift from estrogen to alternative growth signal-mediated proliferation.

Extended Data Fig. 1. Classification of patient tumors as sensitive or resistant to treatment, reflecting changes in tumor size observed at pathology relative to baseline.

Reconstructed trajectories of tumor burden are consistent with results of RECIST 1.1 MRI assessment at day 90 and allow sensitive and resistant tumors to be distinguished at end of treatment (day 180). a, Changes in tumor size during therapy for tumors classified as sensitive or resistant. Tumor growth (y-axis) calculated directly from data as the proportion tumor remaining at end of trial (final observed tumor size at pathology/baseline MRI tumor measurement). Values <1 indicate tumor shrinkage, whilst values>1 indicate an increase in size (Dashed horizontal line = no change in size during trial). A detailed biological response classification was determined by classifying tumors with similar trajectories using a Gaussian mixture model (colors). Sustained or partial responses were grouped and defined as sensitive tumors, whilst those with stable, progressive or rebound disease were classified as resistant tumors. The changes in tumor size are highly significantly different between resistance categories (two-sided ANOVA test: t=4.45, p<0.001). Violins show the distinct distribution of tumor growth observed across patients. Heatmap shows the strong agreement in the end of treatment classification obtained by classifying trajectories of tumor growth vs simple pathology/baseline MRI RECIST assessment of change in size during trial. Number patients (P) with sensitive (S) versus resistant (R) tumors by arm = Letrozole alone:P=11,(S=6, R=5); Intermittent high dose ribociclib: P=12 (S=6, R=6) Continuous low dose ribociclib: P=11 (S=4, R=7). b, Spiderplots show the reconstructed trajectories of tumor size (relative to day 0) during the trial, as inferred using all available clinical measurements of patients’ tumor size. Predicted tumor sizes at day 90 match the RECIST assessments of tumor response (top panels) whilst trajectories of tumor burden distinguish sensitive (shrinking) and resistant (persistent) tumor through to the end of the trial (bottom panels). Number patients (P) with sensitive (S) versus resistant (R) tumors by arm = Letrozole alone:P=11,(S=6, R=5); Intermittent high dose ribociclib: P=12 (S=6, R=6) Continuous low dose ribociclib: P=11 (S=4, R=7). c, Inferred change in tumor size between the start- midpoint (left panel) or start-end (right panel) of the trial, in patient response groups classified by either RECIST assessment at trial midpoint (top row) or the biological response classification from tumor trajectories (bottom row). RECIST assessments distinguish response/non-response at day 90 but not day 180, whilst the biological response classification does distinguish resistance or sensitivity at day 180 (two-sided ANOVA test: MRI day 180 p-value= 0.38 and Biological response day 90 p-value= 0.34). Number patients (P) with sensitive (S) versus resistant (R) tumors by arm = Letrozole alone:P=11,(S=6, R=5); Intermittent high dose ribociclib: P=12 (S=6, R=6) Continuous low dose ribociclib: P=11 (S=4, R=7).

Extended Data Fig. 2. Landscape of tumor and microenvironment of 10 patients with single nucleus isolated by ICELL8 platform.

a, t-SNE plot of 3,484 cells. Cells were classified into cancer cells, normal epithelial cells, immune cells, stromal cells, and unclassified cells, which are indicated by colors and labels. The 3,484 cells are from 7 patients (3 from the Intermittent high dose arm and 4 from the Continuous low dose ribociclib arm. b, Gene copy number profile in cancer cells and neighboring normal cells. Blue color indicates copy number loss and red color indicates copy number gain. c, Expression of marker genes of cancer cells and normal epithelial cells (KRT19, CDH1), stromal cells (FAP, HTRA1), and immune cells (PTPRC). d, Proportion of cancer cells and neighboring normal cells in each patient.

Extended Data Fig. 3. Mutational signature in 24 patients with whole-exome sequencing data.

a, Relative contribution of trinucleotide changes to three de novo mutational signatures identified in 24 patients. b, Relative contribution of each mutational signature to mutations in each patient.

Extended Data Fig. 4. Mutated genes in three frequently altered oncogenic pathways.

Genes are grouped by oncogenic pathway. Presence of gene mutations in each patient is colored as indicated in the legend. Treatment arm and clinical response (Response: sensitive, resistant) are indicated in final two rows of the plot (colors indicated in legend).

Extended Data Fig. 5. Intrinsic subtype of 35 patients with single nucleus isolated by 10x genomics platform and reduced subclonal estrogen receptor (ESR1) expression at end of therapy as correlated with increased basal-like pathway and Creighton endocrine therapy resistance signatures, independent of treatment.

a, Intrinsic subtyping. Each row represents a patient and each column represents an intrinsic subtype at three timepoints. The proportion of cancer cells in each intrinsic subtype was indicated by colors ranging from 0 to 85. Patient samples without cancer cells were indicated by gray. b, Reduced subclonal estrogen receptor (ESR1) expression. Top row shows the ESR1 expression and basal-like (left) and endocrine resistance (right) pathway signatures across subclonal cancer populations with differing MAPK activation (points) and the coloration signifies the treatment received. Fitted lines show the overall trend between ESR1 expression and pathway activity (shaded regions show 95% confidence bands). Bottom row shows the correlation between ESR1 expression and basal-like (left) and endocrine resistance (right) pathway signatures for each cancer subclone present at end of trial, in patients treated with different therapies (colors). Black points and error bars signifies the mean and confidence interval for the correlation between ESR1 and pathway activity under each treatment. Number of cells (n) and patients (P) with sensitive (S) versus resistant (R) tumors by arm = Letrozole alone: n= 46986, P=11, S=6, R=5; Intermittent high dose ribociclib: n= 27790, P=12, S=6, R=6; Continuous low dose ribociclib: n= 34543, P=11, S=4, R=7.

Extended Data Fig. 6. Divergence of JNK and ERK signalling pathway activity during treatment with combination therapy, especially in resistant tumors and heatmaps of the correlation between MAPK gene expression in each treatment arm (columns), showing the dichotomy between JNK and ERK activating genes across treatments.

A, JNK and ERK expression (color=pathway) during treatment (columns) in sensitive and resistant tumors (rows). Pathway trends determined across patients using hierarchical regression (solid lines). Inter-patient variability in pathway activity shown by dashed lines indicating patient specific responses and shaded regions showing confidence intervals of model estimates (JNK ssGSEA pathway=St JNK MAPK and ERK pathway=Biocarta ERK). Number of cells (n) and patients (P) with sensitive (S) versus resistant (R) tumors by arm = Letrozole alone: n= 46986, P=11, S=6, R=5; Intermittent high dose ribociclib: n= 27790, P=12, S=6, R=6; Continuous low dose ribociclib: n= 34543, P=11, S=4, R=7. B, Dendrograms show the collinearity of MAPK gene expression following each endocrine or combination therapies (columns).

Extended Data Fig. 7. Construction of the overall JNK activation phenotype score, utilizing this collinearity of gene expression between ERK and JNK genes

a, UMAP dimension reduction of MAPK genes, showing the bivariate Gaussian distribution of UMAP values, centered around the major axis of phenotypic variation (black line). The frequency of cells found in different parts of the UMAP phenotype space is shown by the color gradient. The major axis of phenotypic variation (the JNK activation phenotype) is identified as the first principle component in the UMAP phenotype space. b, Relationship between the JNK activation phenotype and expression of MAPK genes that are known a JNK activators (red) or ERK activators (blue) across subclonal cancer populations. Loess smooths are added showing the positive relationship between the JNK phenotype score and key JNK activators and the negative association between ERK activators and the JNK phenotype. Number of cells (n) and patients (P) with sensitive (S) versus resistant (R) tumors by arm = Letrozole alone: n= 46986, P=11, S=6, R=5; Intermittent high dose ribociclib: n= 27790, P=12, S=6, R=6; Continuous low dose ribociclib: n= 34543, P=11, S=4, R=7.

Extended Data Fig. 8. Correlation of growth factor receptors expression with estrogen pathway activity (Hallmark estrogen response early) in cancer cells from sensitive and resistant tumors under each therapy.

Strong negative correlations identify genes that are upregulated as estrogen signaling is lost. Specifically, tumors resistant to intermittent high dose and continuous low dose show compensatory activation of FGFR2 and ERBB4 respectively. Number of cells (n) and patients (P) with sensitive (S) versus resistant (R) tumors by arm = Letrozole alone: n= 46986, P=11, S=6, R=5; Intermittent high dose ribociclib: n= 27790, P=12, S=6, R=6; Continuous low dose ribociclib: n= 34543, P=11, S=4, R=7.

Extended Data Fig. 9. Transcriptional heterogeneity of key resistant genes.

a, sensitive and b, resistant tumors. For each patient’s tumor cells, a single-cell phylogenetic tree is shown at the center of circos plot. Cell annotation (timepoint and subclone) as well as expression of key resistance genes (ESR1, CDK6, FGFR2, ERBB4, RORA) are shown as heatmap. Phylogenetic tree of cells were constructed based on the distance between cell gene copy number profile. Subclones were inferred based on gene copy number profile. Zinbwave normalized gene expression were centered and scaled.

Extended Data Fig. 10. Reconstruction of cell cycle, fluctuations in gene expression during the cell cycle, distinct cell cycle phases, frequencies of cells throughout the cell cycle and shifts in gene expression within the cell cycle during therapy.

a, Single cell RNA seq gene expression profiles of cell cycle genes are extracted and used to perform dimension reduction with the UMAP algorithm. Cell cycle states (colors) with differing expression were identified using a Gaussian mixture model and the transitions between these states determined by the shortest distance to travel through each state and return to the original (Traveling salesman route=black line). b, Cells states ordered along the traveling salesman route. c, Example of fluctuations in gene expression of cells around the cell cycle (distance of points from origin=RB1 expression; colors=cell cycle state) Reconstruction of the fluctuation in average gene expression is predicted using a cyclical generalised additive model (black line with shaded confidence bands). d, Reconstructed fluctuations (coloured curves) in expression of genes around the cell cycle are used to classify distinct phases of the cell cycle (annotated by arrows around). Here we show four examples of key cell cycle genes, which influence the classification of cell cycle phases (G0, G1, S/G2). e, The frequency of cells in each stage of the cell cycle (height of bars) was counted and used to examine changes in the fraction of sampled cells in each phases cell cycle phase over time and between treatment and response groups. f, During treatment, the changes in gene expression fluctuations around the cell cycle were examined. Distance of the curve from the origin indicates gene expression and colored curves shows expression at different timepoints. g, Consistent cell cycle stages present across patients. For each patient (subpanel), single cell RNAseq gene expression profiles for cell cycle genes were extracted and the fitted UMAP model used to project cells onto the lower dimensional cell phenotype space (UMAP dimensions 1 and 2). Cell cycle stages (colors) with differing expression, identified using the Gaussian mixture model, were overlaid, showing that all patients have cells that are distributed across the cell cycle phenotype space. The traveling salesman route (black line) shows the transitions between these stages, as determined by the shortest distance to travel through each state and return to the original.

Figure 1.. Landscape of tumor and macroenvironment of early-stage ER+ breast cancer patients in FELINE trial.

a, Schematic diagram of single-cell RNA-seq workflow. The data of 34 patients generated using the 10x genomics platform are shown in Fig.1 b–d. Number of cancer cells (n) and patients (P) with sensitive (S) versus resistant (R) tumors by arm = Letrozole alone: n= 46986, P=11, S=6, R=5; Intermittent high dose ribociclib: n= 27790, P=12, S=6, R=6; Continuous low dose ribociclib: n= 34543, P=11, S=4, R=7. The data of 10 patients generated using ICELL8 platform are shown in Extended Data Fig.2. b, Distinction of different cell types is shown by a t-SNE plot of single cells of 34 ER+ breast cancer patients (color=cell type). c, Gene copy number profile in cancer cells and neighboring normal cells. Blue color indicates copy number loss and red color indicates copy number gain. d, Expression of marker genes of cancer cells and normal epithelial cells (KRT19, CDH1), stromal cells (FAP, HTRA1), and immune cells (PTPRC). e, Proportion of cancer cells and neighboring normal cells in each patient.

Figure 2.. Evolution of genomic mutations in response to endocrine or combination therapy.

a, Frequently mutated genes. Cancer driver genes were ranked based on number of patients carrying large impact somatic mutations. Mutations were counted if they are present in one or both biopsies. The top 10 frequently mutated genes are shown. b, Copy number alteration of key resistant genes in sensitive and resistant tumors. Gene copy number alterations were counted for a patient if one or both biopsies carry the alteration. c, Clonal evolution in response to endocrine therapy or combination therapy. Patient tumors were ranked by initial diversity (Shannon index at Day 0). The height of fishplot scaled to reflect changes of tumor size during treatment. Two biopsies (pre- and post-treatment) were sequenced and used to perform analyses of clonal evolution for each tumor. Day 180 biopsies were used as post-treatment samples except P36, P44, P45, and P46. Only Day 14 biopsies were available for sequencing for these four patients. Copy number alteration (CNA) of ESR1 loss, TP53 loss, and AKT3 gain are labelled for each patient when present at one or both biopsies. A “T” in parentheses indicates CNA is truncal in this tumor defined based on FACETS copy number analysis using WES and inferCNV copy number profiles using scRNAseq. Number of cells (n) and patients (P) with sensitive (S) versus resistant (R) tumors by arm = Letrozole alone: n= 46986, P=11, S=6, R=5; Intermittent high dose ribociclib: n= 27790, P=12, S=6, R=6; Continuous low dose ribociclib: n= 34543, P=11, S=4, R=7. d, Change of tumor heterogeneity in response to endocrine therapy or combination therapy, as measured by Shannon index of overall diversity, and broken down in to the diversity components of subclonal dominance and richness. Violin curves show the between patient variability in the change of heterogeneity during the trial (small points show observed changes in diversity). Hierarchical regression models identified the average change in tumor heterogeneity (large points) during endocrine (blue) or combination therapy (yellow) (uncertainty quantified by 95% confidence interval error bars). Schematic diagrams show the distinction between differences in dominance and richness.

Figure 3.. Accelerated evolution of estrogen independence during combination therapy.

a, Changes in single cell estrogen receptor activity and signatures of estrogen dependent (luminal) and estrogen independent (basal) phenotypes during treatment (columns), in tumors resistant (persistent=red) or sensitive (shrinking=blue) to therapy. Pathway trends across tumors were determined using hierarchical regression (solid lines). Tumor’s pathway trajectories are shown by dashed lines and confidence intervals of model estimates shown by shaded regions. Pathway activity measured using the ssGSEA signatures (Hallmark estrogen response early, Smid breast cancer luminal A up and Smid breast cancer basal up respectively). Number of cells (n) and patients (P) with sensitive (S) versus resistant (R) tumors by arm = Letrozole alone: n= 46986, P=11, S=6, R=5; Intermittent high dose ribociclib: n= 27790, P=12, S=6, R=6; Continuous low dose ribociclib: n= 34543, P=11, S=4, R=7. b, Loss of heterozygosity (LOH) of the estrogen receptor gene (ESR1) associated with reduced average expression in pre-treatment biopsies (two-sided Mann-Whitney test). Violin curves show the distribution of average ESR1 expression in tumors with or without LOH. Points indicate patient specific averages (color indicates treatment and shape signifies tumor response). Box and whisker plots indicate the median and upper/lower quantiles of patient tumor expression, with whiskers signifying the data range. c, Reduction in ESR1 expression accelerated under combination treatments, compared to endocrine therapy (columns: letrozole vs ribociclib treatments). Violin curves show the distribution of single cell expression during each treatment (color). Expression normalized relative to the baseline average (grey dashed line). Hierarchical generalized additive models predict the changes in expression during treatment, accounting for initial patient specific difference (colored curves). Number of cells (n) and patients (P) with sensitive (S) versus resistant (R) tumors by arm = Letrozole alone: n= 46986, P=11, S=6, R=5; Intermittent high dose ribociclib: n= 27790, P=12, S=6, R=6; Continuous low dose ribociclib: n= 34543, P=11, S=4, R=7.

Figure 4.. JNK pathway activation occurs during the emergence of combination therapy resistance and is associated with estrogen independence and increased CDK6 expression.

a, JNK pathway activity increases during treatments (columns) in resistant tumors (red) compared to sensitive tumors (blue) at day 14, and in cancer cells of most tumors by end of treatment (day 180). Pathway trends determined using hierarchical regression (solid lines). Tumor’s pathway trajectories are shown (dashed lines) along with confidence intervals of model estimates (shaded regions) (JNK ssGSEA pathway=St JNK MAPK). Number of cells (n) and patients (P) with sensitive (S) versus resistant (R) tumors by arm = Letrozole alone: n= 46986, P=11, S=6, R=5; Intermittent high dose ribociclib: n= 27790, P=12, S=6, R=6; Continuous low dose ribociclib: n= 34543, P=11, S=4, R=7. b, Heatmap of correlation between MAPK gene expression, showing the dichotomy between JNK and ERK activating genes across treatments (ESR1, ERBB4 and FGFR2 receptors added to indicate their relationship to MAPK genes). Number of cells (n) and patients (P) with sensitive (S) versus resistant (R) tumors by arm = Letrozole alone: n= 46986, P=11, S=6, R=5; Intermittent high dose ribociclib: n= 27790, P=12, S=6, R=6; Continuous low dose ribociclib: n= 34543, P=11, S=4, R=7. c, End of trial expression of CDK6 (ribociclib target gene) in subclonal tumor populations with differing levels of JNK activation (phenotype integrates across MAPK gene expression) and for patients with CDK6 genetic amplification (triangles and blue curves) or normal copy number (circles and red curves). Generalized additive models describe the relationship between JNK signaling activity and CDK6 expression at end of therapy (curves), with shaded regions indicating model confidence bands. Average estrogen receptor (ESR1) expression is shown for each subclonal population with differing JNK activation (point color gradient). Number of cells (n) and patients (P) with sensitive (S) versus resistant (R) tumors by arm = Letrozole alone: n= 46986, P=11, S=6, R=5; Intermittent high dose ribociclib: n= 27790, P=12, S=6, R=6; Continuous low dose ribociclib: n= 34543, P=11, S=4, R=7.

Figure 5.. Activation of ERBB4 and FGFR2 as resistance mechanisms to endocrine and combination therapy.

a, Average end of trial growth factor receptor expression (point color gradient) of ERBB4 (top) and FGFR2 (bottom) is shown for subclonal populations with differing JNK activation (x-axis) and for tumors with genetic amplification (triangles and blue curves) or normal copy number (circles and red curves). Cyclin-dependent kinase 6 (CDK6; ribociclib target gene) expression is shown for each population and generalized additive models (curves) describe the relationship between JNK signaling activity and CDK6 expression at end of therapy, (shaded regions indicate model confidence bands). Number of cells (n) and patients (P) with sensitive (S) versus resistant (R) tumors by arm = Letrozole alone: n= 46986, P=11, S=6, R=5; Intermittent high dose ribociclib: n= 27790, P=12, S=6, R=6; Continuous low dose ribociclib: n= 34543, P=11, S=4, R=7. b, Hierarchical clustering of tumors showing dysregulated ESR1 as well as upregulated ERBB4 and/or FGFR2. For each patient’s tumor biopsies, a hierarchical clustering tree is shown at the center of circos plots. Cell annotation (timepoint and subclone) as well as expression of key resistance genes (ESR1, CDK6, FGFR2, ERBB4, RORA) were shown as heatmap. c, Schematic diagram showing resistance mechinisms driven by upregulation of ERBB4 and CDK6 amplification (red circle signifies amplification) or alternative signaling via FGFR2/RTK’s and JNK signal transduction.

Figure 6.. Cell cycle reactivates during combination therapy follows the loss of estrogen receptor expression, activation of JNK1, repression of the cell cycle inhibitor cyclin-dependent kinase 2A and upregulation of CDK6 during the G1 checkpoint phase.

a, Cell cycle activity of resistant (red) and sensitive (blue) tumors during treatment (columns= regimes). Trend in cell cycle activity, measured by the ssGSEA biocarta cell cycle pathway, are determined by hierarchical regression (solid lines). Tumor specific trajectories are shown (dashed lines) along with confidence intervals of model estimates (shaded regions). Number of cells (n) and patients (P) with sensitive (S) versus resistant (R) tumors by arm = Letrozole alone: n= 46986, P=11, S=6, R=5; Intermittent high dose ribociclib: n= 27790, P=12, S=6, R=6; Continuous low dose ribociclib: n= 34543, P=11, S=4, R=7. b, Visualization of the pseudotime cell cycle reconstruction obtained using the Markov model-based reCAT algorithm. The dynamics of key cell cycle gene expression across stages of the cell cycle (black lines) were recovered from cell specific gene expression (points), using cyclical generalized additive models. The recovered fluctuations of cell cycle gene expression were used to identify three distinct cell cycle phases (colors; G0, G1 and S/G2), using a Gaussian mixture model. Cell cycle stages (clusters of cells) are colored by cell cycle phase and the distance of the point from the origin signifies the cells expression in that stage. Cell cycle orientation is consistent and comparable across circular plots. c, The proportion of S/G2 phase cells (passed the G1 checkpoint) in samples from resistant and sensitive tumors (color) during each treatment (column). Logistic generalized additive models describe trends in S/G2 phase cell frequency over time (solid lines) and heterogeneity across tumors (shaded regions signify confidence bands). Number of cells (n) and patients (P) with sensitive (S) versus resistant (R) tumors by arm = Letrozole alone: n= 46986, P=11, S=6, R=5; Intermittent high dose ribociclib: n= 27790, P=12, S=6, R=6; Continuous low dose ribociclib: n= 34543, P=11, S=4, R=7. d, Changes in ESR1, JNK1, CDKN2A expression around the cell cycle and during treatment (columns). Colored lines show the expected gene expression of cells throughout the cell cycle prior (blue), during (orange) and after treatment (red). The distance from the center of the circle shows gene expression at point in the cell cycle.