Endocrine therapy resistance (ETR) remains the greatest challenge in treating patients with hormone receptor–positive breast cancer. We set out to identify molecular mechanisms underlying ETR through in-depth genomic analysis of breast tumors.
We collected pre-treatment and sequential on-treatment tumor samples from 35 patients with estrogen receptor–positive breast cancer treated with neoadjuvant then adjuvant endocrine therapy; 3 had intrinsic resistance, 19 acquired resistance, and 13 remained sensitive. Response was determined by changes in tumor volume neoadjuvantly and by monitoring for adjuvant recurrence. Twelve patients received two or more lines of endocrine therapy, with subsequent treatment lines being initiated at the time of development of resistance to the previous endocrine therapy. DNA whole-exome sequencing and RNA sequencing were performed on all samples, totalling 169 unique specimens. DNA mutations, copy-number alterations, and gene expression data were analyzed through unsupervised and supervised analyses to identify molecular features related to ETR.
Mutations enriched in ETR included ESR1 and GATA3. The known ESR1 D538G variant conferring ETR was identified, as was a rarer E380Q variant that confers endocrine hypersensitivity. Resistant tumors which acquired resistance had distinct gene expression profiles compared with paired sensitive tumors, showing elevated pathways including ER, HER2, GATA3, AKT, RAS, and p63 signaling. Integrated analysis in individual patients highlighted the diversity of ETR mechanisms.
The mechanisms underlying ETR are multiple and characterized by diverse changes in both somatic genetic and transcriptomic profiles; to overcome resistance will require an individualized approach utilizing genomic and genetic biomarkers and drugs tailored to each patient.
Many patients with estrogen receptor–positive breast cancer develop drug resistance during endocrine therapy treatment. The molecular mechanisms underlying endocrine therapy resistance are not fully understood. We sought to investigate this by utilizing comprehensive and integrated DNA and RNA analysis of samples of breast cancers from 35 patients treated over time with endocrine therapy, whose cancers were responsive, had intrinsic resistance, or developed acquired resistance. We identified known and novel mutations enriched in resistant tumors, as well as unique transcriptomic profiles specific to resistant tumors leading to increased activation of multiple oncogenic signaling pathways. Each patient displayed distinct molecular profiles, showing a variety of mechanisms of resistance, corroborating the highly heterogeneous nature of endocrine resistance. Improving the treatment of resistant tumors will require a personalized approach tailored to patient-specific causes of resistance.
Breast cancer is one of the leading causes of mortality in women. Over 75% of breast cancers are estrogen receptor–positive (ER+) and approximately 65% of these are also progesterone receptor–positive (PR+). ER+ and/or PR+ breast cancer is also known as hormone receptor–positive (HR+), and has the best prognosis among all breast cancers, with 5-year relative survival rates estimated to be around 90% (1). Good outcomes are largely attributed to success in therapeutically targeting HRs. However, not all HR+ patients respond to endocrine therapy due to a variety of reasons including molecular heterogeneity (15%–20%: intrinsic resistance) as well as cellular phenotypes that develop and evolve during treatment (30%–40%: acquired resistance; ref. 2). Breast cancers exhibit unique somatic mutations, copy-number alterations (CNA), and transcriptomic profiles, which are particularly frequent among the HR+ subgroup (3–5) and can contribute to endocrine resistance.
Commonly used endocrine agents in clinical practice that target the ER include the selective ER modulator tamoxifen and the selective ER degrader fulvestrant (6). An alternative option is to reduce levels of circulating estrogen using the aromatase inhibitors (AI: letrozole, anastrozole, and exemestane); these work by inhibiting aromatase, and thus blocking estrogen production in postmenopausal women. They can be used in premenopausal women when combined with ovarian suppression or ablation. Despite there being a variety of agents with different modes of action, resistance can develop to each of these treatments, and is inevitable in the metastatic setting.
Extensive studies have investigated resistance mechanisms to endocrine therapies (7–9). These fall mainly into the following categories: changes in the estrogen signaling pathway, activation of growth factor signaling pathways, and cell-cycle dysregulation. Changes in the estrogen signaling pathway represent the most common resistance mechanisms where genetic and/or epigenetic alterations in ER or ER-associated transcription factors and coactivators have been identified (10–13). Specifically, somatic mutations in ESR1 following AI treatment alter ER signaling and result in endocrine insensitivity due to constitutive activation (14). Growth factor signaling pathway activation including EGFR/HER2, FGFR, MAPK, and PI3K has also been associated with resistance (15–18). Cell cycle regulator changes including amplification of cyclins, cyclin-dependent kinases (CDK), MYC, and loss of RB can uncouple the cell cycle from estrogen-mediated entry into G1/S and represent a further resistance mechanism (19–21).
Breast cancer is a highly heterogeneous disease and the full range of mechanisms underlying endocrine resistance has not as yet been fully characterized. To best identify and characterize known and postulated resistance mechanisms, studies using patient-derived materials utilizing comprehensive DNA and RNA sequencing (RNA-seq) can provide optimal insight, however, there have been few comprehensive studies addressing this using a multi-omics approach. Here, we present an integrated analysis of both DNA sequencing and RNA-seq of a cohort of 35 patients with ER+ breast cancer treated with primary endocrine therapy, who had tumor biopsies performed before and during treatment with endocrine therapy, and at the time of recurrence in those with progressive disease. Clinical response was defined according to institutional standards by monitoring changes in tumor volume by 3D ultrasound in the neoadjuvant setting and conventional monitoring for recurrence in the adjuvant setting. Twelve patients received two or more lines of endocrine therapy agents, with second-line agents being initiated at the time of development of resistance to first-line agents.
The aim of this study was to delineate the mechanisms of endocrine resistance in these patients. We have examined mutation and copy number landscapes, and transcriptional profiles, and found genetic and transcriptomic changes specific to resistant tumors. These results highlight the genetic diversity of the resistance process.
Materials and Methods
Patients, consent, and tissue processing
Tumor tissue was obtained from 35 patients with ER+ breast cancer who consented to participate in this research study at the Edinburgh Breast Unit and University of Edinburgh (Edinburgh, UK). Written informed consent was obtained from all patients. Ethical approval for the study was granted under the Lothian NRS BioResource approval number 20/ES/0061 and the study was conducted in accordance with Scottish Common Law. Patient characteristics are presented in Table 1 and Supplementary Table S1. All patients received neoadjuvant or primary endocrine therapy, which was continued in the adjuvant setting for those patients who eventually underwent surgery. In total, 23 patients were treated with a single endocrine agent only, while 12 patients received two or more lines, with second-line agents being initiated at the time of resistance to first-line agents in those patients who did not wish surgery or were deemed unfit for surgery at that time. Neoadjuvant treatment in this cohort ranged from 3.6 months to 61.2 months with a median duration of 7.8 months. Of 34 patients who underwent surgery, 23 had a wide local excision and 11 had mastectomies. Six patients had surgery at the time of cancer recurrence while receiving adjuvant endocrine therapy with samples at surgery being obtained for analysis (Fig. 1).
|Age at diagnosis||75 y (mean); 43–95 y (range)|
|ER+ (Allred 7)||7|
|ER+ (Allred 8)||28|
|Lymph node status|
|Unknown (no surgery)||1|
|Wide local excision||23|
|Age at diagnosis||75 y (mean); 43–95 y (range)|
|ER+ (Allred 7)||7|
|ER+ (Allred 8)||28|
|Lymph node status|
|Unknown (no surgery)||1|
|Wide local excision||23|
Abbreviation: y, years.
Neoadjuvant clinical response was assessed by monitoring dynamic changes in 3D tumor volume, determined by repeat ultrasound measurements recorded during the neoadjuvant treatment window and performed by a single sonographer (J.M. Dixon). Clinical response classification was determined using RECIST 1.1 criteria. All patients had long-term follow-up and for those who had surgery (n = 34), adjuvant response was determined by monitoring for cancer recurrence. Of the 35 patients, 13 patients remained endocrine therapy responsive throughout treatment, 3 patients had intrinsic resistance to neoadjuvant therapy, and 19 patients acquired resistance: 13 during neoadjuvant therapy and 6 developed a recurrence on adjuvant endocrine therapy following previous neoadjuvant therapy. Of the 6 who developed a recurrence, 3 went on to develop further recurrences on subsequent adjuvant endocrine agents: each had surgery followed by a new adjuvant endocrine therapy agent at the time of each recurrence. Overall, 13 patients remained sensitive to the first endocrine agent (Supplementary Table S1). All patients received clinical follow-up for 10 years from diagnosis or until death. No responsive patients developed a recurrence in that time, and none died from breast cancer. For each patient, tumor core biopsies were taken during the neoadjuvant treatment window using ultrasound-guided 14-guage needle core biopsy. On average, three sequential tumor tissue samples were taken from each patient (range 2–8 samples). Biopsies were taken at diagnosis (pretreatment), at 2 to 4 weeks on-treatment, and obtained at progression or every 3 to 6 months in those patients who declined or were not deemed fit enough for surgery at the time and were instead managed by longer-term primary endocrine therapy. A relative increase in tumor volume of 20% or greater prompted surgery or rebiopsy and initiation of a second- or third-line endocrine agent. Tumor tissue was collected at surgery and when recurrence developed. All surgeries were performed to clear margins by a single surgeon (J.M. Dixon). Wherever feasible, multiple tumor samples were taken from each recurrence. In total, sufficient RNA and DNA was successfully extracted from 32 primary tumors and 137 tumor specimens during treatment and was used for DNA whole-exome sequencing (WES) and RNA-seq.
For 17 patients, fresh-frozen tissue was available. Core biopsy and surgical tissues were snap-frozen in liquid nitrogen immediately after the procedure. For the remaining 18 patients, formalin-fixed, paraffin-embedded (FFPE) tissues were analyzed. For fresh-frozen tissues, DNA and RNA were isolated using QIAGEN RNAeasy and DNAeasy kits, respectively, according to the manufacturer's protocols. For FFPE tissues, RNA and DNA were isolated with QIAGEN AllPrep FFPE Kit, according to the manufacturer's protocols. RNA quality was assessed with an Agilent BioAnalyzer RNA 6000 Nano Kit.
DNA was prepared for sequencing using the Agilent Technologies SureSelect XT library protocol. Fresh-frozen tumors were processed according to the manufacturer's protocol for 3 μg input, while FFPE tumors were processed with the low-volume input according to the manufacturer's protocol for 200 ng input. DNA libraries were captured and amplified with Agilent Technologies SureSelect Human All Exon, version 5 for fresh-frozen tissues or version 6 for FFPE tissues, according to the manufacturer's protocol. The quality of both the DNA libraries and DNA exome capture and concentration were quantified with Agilent TapeStation DNA 1000 and High Sensitivity D1000, respectively. Paired-end sequence data (2 × 100 bp) were generated using the Illumina HiSeq 2500 for each tumor or normal sample, with three samples per lane. Illumina reads were mapped to the hg19 reference sequence with Burrows-Wheeler Aligner (BWA) 0.7.9a (22), realigned with ABRA, version 0.96 (23), and processed by biobambam2 (24). Viral alignments were counted with Samtools (25) and BEDTools-Version- 2.15.0. Picard 1.92 (26) was used to calculate sequencing metrics. ISAAC (27) and Freebayes were used to call germline mutations with quality scores above 30. SnpSIFT, version 1.3.4, band SnpEFF (28) was used to annotate alterations with population-level frequencies. CADABRA SomaticLocusCaller was used for further filtration. Somatic variants were called with STRELKA (29) using strelka_config_bwa_default.ini.
Minor allele frequencies of highly variable SNPs called by Freebayes in the general population were used for sample identity. All tumor normal pairs had an expected 89% to 100% identity from the same patient.
DNA CNAs were identified with SynthEx (30) using 50,000-bp– sized bins and K-nearest neighbors (KNN) = 4 from the pool of 31 available normal tissues. Briefly, the ratio of on-target and off-target exome reads of the tumor were compared with a normal tissue selected from the dataset by the highest degree of similarity by KNN based on library size and fold enrichment. Segment-level ratios were calculated and log2-transformed. Copy number levels greater than 0.25 were considered gains, and levels below –0.2 were considered losses. Using Ensemble hg19 gene annotations, genes were mapped to segments in each sample and a gene by sample copy number matrix was constructed. Specifically, a gene that totally falls into a segment was assigned the copy number of that segment.
Fresh-frozen RNA ribo-zero libraries were prepared with rRNA removed from total RNA using Epicentre's Ribo-Zero rRNA Removal kit (catalog no. RZH11042). Thirty to 100 ng Ribo-Zero RNA was used for the construction of the library using the Illumina TruSeq RNA Sample Prep Kit (catalog no. RS-122–2001) and followed the manufacturer's instruction, except for omitting the purification step before fragmentation. FFPE RNA was prepared with the Illumina TruSeq FFPE RiboZero Gold protocol according to the manufacturer's instructions. RNA libraries were sequenced as 2 × 50 bp paired-end reads with two samples per lane on Illumina HiSeq 2500 sequencers. Reads were aligned with STAR2.4.2a (31), and gene values were quantitated with Salmon (32). Raw gene counts were upper quartile normalized, filtered to genes that present in over 70% of samples, and log2-transformed. To correct the batch effect between fresh frozen and FFPE samples, log2-transformed gene expression values for fresh-frozen and FFPE sample sets were separately median centered and column standardized, then merged to a uniform dataset for downstream analysis (33, 34).
PAM50 subtyping was implemented as described in Fernandez-Martinez and colleagues (35). Briefly, batch-corrected gene expression data were normalized to ER+ samples in the PAM50 training set (36). Correlation to each subtype centroids was calculated, and subtypes were called according to the nearest centroid.
Computational and statistical analysis
Gene expression signatures
A panel of 654 previously published gene expression signatures were used to fully characterize cancer expression phenotypes (Supplementary Table S2). These 654 signatures were obtained from multiple publications or gene set enrichment analysis (GSEA; ref. 37) and were summarized in Xia and colleagues (26). Signature scores were calculated in a way consistent to how they were derived, noting that the large majority were median expression of a predetermined set of genes (38).
Linear mixed model of CNA, gene expression, and gene signatures
To identify differential CNAs, gene expressions, and gene signatures, biological replicates in CNA data, batch corrected gene expression, data, or gene signatures data were first collapsed into pseudo-samples by taking the mean value of multiple biological replicates (for example, sample 020–4-1A, 020–4-1B, 020–4-1C, 020–4-1D, and 020–4-2A were collapsed into 020–4 by averaging each CNA/gene/gene signature values). Patients were then filtered to those with documented acquired resistance and those that had at least one sensitive tumor and at least one resistant tumor, resulting in N = 13 patients with 38 samples for CNA analysis and N = 17 with 53 samples for gene expression, and gene signature analysis. For gene expression, genes were further filtered to those with Entrez ID. Each CNA, gene, or signature was then tested for differential expression in the resistant tumors versus sensitive tumors, with the patient taken into account as a confounding variable using lme4 package in R (39): lmer[value ∼ sensitive/resistant + (1|patient)]. Permutation-based FDR was calculated by permuting the tumor sensitive/resistant labels 100 times.
Hierarchical clustering of gene expression, and gene signatures
For gene expression, hierarchical clustering was done using correlation distance metric and complete linkage. Clustering for gene signatures was done using Euclidean distance metric and complete linkage. All clustering was done using R package pheatmap.
Computational reinterrogation of somatic mutations in related tumors
Low read coverage or low tumor cell purity can cause the rigorous somatic mutation caller to miss mutations (23–40). Thus, all of the high-confidence somatic mutations from every tumor taken from 1 patient were reinterrogated within the same tumors from that same patient. First, all of the somatic mutations from the tumors within a patient were collapsed into one file, excluding any guanine-to-adenine or cytosine-to-thymine mutations from FFPE samples. For each mutation from a single patient, we then counted the mutant and reference alleles at that position from the original BAM file of each tumor from that patient. Variant allele frequencies (VAF; alternate counts/total read counts) were recalculated from the new calls. All mutations from the data set were interrogated in the normal sequence for all tumors in this data set to account for false-positives. Mutations with VAFs of greater than 20% in at least two normal tissues from unrelated patients were excluded from future analyses.
Count table for mutations
We filtered patients to those with acquired resistance and 12 patients had at least one sensitive and one resistant tumor sample with successful DNA sequencing. For each gene in the Pan-Cancer 299 significantly mutated genes, we divided the 12 patients into the following four categories according to mutation status of the gene: (i) patients with mutation present in both sensitive and resistant tumors, (ii) patients with mutation present in resistant tumors only, (iii) patients with mutation present in sensitive tumors only, and (iv) patients with mutation present in neither sensitive or resistant tumors. All sensitive tumors and all resistant tumors from each patient were considered as a whole.
Mutational signatures were calculated using R package DeconstructSigs v1.8.0 to determine the weights of known mutational processes (COSMIC mutational signatures v2 – March 2015: https://cancer.sanger.ac.uk/signatures/signatures_v2/) in each tumor sample. We then focused on four signatures that are known to be important in breast cancer including Signature 1: age-related; Signature 2: and Signature 13: apobec-related; and Signature 3: homology directed repair (HDR)-related.
All statistical analyses were performed using R v. 3.5.2 in RStudio.
DNA WES data are available at the EMBL-EBI ArrayExpress database under the accession number E-MTAB-10080. RNA-seq data are available at the EMBL-EBI ArrayExpress database under the accession number E-MTAB-9917.
We examined the clinical features and molecular subtypes of all 169 samples, representing 35 patients (Supplementary Table S3). This cohort had a median age of 75 years at the diagnosis of breast cancer and a median follow-up of 4 years (Table 1). Responses to endocrine therapy were categorized as intrinsic resistance (nonresponse, n = 3), acquired resistance (initial response followed by tumor regrowth, n = 19), or ongoing sensitivity (initial response that was maintained or no evidence of recurrence, n = 13; Fig. 1; Supplementary Fig. S1). Thirty-three patients received letrozole as initial neoadjuvant therapy, 1 received anastrozole, and 1 received 2 weeks of fulvestrant prior to surgery and subsequent adjuvant tamoxifen, as part of a clinical study. Of the 33 who received primary letrozole, 3 were changed to alternative endocrine therapies (2 received tamoxifen and 1 received anastrozole) due to side-effects with letrozole. All patients had ER+ tumors (Allred score 7–8; Table 1) by clinical pathologic report and 3 had HER2+ tumors by IHC 3+; Fig. 1). We applied the PAM50 subtype predictor to determine their intrinsic molecular subtype (35, 36). Of 32 primary tumors with successful RNA-seq, 14 were Luminal A (LumA); 14 were Luminal B (LumB), 3 were HER2-enriched, and 1 was Normal-like. Among the three HER2-enriched primary tumors, 2 were HER2+ by IHC and 1 was HER2− by IHC. Of note, the three intrinsically resistant primary tumors had subtypes of LumA, LumB, and Normal-like. This Normal-like primary tumor had very close correlation to LumA centroid and to the Normal-like centroid, which is not uncommon for many low cellularity LumA tumors (4). Of all 167 tumor specimens subtyped including the primary tumors, 79 were LumA, 50 were LumB, 17 were HER2-enriched, 1 was Basal-like, and 20 were Normal-like. Only 8 patients had consistent subtypes across all samples from that patient, and there were cases when even two specimens of the same tumor showed different subtypes, indicating that spatial and temporal tumor heterogeneity does occur.
Genomic landscape of endocrine resistance
WES revealed many somatic mutations and especially many CNAs. Previous work has shown that low-frequency subclones present at 1% to 5% in the primary tumor can be enriched to greater than 40% in related metastases (41). Despite this some have used cutoffs of VAF to exclude low-frequency variants (42, 43). To minimize false-positives while maintaining high sensitivity, we followed a computational reinterrogation approach described by our group previously (34). Briefly, high-quality somatic variants called from the multiple samples of the same patient were combined into a single file per patient. The combined file containing all variants from the same patient were then reinterrogated across all samples from that patient. This method greatly improved detection of variants of low VAF and increased the percentage of shared variants across samples from each patient (Supplementary Fig. S2). Using the reinterrogated variants set, we explored the somatic mutation landscapes (Supplementary Table S4).
Genes previously shown to be recurrently altered in breast cancers were prevalent in this cohort. PIK3CA had the highest mutation rate with occurrence in 15 out of 35 patients. The other frequently mutated genes in our set included GATA3, CDH1, TP53, and ESR1, with mutation occurrences over 7 patients. The frequency in gene alterations reflected the multiple timepoints for a patient and/or multiple biopsies at a given timepoint (Fig. 2A). Of note, PIK3CA hotspot mutations were observed recurrently in all 3 patients that had de novo endocrine resistance, providing supporting evidence for previous findings that demonstrate hyperactivation of the PI3K pathway can promote endocrine therapy resistance (ETR; Supplementary Fig. S3; ref. 44). Mutations in the ESR1 gene have been extensively linked to ETR with the missense mutation D538G being the most prevalent (45, 46); this mutation confers ligand-independent constitutive activation. Consistent with previous findings, ESR1 D538G had the highest frequency in our cohort, showing a pattern of enrichment in resistant tumors (Fig. 2B). It was observed in none of the intrinsic resistant patients and in 5 acquired resistance patients during treatment with aromatase inhibition. Of these, 1 patient had the variant in the primary tumor with a very low VAF which increased in frequency in the subsequent endocrine-resistant tumor; 1 patient did not have the variant in primary tumor and acquired the variant in the subsequent resistant tumor; and the remaining 3 patients had the variant in resistant tumors with unknown status in primary tumors (DNA sequencing not available). Interestingly, we also found an ESR1 E380Q variant present at diagnosis in 1 responsive patient with ongoing response to letrozole. The VAFs gradually decreased in subsequent samples, in line with this mutants’ estradiol hypersensitivity nature (47). GATA3 was the second highest mutated gene, with a range of detected variants (Fig. 2C). GATA3 is an ESR1-cooperating transcription factor with frame-shift mutations being most commonly reported, however, there have been few studies of its potential role in endocrine resistance (48, 49). Recently Takaku and colleagues demonstrated that a specific type of GATA3 mutations in the second zinc-finger region reprogrammed progesterone receptor signaling (50). Accordingly, we saw mostly frame-shift mutations in our cohort, with four variants residing in the second zinc finger region and two acquired in resistant tumors (Fig. 2C).
To identify genes enriched in resistant tumors genome wide, we focused on the patients with DNA sequencing matched pre- and on-treatment samples who acquired resistance in the neoadjuvant setting (N = 12) and performed a matched mutation analysis (i.e., primary vs. matched later resistant samples per patient). We constructed count tables for each gene summarizing its mutation status in any of the sensitive and resistant tumors for each patient. Using a previously identified Pan-Cancer significantly mutated gene list (51), we found a variety of driver genes that were acquired in resistance in at least 1 patient (Supplementary Table S2). ESR1 appeared at the top of the list with 2 patients having the mutation acquired in resistant tumors. We identified 38 such genes including GATA3.
DNA CNAs including chromosome arm and subarm level changes were called using SynthEx (Supplementary Table S5; refs. 30). The copy number gain and loss frequencies were comparable with The Cancer Genome Atlas (TCGA) ER+ cohort (4). For example, we saw frequent 1q gain and 16q loss, noting that 16q loss is specific to luminal/ER+ breast cancers (Fig. 2D). Across the whole dataset, there were some differences between copy number landscapes of sensitive and resistant tumors, however, using a linear mixed model comparing copy number values between sensitive (N = 18 from 13 patients) and resistant tumors (N = 20 from 13 patients, see Methods for detail), and accounting for patients in the acquired resistance patient group, no significant regions of loss or gain were identified. This may be due to our small sample size.
We calculated the weights of known DNA mutational signatures (COSMIC v2) for each sample using DeconstructSigs. We then investigated both known DNA copy number changes that confer endocrine resistance, i.e., ERBB2 and FGFR1, and mutational signatures that are important in breast cancer, i.e., age-, apobec-, and HDR-related signatures in each patient that acquired endocrine resistance (Supplementary Fig. S4). Results demonstrated that few patients showed potential resistance mechanisms related to investigated CNA and mutational signatures. For example, patient 020 showed acquired ERBB2 amplification in resistant samples and patient 022 showed acquired FGFR1 amplification in one of the resistant samples. These selected genetic features also showed the vast genetic heterogeneity among patients and spatial biological replicates.
Differential transcriptional program in endocrine-resistant tumors compared with paired sensitive tumors
To identify genes differentially expressed in resistant versus sensitive tumors, we used linear mixed models to compare matched sensitive tumors with resistant tumors accounting for patients using the RNA-seq data. Briefly, biological replicates of the same tissue were first collapsed to construct a pseudo-sample by taking mean expression values across samples for each gene. We used RNA-seq data from 19 patients with acquired resistance and did matched comparisons of the untreated primary and the same tumor after acquisition of resistance. The t-statistic derived from the linear mixed model defines how consistently a gene is altered comparing resistant and sensitive tumors from each patient. The labels indicating treatment response were randomized 100 times to calculate FDRs.
We found 316 upregulated genes and 123 downregulated genes in resistant tumors that have FDR < 5% and fold change > 4 (for upregulated genes) or < 0.25 (for downregulated genes; Supplementary Table S6). Hierarchical clustering of these significantly expressed genes across the whole dataset encompassing all patients and samples largely separated sensitive and resistant tumors (Fig. 3A; Supplementary Table S7). Primary tumors and biopsy samples from tumors that retained endocrine sensitivity clustered together and resistant tumors also clustered together. Gene ontology (GO) analysis revealed that the upregulated genes in acquired resistant tumors were associated with cell proliferation, with GO terms chromatin assembly, nucleosome assembly, and DNA packaging (Fig. 3B). We then checked a known proliferation score and it was higher in resistant tumors among the 19 patients (Fig. 3C; ref. 52).
To fully characterize tumor transcriptional portraits, we used a panel of 654 previously published gene expression signatures measuring a variety of tumor phenotypes including amplicon signatures, oncogenic pathways, proliferation, and the tumor microenvironment (Supplementary Table S2; ref. 26). We calculated gene expression signature scores for each tumor and performed similar differential expression analysis using these signatures features. Among the acquired resistance tumors compared with their matched untreated primary, we found 37 upregulated gene signatures with FDR < 0.01 and fold change >1.15, and 3 downregulated gene signatures with FDR < 0.05 and fold change < 1 (Supplementary Table S8). Hierarchical clustering of differentially expressed signature scores again separated resistant and sensitive tumors (Fig. 4A; Supplementary Table S9). In particular, PAM50 subtype correlation to LumA centroids signature was one of the 3 downregulated signatures, concordant with the fact that resistant tumors were enriched with HER2-enriched and LumB subtypes. Among the upregulated signatures include two amplicon signatures, i.e., 16q23 and 8q, hypoxia, and several oncogenic pathway signatures; both estrogen and HER2 signaling pathways were also upregulated. A signature composed of GATA3-induced genes was also upregulated in resistant samples. In addition, AKT, RAS, and p63 pathways showed higher expression in resistant tumors, consistent with previous findings that growth factor signaling pathways play important roles in endocrine resistance (Fig. 4B; ref. 53). These results demonstrate resistant tumors harbor many transcriptional program changes compared with sensitive tumors.
Molecular portraits of selected patients
Finally, we sought to integrate our DNA and RNA analysis for each patient to identify the mechanisms of resistance in individual patients. Patient 006 was responsive to letrozole treatment, and the tumor volume changes clearly showed a good response to treatment (Supplementary Fig. S5A). VAF plots showed a decreased number of variants and lower VAF likely due to low tumor purity, although cancer cell cellularity remained high as assessed by a pathologist (Supplementary Fig. S5B). This patient had ESR1 E380Q hypersensitivity variant in the primary tumor and the following two samples had a lower VAF in later tumor specimens, suggesting selective death of these endocrine-sensitive cells (Supplementary Fig. S5C). E380Q is characterized as hypersensitive to estrogen ligand, which likely contributed to this tumor's responsiveness. Accordingly, proliferation rate was initially high in the primary tumor, which subtyped as LumB, and decreased in subsequent LumA tumor specimens from this patient (Supplementary Fig. S5D).
Patient 020 acquired resistance and had only received letrozole treatment (Fig. 5A). The primary tumor and the second specimen were taken while the patient was responding, while the third and the fourth specimens were collected when the cancer was growing and thus was resistant. A number of somatic mutations were detected (Fig. 5B), and there was a clear clonal shift from the second sensitive specimen to the fourth resistant specimen, from which there were five biological replicates (Fig. 5C). This patient gained both ESR1 and GATA3 mutations along with many other mutations including CHD4, CDH1, and PIK3R1. Gene expression showed a decreased proliferation rate in the second sensitive tumor that increased in the subsequent resistant tumors. Molecular subtype switched from primary tumor of LumB to LumA and back to LumB (Fig. 5D). Copy number landscapes showed increased genome instabilities in resistant specimens as well (Fig. 5E).
Finally, patient 028 had high ERBB2 mRNA expression and had a HER2-enriched gene expression subtype but lacked clinical HER2 status. This patient had multiple lines of endocrine therapy including fulvestrant, tamoxifen, and letrozole. Detected somatic variants showed clear clonal shifts (Fig. 6A and B), and acquired mutations included FGFR3 and GABRA6. There was also expansion of an existing clone in primary tumor containing mutations of PTEN, PIK3R1, and TP53, which are all recurrent variants. In terms of gene expression, all samples were classified as HER2-enriched subtype, with high ERBB2 and FGFR4 mRNA expression. Recently Garcia and colleagues showed that high FGFR4 expression plays an important role in HER2-enriched breast cancers and was associated with disease progression (54). Proliferation was lowest in the primary tumor and increased in later resistant tumors (Fig. 6C). The resistant tumor samples showed higher levels of copy number gains and losses (Fig. 6D).
ETR remains the biggest barrier to prolonged survival and cure for patients with ER+ breast cancer. There have been many studies investigating mechanisms of resistance and several have been defined, including acquired ESR1 mutations (14, 53). Here we sought to understand more about the development of endocrine resistance in a series of ER+ patients during hormonal therapy treatment utilizing integrated DNA sequencing and RNA-seq. We identified mutations that were enriched in resistant tumors, and we identified transcriptional profiles specific to tumor resistance, thus providing a molecular mechanism of resistance for some patients. This study is unique in its ability, through repeated and serial tissue acquisition during different phases of endocrine sensitivity within individual patients, to examine real-time alterations contributing to endocrine resistance.
Our work confirmed previous findings related to endocrine resistance mechanisms. Mutations enriched in resistant tumors include ESR1, GATA3, and FGFR3. Multiple signaling pathways were activated in resistant tumors including HER2, AKT, and RAS signaling, leading to increased proliferation rates within resistant specimens. Due to our sample size, larger cohort studies are needed to confirm mutations we identified that were only acquired by 1 patient, and increased power may allow detection of recurrent copy number changes.
Studies of sensitive breast cancers showed that clones that are sensitive to endocrine therapy can disappear under selective drug pressure. This was clearly seen in the patient with a sensitizing ESR1 E380Q variant that was present in the primary tumor and gradually disappeared in later samples. The genomic complexity seemed to reduce in sensitive cancers during treatment with a possible explanation being reduced cancer cellularity although the samples analyzed contained significant amounts of tumor on histological assessment. It is well known that multiple clones can exist within a primary cancer. Our observations show some clones are eliminated or reduced while others present at a low VAF at diagnosis increase and new clones with new mutations not detectable evolve possibly under the pressure of drug treatment.
A powerful aspect of this study is our unique series of patients that included patients who remained sensitive during treatment, and others with acquired resistance, with serial sampling over time from the same tumor in vivo. This study shows very clearly that there is no one mechanism that explains resistance in even a majority of patients. Each patient was unique in genetic and transcriptomic profiles. The most common ESR1 mutation only accounts for the resistance mechanism for 2 out of the 19 patients with acquired resistance by paired sensitive/resistant tumor analysis, and suggests that, although an important mechanism of resistance, it accounted for only a small number of resistant tumors in this study.
We collected multiple biological replicate samples from the same resistant tissue to assess tumor spatial heterogeneity. We observed clear heterogeneity between biological replicates in genetic and transcriptomic features. The biological replicates had both shared and unique mutations (Figs. 5C and 6C). As we used computational reinterrogation to recover low-frequency variants in related specimens, the unique mutations specific to each replicate are likely truly unique clones. Overall, gene expression profiles were highly correlated, and hierarchical clustering showed that replicates tended to cluster together (Fig. 3A). Subtle differences in gene expression were observed (Fig. 5D), with these findings suggesting that spatial heterogeneity (i.e., intratumor heterogeneity) could play a role in resistance to endocrine therapy. This study has confirmed some known mechanisms of resistance and identified some novel mechanisms. Most importantly it has highlighted the heterogeneity that underlies breast tumor biology and shows that endocrine resistance is complex and that more studies of patient-derived breast cancers are required. Taken together, it is clear that multiple molecular mechanisms underlying ETR exist and there is a need to understand the mechanism in individual patients to effectively combat resistance in each patient.
C. Martinez-Perez reports grants from Breast Cancer Now during the conduct of the study. J.S. Parker reports other support from Reveal Genomics outside the submitted work; in addition, J.S. Parker has a patent for PAM50 with royalties paid from Veracyte. C.M. Perou reports grants from Breast Cancer Research Foundation during the conduct of the study, as well as personal fees from Bioclassifier LLC outside the submitted work; in addition, C.M. Perou has a patent for U.S. Patent No. 12,995,459 issued, licensed, and with royalties paid from Bioclassifier. A. Turnbull reports grants from Breast Cancer Now and Lyda Henderson Fund at the Edinburgh Lothians Health Foundation during the conduct of the study. No disclosures were reported by the other authors.
Y. Xia: Formal analysis, investigation, visualization, methodology, writing–original draft, writing–review and editing. X. He: Methodology, writing–review and editing. L. Renshaw: Resources, data curation, writing–review and editing. C. Martinez-Perez: Data curation, methodology, writing–review and editing. C. Kay: Methodology, writing–review and editing. M. Gray: Methodology, writing–review and editing. J. Meehan: Methodology, writing–review and editing. J.S. Parker: Formal analysis, investigation, visualization, methodology, writing–review and editing. C.M. Perou: Conceptualization, resources, supervision, funding acquisition, writing–original draft, writing–review and editing. L.A. Carey: Conceptualization, resources, supervision, funding acquisition, writing–original draft, writing–review and editing. J.M. Dixon: Conceptualization, resources, data curation, supervision, funding acquisition, writing–original draft, writing–review and editing. A. Turnbull: Conceptualization, formal analysis, supervision, investigation, visualization, methodology, writing–original draft, writing–review and editing.
We would like to acknowledge the contribution and support of our friend and colleague Dr. Andrew Sims, University of Edinburgh, who sadly passed away prior to submission of this journal article. We are deeply grateful to the patients who allowed us to sample and use their tumor samples for this study. Thanks to Breast Cancer Now who funded the Translational Oncology Research Group at the Edinburgh Cancer Research Centre (J.M. Dixon, A. Turnbull). Funds from a generous grant from the Lyda Henderson Fund held at the Edinburgh and Lothians Health Foundation also supported this research. Funds from a generous grant from the Chewning Family to L.A. Carey also supported this research. This work was also supported by Breast Cancer Research Foundation (C. Martinez-Perez, L.A. Carey). Samples were collated and identified with help from the Edinburgh Breast Unit team and Cancer Research UK who fund the Edinburgh Tissue Group.
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.
Note: Supplementary data for this article are available at Clinical Cancer Research Online (http://clincancerres.aacrjournals.org/).