Cohort characteristics

A cohort of 736 fresh frozen metastatic lesions from unique BC patients was available for analyses (Additional file 1: Figure S1). The ER status was known for the primary disease in 700 cases and showed 568 ER + cases (81%). The biopsies of the metastases were mainly taken from the liver (45.2%), lymph nodes (19.6%) and bone (12.2%). Using the WGS data of this cohort, somatic single- and multiple-nucleotide variants (SNV/MNV), insertions/deletions (InDels), copy number variants (CNV) and structural variants (SV) were obtained. Based on CNV data of ERBB2, 68 MBC cases were considered as HER2-amplified (HER2+, see methods and Additional file 1: Figure S2).

The tumor mutational burden (TMB, the number of SNV/MNV/InDels per Megabase) of the HER2+ cases ranged from 0.65 to 32.51 with a median of 4.37 (95% confidence interval (CI) of 3.54–5.17). The number of SVs ranged from 78 to 2181, showing a median of 599 (95% CI 494.33–706.01). In line with the literature [25], the TMB and number of SVs were significantly higher in HER2+ versus HER2− tumors (MWU p = 0.0024 and p = 6.977e−12, respectively); however, this observation was partly driven by ER status as differences were less pronounced when only investigating ER-negative tumors (ER−/HER2+ vs ER-/HER2: MWU p = 0.64 and p = 0.011 for TMB and SV, respectively). However, when considering a cutoff of 10 mutations per MB, which defined as TMB-high in the literature [25], a small proportion of 6% of samples adhered to this definition. This was not significantly different from the proportion of TMB-high samples in the ER+- and triple-negative subtypes (TMB > 10/MB in 11% and 12%, respectively, Chi-square test p = 0.37).

Metastatic HER2-positive tumor tissue is comparable to HER2-positive primary breast cancer

To evaluate potential differences between unpaired primary and metastatic HER2+ tumor tissue, several genomic characteristics were compared: the TMB and number of SVs, putative driver mutations (here defined as amino acid changing events), mutational signatures (including Single and Double Base Substitutions, Insertion/Deletion (Indel) and SV signatures, respectively, SBS, DBS, ID and SVsig). The frequencies of 68 HER2+ MBC cases were compared to the 73 HER2+ PBC cases from the BASIS cohort [24] (showing a balanced proportion of 37% ER-negative cases for both cohorts, Fisher exact test p = 0.95).

In metastatic lesions, the median number of somatic nucleotide variants (including Indels) was 11,845.5 (interquartile range (IQR): 7364.5–19,398) and the median number of SVs was 573.5 (IQR: 339–912.5) which both were significantly higher than the number of events in HER2+ PBC (median, IQR and MWU for somatic nucleotide variants were 4702 (2935–7481) p = 1.318e−9 and for SVs were 172 (88–329), p = 2.014e−14; Fig. 1a).

Fig. 1

Differences between MBC and PBC. Levels and frequencies in unpaired PBC (green) and MBC (orange) for a total number of variants, being single-nucleotide variants and indels, and number of SVs, b % of relative contribution of the DBS2 mutational signature in HER2+ (left) and HER2− (right), showing that enrichment of the DBS2 signature in metastatic tumors is subtype-specific. p values derived from FDR-corrected MWU tests

Due to this significantly higher overall TMB in MBC, differences in mutation frequencies of individual genes between unpaired MBC and PBC were investigated taking TMB into account in a multivariate model. After correcting for multiple testing, besides a higher TMB (regression coefficient (95%CI) 0.027 (0.009–0.044), p = 0.003) a higher frequency of TP53 mutations was observed in MBC tissues (63% vs 36%, false discovery rate (FDR) corrected Fisher exact test p = 0.028), regardless of ER status (regression coefficient (95%CI) 0.266 (0.109–0.423), p = 0.001). No other significant differences in mutation frequency were observed. Next, a post hoc analysis of TP53 mutation frequency differences between MBC and PBC in other clinical subtypes—ER positive, HER2 negative and triple negative—also showed enriched TP53 mutations in MBC compared to PBC (for both subtypes MWU p < 0.0001), indicating enrichment of TP53 mutations was not specific for the HER2+ subtype. Of note, ERBB2 mutations were numerically higher in HER2+ MBC (9 out of 68) compared to HER2+ PBC (2 out of 73, uncorrected Fisher exact p = 0.027), but this failed to reach statistical significance after correcting for multiple testing (p = 0.35).

Next, differences in predefined COSMIC mutational signatures [15] between HER2+ MBC and HER2+ PBC were compared for ER+- and ER-negative cases separately. Out of all possible signatures, those with at least 10% contribution in at least 10% of the MBC samples were evaluated (n = 26, see Additional file 1: Table S1). Only DBS2 appeared to be enriched in HER2+ MBC compared to HER2+ PBC regardless of ER status (FDR corrected MWU p = 0.0057, Fig. 1b). Though significant higher, the median contribution DBS2 in MBC samples of 13.5% (vs 7.4% in PBC) appeared modest. In summary, besides a general increase in the overall number of somatic changes, no evidence was found for specific somatic changes associated with progression of HER2+ disease.

Finally, potential genomic scarring due to anti-HER2 treatment administered prior the biopsy was evaluated. HER2+ samples with prior anti-HER2 therapy (n = 32) were compared to HER2+ samples with no prior treatment (n = 20). Analyzing mutation frequencies and mutational signatures, with and without stratification for ER status, showed no significant differences after multiple testing correction.

Genomic features of therapy resistance

Next, WGS data were explored for factors that could predict response to HER2-targeted therapy. To this end, genomic characteristics of the metastasis were associated by Cox regression analysis to response in patients who received anti-HER2 therapy after tissue biopsy (n = 69, of which 44 ER+, 22 ER negative and 3 with unknown primary ER status. These patients were not necessarily HER2+ according to our genomic pipeline). For 33 patients (48%), this was the first line of therapy for metastatic disease. Treatment following biopsy included trastuzumab-based therapy (n = 66), T-DM1 (n = 2) and lapatinib (n = 1). Survival was defined from the start of treatment to radiological progression. Otherwise, patients were censored at the last radiological evaluation. Items tested in univariate Cox regression analysis (95 items in total) were the number of mutations (all somatic nucleotide changing variants), number and type of structural variants, WGS-estimated genome ploidy, mutational signatures (those with at least 10% contribution in at least 10% of samples), whole-genome duplication (WGD) and chromothripsis status, genes mutated (amino acid changing variants) and recurrent CNV regions (identified using GISTIC). Furthermore, ER and HER2 status, as well as the number of prior therapy lines were included (full list in Additional file 1: Table S1) In total, 13 items were significantly associated with progression-free survival time in univariate analyses (PFS, Cox regression p < 0.05), of which 5 remained significant in a multivariable model following forward selection: the number of lines of prior treatment, mutations of PIK3CA or CDK12, gain of chromosome 8p11.23 and a high contribution of gene signature DBS3 (Table 1).

Table 1 Univariable and multivariable Cox model for progression-free survival

A single risk score was calculated by combining the five items in the model (see “Methods” section), and the optimal cut point of this score was determined to classify patients as high or low risk. These groups were significantly associated with PFS (log rank p < 0.001) on anti-HER2 therapy (Fig. 2a).

Fig. 2

Multivariable model prediction of PFS. Kaplan–Meier progression-free survival curves of patients grouped by risk score. a Patients from our own cohort, labeling samples with a risk score ≥ 2.25 as high risk (see “Methods”) using five features: the number of lines of prior treatment, mutations of PIK3CA or CDK12, gain of chromosome 8p11.23 and a high contribution of gene signature DBS3. b Validation cohort, using PIK3CA, CDK12 mutation and Amp-peak 6 gain in the model as the only available features for the validation samples. Since these 3 features were mutually exclusive in this cohort, if any events of these 3 features were present the sample was assigned as high risk, the sample was assigned as high risk. p values derived from log-rank test

To validate these findings, publicly available data [11] documenting response to anti-HER2 therapy were used. From this validation cohort, we selected all metastatic breast cancer samples, which were sequenced by the MSK-IMPACT panel, covering 341 cancer-associated genes. Unfortunately, because of this targeted sequencing, DBS3 could not be reliably estimated, and furthermore, the number of prior therapy lines was unknown. However, the remaining three items (PIK3CA and CDK12 mutation status, and gain of 8p11.23) could still be combined (Fig. 2b) and again showed a clear association with PFS, also when we applied this score of the three parameters to our own cohort (log-rank test p < 0.001).

Establishing a HER2-driven expression profile

For a subset of 366 MBC samples (of which 49 HER2+), RNA sequencing was performed. To further investigate differences between HER2+ cases, a HER2-driven expression profile of 878 genes was established, which were identified by comparing HER2+ versus HER2− cases (see “Methods”, Additional file 1: Table S2). Of note, expression of ERBB2 itself and other chromosome 17 genes (which are potentially co-amplified with ERBB2) were excluded from this gene list. Gene expression values were used to create a correlation matrix of sample versus sample which was subsequently used for hierarchical clustering (Fig. 3, Additional file 1: Figure S3). The first (red) and last (blue) of the four cluster groups consisted of mainly ER+- and HER2-negative samples. The majority of the samples in cluster 2 (magenta) were ER negative, HER2-negative samples. Cluster 3 (yellow) was significantly enriched for HER2+ samples (47 out of 49 HER2+ samples, Fisher exact p = 1.94e−22) and prior anti-HER2 therapy treated cases (38 out of 49, Fisher exact p = 2.97e−12). The enrichment of HER2+ cases in cluster 3 was also evident when considering the PAM50 derived molecular subtypes [26] (62 out of 70 PAM50 HER2 subtypes are in cluster 3, Chi-square p = 8.24e−56). Remarkably, 16 of the 18 mutated ERBB2 samples in this dataset were also in cluster 3 (Additional file 1: Table S3, Fisher exact p = 1.17e−5) even though 10 of these 16 were HER2− (CN below threshold).

Fig. 3

Integrative analysis of hierarchical clustering using HER2-associated genes. A gene expression-based correlation matrix was used to cluster 366 MBC cases (see Additional file 1: Figure S3 for the accompanying heatmap). Relationships with ERBB2 log2-CN (dashed line indicates CN of 2, red indicate samples labeled as HER2+ in our analysis), ERBB2 mutation (amino acid changing events in: Ex, Extracellular domain, Tk, Tyrosine Kinase domain), PAM50 molecular subtype estimation, anti-HER2 therapy (AHT) given prior to biopsy and the ER status are shown

To verify that the 878 genes used for the clustering indeed represent a HER2-driven expression profile, independent datasets of publicly available expression data of PBCs (n = 867) were used. Samples were similarly clustered using the HER2-specific genes and the three main clusters were associated with the PAM50 molecular subtypes and ERBB2 expression. This again showed a clear association of the HER2-driven expression profile and HER2+ samples (i.e. ERBB2 expression was significantly higher in the HER2-enriched cluster containing 102 out of the 112 HER2 molecular subtypes (Additional file 1: Figure S4a–c). In addition, the HER2-specific genes were analyzed using IPA (Ingenuity Pathway Analysis), focusing on the ‘upstream regulators’ part, as defined by IPA. This showed that among the genes in the HER2-driven expression profile, several overlapped with target genes of known kinases, growth- and transcription factors, the top being ERBB2 (p = 2.65e−5, activation z-score 2.339) with 63 target genes in the HER2-specific gene list (see Additional file 1: Table S4 for activated regulators with p < 0.001). Other notable significantly activated upstream regulators are NRG1 and EGF, both ligands for the ERBB family of genes, and RAF1, involved in the MAP kinase pathway.

Underlying biological processes of the HER2-driven expression profile

Combining these results, the 116 samples from cluster 3 (Fig. 3) are clearly driven by common pathways which are specific for, but not exclusively linked to, HER2−positivity since ERBB2 is not always amplified or over-expressed (Additional file 1: Figure S5a). Since the original ER status was derived from pathology reports of the primary disease, ESR1 expression of the metastatic samples was evaluated as well. Surprisingly, 42 samples in this dataset (11.5%) show a low expression of ESR1 while the ER status was reported as positive in primary disease, and 36 out of these 42 samples are located in cluster 3 (Fig. 4a). Of these 36 samples, 25 do not show ERBB2 amplification (median CN 2.01, 95% CI 1.90–2.62). In contrast, in the PBC BASIS cohort [24] ESR1 expression levels corresponded with ER pathology status (Additional file 1: Figure S3b), suggesting that the observed ER pos-to-low cases in MBC may be an acquired trait. The ESR1 module score, a previously described gene signature associated with an active ER pathway [27], corroborates this change in ER. The median level of the ESR1 module was lower in cluster 3 versus cluster 1 and 4 (Additional file 1: Figure S5b). This was specifically observed in the ER pos-to-low cases than in the concordant ER+ cases and higher than in the concordant ER-negative cases (Kruskal–Wallis (KW) p = 2.39e−37, Fig. 4b).

Fig. 4

Overview of cluster 3 samples, showing a shift toward increased MAPK pathway expression. a ESR1 expression in the 4 clusters, red indicates samples labeled as ER+ in their primary disease, blue ER negative in primary disease. b The ESR1 module score over samples, grouped by the meta ER status: a combination of MBC ESR1 expression and reported ER status in primary (i.e. ‘pos.to.low’ indicate samples with a reported ER+ status in primary but with low ESR1 expression in the metastatic lesion and vice versa for ‘neg.to.pos’). c, d MAPK signature expression over cluster groups (c) and meta ER status (d). e Summary of characteristics of cluster 3 samples. Top 2 tracks show log2 CopyNumber (dotted line indicates CN of 2) and mutation of ERBB2. The third track indicates patients who received anti-HER2 therapy prior to biopsy. Track 4 shows ER status: brown indicate cases of which the primary tumor was ER positive but the metastatic biopsy has low ESR1 expression, i.e. ER pos-to-low. Green: ER neg-to-high, orange: concurrent ER pos, gray: concurrent ER neg. Bottom track shows expression levels of the MAPK signature, and dotted line indicates median

Although the absolute percentage of ER+ cells of the primary tumor is unknown and may explain this observation, an attractive alternative hypothesis for these originally ER+ samples with low ESR1 expression, is that the tumor cells may have switched from ER-dependent growth to using alternative pathways. Since the majority (86%) of these samples cluster together with samples previously exposed to anti-HER2 treatment, the MAPK pathway is of particular interest since it has been linked to both anti-HER2 therapy and endocrine resistance [11, 28]. Somatic mutations labeled as pathogenic (according to OncoKB) in either ERBB2, EGFR, NF1, KRAS, BRAF and MAPK2 were first used as biomarker for activated MAPK signaling [11]. Of the 23 cases with such a mutation, 12 were located in cluster 3: 9 ERBB2 and 3 EGFR. All 6 KRAS/BRAF-mutated cases were located in cluster 1 (along with 3 EGFR and 2 ERBB2 cases) and that cluster as a whole showed a lower MAPK expression signature (Fig. 4c). This is in line with results by Wagle et al. [27], showing a disconnection between RAS/RAF mutation status and MAPK pathway activity as determined by a MAPK expression signature.

Using this established MAPK signature, we found significantly higher expression in cluster 3 as a whole (KW p = 6.25e−10, Fig. 4c), but this was specifically observed in ER pos-to-low cases (KW p = 7.65e−5, Fig. 4d). A post hoc analysis showed a significantly higher MAPK signature (MWU p = 0.00088) in the ER pos-to-low samples compared to concordantly ER+ cases (i.e. primary ER+ status with high ESR1 expression in the metastatic lesion). In PBC, the expression of the MAPK signature was significantly higher in both the HER2 cluster and the luminal-driven cluster, compared to the basal cluster (Additional file 1: Figure S4d, MWU p = 4.76e−9 and p = 3.20e−21). Considering results in MBC, which showed the MAPK signature most active in cluster 3 (Fig. 3d), this suggests that in MBC MAPK activity is mostly sustained in HER2-driven samples.

As summarized in Fig. 4e, in our metastatic BC cohort 116 samples show a HER2-driven expression profile. Forty-seven of these had a high ERBB2 CN, while of the remaining 69 cases, 10 had a mutation in ERBB2 and 9 had received prior anti-HER2 therapy (3 cases were both mutated and had received prior anti-HER2 therapy), but were progressive on this therapy. Of the 53 ERBB2 low CN without an ERBB2 mutation and not having received prior anti-HER2 therapy, 21 potentially switched from being ER-driven (ER+ according to status in primary) to HER2-driven disease (metastatic biopsy had low ESR1 expression but co-clustered with HER2 cases) and a further 32 samples may be driven by MAPK (above median MAPK signature expression), with 12 samples showing both phenotypes (low ESR1/high MAPK).