Graphical abstract graphic file with name fx1.jpg [37]Open in a new tab Highlights * • Deep (phospho)proteomics profiling of pre-treatment breast cancer needle biopsies * • Accurate quantification of HER2 and ER abundance improves outcome predictions * • Low tumor infiltration of immune cells is associated with poor treatment response * • Expression levels of HER2, HER4, ER, IGF1R, and Kalirin predict treatment response __________________________________________________________________ Resistance against targeted therapy in HER2^+ breast cancer patients is prevalent. Debets et al. perform (phospho)proteomics profiling of early-stage HER2^+ breast cancer biopsies to identify signatures predictive of treatment response. Resistance is multifaceted, and combining receptor protein quantification with corresponding signatures of activation enables better prediction of patient response. Introduction Invasive breast cancer (IBC) is a highly heterogeneous disease that relies on subtype classification to prognosticate the disease course and to select treatment strategies.[38]^1 Current IBC patient classification focuses on the expression of three receptor proteins: the estrogen receptor (ER), the progesterone receptor (PR), and the human epidermal growth factor receptor 2 (HER2). Subtype-specific treatment strategies are directed at these receptors and their downstream signaling pathways because these are considered to drive tumor progression. The efficacy of these subtype-driven therapeutic interventions is limited by our ability to further classify IBC within the main subgroups and identify tumors that are truly biologically driven by these targeted receptors. Approximately 15% of all IBC tumors overexpress HER2.[39]^2 The current standard neo-adjuvant treatment for HER2^+ IBC applies targeted therapy using the monoclonal HER2-directed antibodies trastuzumab (TTZ) and pertuzumab (Ptz) in combination with conventional chemotherapy drugs, such as taxanes and carboplatin.[40]^3^,[41]^4^,[42]^5^,[43]^6 The introduction of these targeted therapies has improved the clinical outcome for this group of patients considerably, yet treatment resistance, both intrinsic and acquired, occurs because 20%–40% have only partial response to neo-adjuvant treatment.[44]^7 A full understanding of the molecular mechanisms underlying treatment resistance in vivo is lacking. Although there has been progress in identifying genomic and transcriptomic features that predict treatment success of HER2-targeted therapies, no biomarkers are currently used in clinic.[45]^8 The HER2^+ subgroup represents a remarkably heterogeneous population of tumors.[46]^9 Poor treatment response in some patients has been attributed to this heterogeneity, because it represents a mismatch between tumor biology and applied therapeutics. For example, HER2^+ tumors co-express varying levels of ER, wherein ER expression correlates with treatment outcomes for HER2-targeted therapy.[47]^10 Therefore, enhanced patient stratification could improve clinical outcomes by identifying more suitable therapeutic strategies in some cases and reduction of overtreatment in others. Deep tumor profiling is required to enhance IBC classification and uncover the in vivo biological complexity and diversity of treatment resistance in HER2^+ tumors. Delineation of the dynamic and complex cellular networks in individual tumor samples opens the door to precision oncology. In the current study, we perform proteomics and phosphoproteomics profiling of 45 pre-treatment biopsies of patients with early-stage HER2^+ IBC (37 ER^+ and 8 ER^− cases) to identify molecular signatures predictive of treatment response to neo-adjuvant carboplatin, paclitaxel, TTZ, and Ptz. We demonstrate the feasibility of microscale clinical proteomics; we present deep (phospho)proteomics profiling with high data quality using very limited sample input (half of a 14G needle biopsy). Our data show that IBC subtype classification by accurate quantification of IBC biomarkers combined with the assessment of associated biological features improves treatment outcome prediction. We demonstrate an enhanced IBC classification scheme where signatures of biological activity are used in combination with receptor protein expression data to identify tumors that are HER2 or ER driven. Furthermore, we identify multiple cellular mechanisms that precondition tumors to resist therapy: unfolded protein response (UPR) induced cellular dormancy, a metabolic switch toward oxidative phosphorylation (OXPHOS), and reduced numbers of tumor-infiltrating leukocytes (TILs). Together, the identified resistance mechanisms constitute to a strong signature associated with treatment success in this dataset. Results To gain insight into molecular signatures predictive of treatment response for paclitaxel, carboplatin, TTZ, and Ptz (PTC-Ptz), we analyzed 45 treatment-naive early IBC needle biopsies using our microscale (phospho)proteomics workflow ([48]Figure S1). All patient biopsies included in this study were classified as HER2^+ with a minimum tumor grade of 2 and tumor cell density of 60% or higher. After needle biopsies were taken, patients completed multiple drug treatment cycles according to the regimen of the TRAIN2 study.[49]^6^,[50]^11 Finally, treatment response was determined at surgery. The systemic treatment response was categorized as pathological complete response (pCR = ypT0/isypN0) (n = 26), near pCR (npCR) (n = 8), or No pCR (n = 11). If pCR was not achieved, samples were referred to as npCR if <10% of tumor remained and as treatment resistant (No pCR) if >10% of tumor remained. A complete list of patient biopsy details is provided in [51]Table S1. We obtained deep proteome and phosphoproteome coverage of patient samples using tandem mass tag (TMT) isobaric labeling of digested proteins from limited tissue material (less than half of a 14G needle biopsy, ca. 6–12 mg of fresh frozen tissue). Labeled samples were multiplexed into separate TMT 10-plex sets all sharing an identical pooled reference channel, enabling accurate quantification across TMT sets. Protein and phosphopeptide abundance values were normalized to this shared pooled reference channel within each TMT set. Altogether, 11,088 protein groups and 37,696 phosphopeptides were quantified across the 45 breast cancer biopsies, of which 9,340 protein groups and 11,234 phosphopeptides were identified in at least 75% of the samples ([52]Figures S1B and S1C; [53]Tables S2 and [54]S3). We observed strong longitudinal reproducibility of our peptide and protein quantification demonstrated by the high correlation of the pooled TMT reference channel measurements across TMT sample sets ([55]Figure S1D). The quantitative reproducibility was further determined by replicate samples that were measured across different TMT sets. These replicates showed high correlation in both their proteome and phosphoproteome measurements ([56]Figure S1E) and tight grouping during unsupervised clustering of our datasets ([57]Figure S1F). Finally, we performed unsupervised clustering to confirm that our data did not suffer from any TMT batch effects or clustering bias based on tumor percentages ([58]Figure S1F). Together, these data demonstrate the technical strength and quantitative reproducibility of our microscale (phospho)proteomics workflow. Enhanced IBC subtype classification improves treatment outcome prediction Low HER2 activity contributes to treatment resistance The efficacy of HER2-targeted therapy is partially attributed to the inhibition of HER2 signaling in tumors that are dependent on this pathway for cell growth and proliferation (i.e., HER2-driven tumors).[59]^12 Low activity of HER2 signaling has therefore been proposed as a resistance mechanism against HER2-targeted therapy.[60]^13 Accordingly, we anticipated that identifying tumors that are not HER2 driven would enable detection of poor treatment response in patients. To identify HER2-driven tumors, we first examined HER2 quantification within the patient biopsies. Current HER2 classification is based on immunohistochemistry (IHC) scoring supplemented by in situ hybridization if IHC provides an ambiguous score (2+). Tumors with a 2+ or 3+ score are eligible for HER2-targeted therapy. We compared our proteomics measurements with this classification and found that the IHC scores correlated with HER2 protein levels in general but were inconsistent for some tumors ([61]Figure S2A). We found that a 2+ IHC score could not sufficiently predict treatment response ([62]Figure S2B), whereas protein quantification of HER2 by proteomics measurement showed a strong correlation with pCR status ([63]Figure 1A). HER2 was significantly lower among the No pCR samples in our proteomic data (p < 0.05) ([64]Figure 1A). To validate our mass spectrometry (MS)-based measurements of HER2 abundance, we performed western blot analysis on a subset of patient samples and found a strong correlation (R = 0.89, p = 3.7e−6) in relative protein expression levels ([65]Figure S2C). Together, these data indicate that a more quantitative readout for HER2 protein levels, compared with current IHC semi-quantitative classification, may better discern between treatment outcome groups. Figure 1. [66]Figure 1 [67]Open in a new tab Pseudo-HER2^+ signature is associated with treatment resistance (A) Boxplot of HER2 protein expression by treatment outcome group. Protein abundances are normalized to the pooled reference channel. ∗Indicates p value < 0.05 in unpaired students t test. (B) Heatmap of unsupervised clustering of HER2 and adjacent genes PGAP3, STARD3, and GRB7 (data are Z scored). (C) Volcano plot comparing phosphosites between true HER2^+ subgroup vs. pseudo-HER2^+ subgroup. All phosphosites belonging to the HER2 protein are indicated in red. Difference values are log2 ratios between subgroups. (D) GSEA of HER2-related Reactome pathways and Gene Ontology. NES (normalized enrichment score) is shown on the x axis; the dot size represents the –log10(FDR). Pseudo-HER2 samples are shown in black. (E) Left: pCR rate (%) among the pseudo-HER2^+ subgroup (n = 12) compared with the full dataset (n = 45). Right: number of No pCR samples (%) in pseudo-HER2^+ subgroup compared with IHC 2+ subgroup. (F) Heatmap of unsupervised clustering of HER2 family members (EGFR, HER3, HER4), ER, and IGF1R expression levels groups the samples into two main clusters. Left cluster shows enrichment of pseudo-HER2^+ patients with No pCR and high levels of HER3, HER4, IGF1R, and ER. (data are Z scored). (G) Boxplot of HER4 protein expression by treatment outcome group (FDR = 0.02). ∗p < 0.05; ∗∗p < 0.01. To support the HER2 protein levels measured in our dataset, we evaluated additional features that indicate HER2 activity. Herein, we identified a subset of samples, referred to as “pseudo-HER2^+,” which do not appear HER2 driven despite being eligible for HER2-targeted therapy. First, we evaluated the expression of genes adjacent to HER2 in the genome (PGAP3, STARD3, and GRB7) as a proxy for regional genome activity. Pseudo-HER2 samples showed low protein expression of these cis genes, indicating reduced genome activity and consequently low HER2 expression, even in patients possessing a high number of HER2 genome insertions ([68]Figures 1B and [69]S2D). Expression levels of TOP2A, however, did not correlate to HER2 expression levels and were similar for the pseudo-HER2^+ subgroup ([70]Figures S2E and S2F). Second, we observed that the pseudo-HER2^+ subgroup displayed a lower activation state of HER2 as demonstrated by significant downregulation of HER2 phosphorylation sites (false discovery rate [FDR] < 0.05) ([71]Figures 1C and [72]S2G). Finally, we found that HER2 downstream signaling was decreased within the pseudo-HER2^+ tumors ([73]Figure 1D). We next assessed the pCR status of the patients in our pseudo-HER2^+ group and found that pseudo-HER2^+ classification was a far stronger predictor for treatment resistance compared with the HER2 IHC score ([74]Figure 1E). These data indicate that patient stratification based on a multi-level assessment of HER2 features improves predictions for patient response to treatment. We next set out to explore the differences between the pCR and No pCR patients within the pseudo-HER2^+ subgroup. We found that pseudo-HER2^+ tumors with high expression levels of HER2 family members (HER3 and HER4) and alternative hormone receptors (ER and insulin growth factor 1 receptor [IGF1R]) were generally treatment resistant ([75]Figure 1F). This finding suggests that these receptors could provide a compensating mechanism, bypassing HER2 dependency in No pCR patients within the pseudo-HER2^+ subgroup. For HER3, this is in contrast with previous research, in which low HER3 levels were linked to treatment resistance.[76]^14 Notably, HER4 was also found to be significantly upregulated (FDR = 0.02) among the treatment-resistant tumors within the full dataset (regardless of the pseudo-HER2^+ status) ([77]Figure 1G). This was accompanied by increased downstream HER4 signaling ([78]Figure S2H) and suggests a more widespread role of HER4 in treatment resistance. To validate our findings in a larger BC clinical dataset, we analyzed biopsy data from the I-SPY2 neoadjuvant trial (NIH identifier: [79]NCT01042379), which contains pre-treatment biopsies matched to patient outcomes for 10 different treatment arms.[80]^15 We extracted reverse-phase protein array (RPPA) and mRNA data for 43 HER2^+ patients who received paclitaxel + pertuzumab + trastuzumab (PPT) treatment, of which 25 patients achieved pCR. We first examined HER2 RPPA protein expression levels and found that, in agreement with our observations, HER2 abundance and HER2 phosphorylation sites were significantly lower in No pCR patients compared with pCR patients ([81]Figure S2I). We next evaluated the “pseudo-HER2” signature in patient mRNA data and found that biopsies possessing low expression of HER2 adjacent genes (PGAP3, STARD3 GRB7) were enriched for No pCR treatment outcomes ([82]Figure S2J), supporting our proteomics findings. Finally, we examined the expression of HER3 and HER4 family members and found significant mRNA upregulation of HER4 in No pCR patients ([83]Figure S2K). HER4 upregulation was even more pronounced when we evaluated treatment response within the “pseudo-HER2” patient subset ([84]Figure S2K). Interestingly, we did not find significant HER4 upregulation in RPPA protein data. We also did not observe any HER3 upregulation in No pCR patients in mRNA or RPPA measurements ([85]Figure S2M). High ER signaling is associated with treatment resistance, especially in combination with high IGF1R activity In addition to the determination of the HER2 status, IBC tumors are also classified based on IHC staining of ER. A negative ER IHC status has been associated with better treatment response.[86]^10 Accordingly, we found a higher pCR rate among ER tumors in our dataset ([87]Figure S3A). Positive ER classification by IHC scoring, however, did not show an enrichment of treatment-resistant patients ([88]Figure 2A). ER^− tumors are likely responsive to treatment, whereas the ER^+ population has a heterogeneous treatment outcome. We found that ER protein expression levels measured by MS were generally associated with the IHC status (FDR < 0.05). However, there was a large overlap in ER protein expression levels between the ER^− and ER^+ tumors ([89]Figure 2B). To validate our MS quantification, we performed anti-ER WB analysis on a subset of patient biopsies and found a strong correlation in our protein quantitation (R = 0.86, p = 1.5e−5) ([90]Figures 2C and [91]S3B). Together this demonstrates that even though a negative ER IHC status is generally predictive of pCR, ER^+ tumors have a mixed treatment outcome and thus treatment-resistant tumors cannot be detected by this classification. Furthermore, for some tumors the ER protein expression levels did not reflect the IHC classification. These observations indicate that further discrimination between ER^+ patients is needed to better identify non-responsive patients. Figure 2. [92]Figure 2 [93]Open in a new tab High ER signaling predicts poor treatment response, especially in combination with high IGF1R activity (A) pCR rate (%) for the ER-dependent subgroup (n = 13), tumors with an ER + IHC status (n = 37), and the full dataset (n = 45). (B) ER protein expression (determined by proteomics) by ER IHC status. ∗ FDR < 0.05 (C) Correlation plot of ESR protein abundance as established by western blotting (WB) compared with MS. (D) Heatmap showing increased ER signaling and increased levels of FOXA1 and GATA3 among a subgroup of tumors that were classified as ER dependent. (data are Z scored). (E) Correlation plot of the ER and IGF1R expression levels. Point color represents the treatment outcome. Pearson correlation: 0.78. (F) Heatmap of unsupervised clustering of IGF1R and phospho-IGF1R-T1366 groups four patients together with the highest expression levels. These tumors were resistant to treatment and ER dependent. (data are Z scored). (G) PTM-SEA analysis NESs averaged for the patients with a high IGF1R expression (in black) compared with the rest of the data. An increased enrichment score was found for SGK1 kinase activity and IGF1 perturbation among the tumors expressing high levels of IGF1R. We hypothesized that among the ER^+ tumors a subset of samples was truly ER dependent, resulting in a poor treatment response. To establish whether ER-dependent signaling could be discerned within the ER^+ tumors, we first analyzed the protein expression levels of ER and ER-associated proteins GATA3 and FOXA1. The results revealed that ER and GATA3 protein expression were significantly upregulated among the treatment-resistant tumors within the ER^+ subgroup ([94]Figure S3C) (p = 0.007 and 0.046, respectively). FOXA1 contrarily did not significantly differ between the outcome groups ([95]Figure S3C). Second, we quantified ER-dependent signatures in the tumor proteome using gene set enrichment analysis (GSEA). Combining the results together identified a subgroup of tumors displaying an ER-dependent signature, characterized by high expression levels of ER, FOXA1, and GATA3 and enrichment of ER-related signaling ([96]Figure 2D). This classification of tumors, based on accurate ER protein quantification by proteomics and ER-driven biological features, was a better predictor of treatment resistance than the ER IHC status ([97]Figure 2A). ER-dependent tumors showed >50% increase in No pCR patient enrichment compared with IHC ER^+ tumors. We next set out to investigate the differences between the pCR and No pCR tumors within the ER-dependent subgroup. We found that the correlation between the ER and IGF1R expression levels was very high and, in agreement with previous findings,[98]^16 that high expression levels of IGF1R and ER were correlated with treatment resistance ([99]Figure 2E). Unsupervised clustering identified a subset of four tumors that were ER dependent, treatment resistant, and exhibited the highest IGF1R and phospho-IGF1R levels ([100]Figure 2F). These elevated (phospho-)IGF1R expression levels were linked to increased IGF1R downstream signaling, as demonstrated by increased phosphorylation activity caused by IGF1 (an important IGF1R ligand[101]^17) and high SGK1 activity (a kinase downstream of IGF1R) ([102]Figure 2G). Combined, these data suggest that increased IGF1R activity within the ER-dependent subgroup may contribute to treatment resistance. This is in line with previous research that has suggested interplay between ER and IGF1R[103]^18^,[104]^19 and a correlation with treatment resistance.[105]^16 Our data revealed that patient stratification based on the quantitative assessment of ER by proteomics, combined with the analysis of ER-associated proteins (GATA3 and FOXA1) and ER-related signaling, identified a subset of tumors with an ER-driven signature. This signature was a better predictor for treatment resistance than ER IHC scoring in this dataset. Tumors with an ER-dependent signature and high IGF1R activity had an especially poor treatment response. To validate these findings in an independent patient cohort, we analyzed ER protein expression and ER, IGFR1, FOXA1, and GATA2 mRNA expression levels in the I-SPY2 neoadjuvant dataset.[106]^15 In alignment with our MS results, we found that patient tumors possessing high expression of ER, FOXA1, and GATA2 mRNA showed enrichment in No pCR response ([107]Figure S3D). We also saw a strong correlation between ER and IGFR1 mRNA expression levels ([108]Figure S3E). Lastly, we found that within hormone-positive tumors, ER protein levels were higher in No pCR patients than pCR patients ([109]Figure S3F). Treatment-resistant tumors are preconditioned to evade therapy Next, we performed an exhaustive analysis to identify further associations between molecular signatures and patient pCR status. We focused on those proteins showing statical significant difference within the patient samples for further bioinformatics analysis and literature mining. This resulted in the observation of several important biological processes observed in subsets of the patient samples that could be involved in therapy resistance. Increased OXPHOS activity preconditions tumors to resist therapy To determine whether treatment-resistant tumors are metabolically preconditioned to evade drug therapy, either targeted treatment or chemotherapy because these individual responses cannot be distinguished in our dataset, we looked for upregulation of metabolic pathways known to be associated with treatment resistance. Our data revealed a distinct upregulation of proteins involved in OXPHOS among two treatment-resistant tumors ([110]Figure 3A). GSEA showed a stark enrichment of OXPHOS activity among these same tumor samples ([111]Figure 3B), but not in any of the other tumors ([112]Figure S4A). OXPHOS metabolic reprogramming has been recognized as an emerging hallmark of cancer and is associated with acquired drug resistance.[113]^20^,[114]^21^,[115]^22 Although our observation is in a very small number of patients, the upregulation of the OXPHOS metabolic profile is exceptionally strong and unambiguous. Figure 3. [116]Figure 3 [117]Open in a new tab Increased OXPHOS activity preconditions tumors to resist therapy (A) Heatmap of unsupervised clustering of proteins involved in OXPHOS Reactome pathway grouped two patients together. They were treatment resistant and showed clear increase in OXPHOS proteins. (data are Z scored). (B) Volcano plot of GSEA of hallmark signatures of the two tumors marked with increased levels of OXPHOS. NESs plotted against the –log10(FDR). Hallmark term OXPHOS indicated in the box. (C) PTM-SEA analysis of AMPKA kinase activity. NESs are plotted on the x axis, dot size represents the FDR score, and color represents the metabolic switch. (D) Dotplot showing the log2 normalized phosphosite and protein abundance. Dot color represents the metabolic switch. The OXPHOS signature was accompanied by the highest expression levels of ER ([118]Figure 3A) and increased AMP-activated protein kinase (AMPK) activity ([119]Figure 3C). This agrees with previous research showing that OXPHOS can be activated by ER-mediated activation of AMPK in response to glucose deprivation.[120]^23 We also find increased phosphorylation of known AMPK targets Transcription Factor A (TFAM) and Acetyl-CoA carboxylase 1 (ACACA) ([121]Figure 3D). AMPK-mediated activation of TFAM and ACACA has been shown to stimulate OXPHOS and mitochondrial biogenesis in vitro and in vivo. 22[122]^24,[123]^25 The increased expression levels of mitochondrial proteins found in the two samples with high OXPHOS activity could be the result of increased mitochondrial biogenesis ([124]Figure S4B). Together these data suggest that ER-mediated activation of AMPK may be responsible for the increased OXPHOS found in these tumors. UPR-induced cellular dormancy preconditions tumors to resist therapy A significant decrease of the entire ribosomal machinery was found among five tumor samples, of which four were treatment resistant ([125]Figure 4A). These tumors showed decreased expression levels of proteins involved in translation, indicative of a dormant cell type ([126]Figures 4B and [127]S5). We hypothesized that this non-proliferative cell type may prove less responsive to cell-cycle-dependent chemotherapy, resulting in treatment resistance. Figure 4. [128]Figure 4 [129]Open in a new tab UPR-induced cellular dormancy preconditions tumors to resist therapy (A) Heatmap of unsupervised clustering of all ribosomal proteins grouped five patients together with the lowest abundance of these proteins. Four of these were resistant to therapy. (data are Z scored). (B) A t test was performed on the proteins within the dataset comparing the five tumor samples with the low-ribosomal proteins against the rest. Difference is plotted on the x axis, and −log10(FDR) on the y axis. Proteins in yellow are involved in translation initiation (Reactome pathway). (C) Boxplot of STC2 comparing the tumors expressing low amounts of ribosomal proteins vs. the rest. ∗∗p < 0.05. Dot color represents the treatment outcome. (D) A t test was performed on PTM-SEA kinases NESs, including all kinases that were found enriched (FDR < 0.05) in at least one of the low-ribosomal tumors. PDK1 was clearly upregulated among the low-ribosomal tumor samples. Among the tumors with this dormant cell type, we found increased levels of Stanniocalcin-2 (STC2) ([130]Figures 4B and 4C) (p < 0.05). Previous research has shown that a rapid upregulation of STC2 is associated with the UPR.[131]^26^,[132]^27 UPR is triggered by the accumulation of misfolded proteins within the endoplasmic reticulum and aims to reinstate cellular homeostasis by clearance of the endoplasmic reticulum and constriction of protein synthesis. UPR-induced cellular dormancy has been associated with chemotherapy resistance.[133]^28 The upregulation of STC2 among the dormant tumors suggests that UPR may be activated and responsible for the non-proliferative cell type. We also found increased PDK1 activity among the dormant tumors ([134]Figure 4D). Upregulation of PDK1 has been shown to circumvent endoplasmic reticulum stress-induced apoptosis in tumor cells.[135]^29 Hence the increased PDK1 activity found in these tumors could provide a pro-survival mechanism upon endoplasmic reticulum stress. These data indicate that increased UPR could result in a dormant cell type in a subset of tumors that were largely treatment resistant. UPR and endoplasmic reticulum stress have been linked to acquired therapy resistance against many chemotherapeutics.[136]^30^,[137]^31^,[138]^32 Our data suggest that these processes are not merely secondary resistance mechanisms but can also play a role in de novo treatment resistance. Low immune cell infiltration in the TME indicates poor treatment response Evasion of immune system clearance is a prerequisite for successful tumor progression, and excluding lymphocytes from the tumor microenvironment is one mechanism that enables immune escape. Accordingly, we anticipated that patients who respond poorly to therapy possess low levels of immune cell infiltration in their tumor biopsies. Using a validated selection of immune cell markers,[139]^33 we found that >80% of No pCR tumors had depleted levels of immune cells ([140]Figures 5A and 5B). Furthermore, we found that this “low immune infiltration” patient cohort showed strong downregulation of immune-related signatures within their proteome, whereas the opposite was true for “high immune infiltration” patients ([141]Figure 5C). This demonstrates that both a curated list of specific cell markers and a global analysis of proteome signatures identify the vast majority of non-responders as being immune depleted. We next validated this observation in the I-SPY2 neoadjuvant dataset by clustering patient mRNA data based on the Danaher et al.[142]^33 immune cell marker panel. Here we found that the cluster of patients showing low levels of TILs was composed of 70% No pCR patients ([143]Figure S6A). Figure 5. [144]Figure 5 [145]Open in a new tab Low immune cell infiltration and high KALRN levels correlate with poor treatment response (A) Unsupervised clustering of validated leukocyte markers shows a distinct grouping of tumors that contain reduced immune cells (Low TIL); this group includes nine No pCR patients. Data represent scaled abundances. (B) Downregulation of T cell marker proteins in the low TIL patient subset. Boxplot of T cell protein expression in all samples, grouped by clusters defined in (A). ∗∗p < 5e−9. Dot color represents the treatment outcome. (C) Downregulation of immune-related pathways in low-TIL sample subset. Dots represent individual GSEA pathways that were significantly enriched (FDR < 0.05) in at least 75% of the samples in each group. Black dots represent immune-related GSEA terms; x axis represents average NES for each group. (D) Upregulation of Kalirin in No pCR tumors within the low-TIL patient group. Volcano plot comparing protein abundance between low-TIL vs. high-TIL subgroup. Vertical bars indicate 2-fold change in protein abundance. Blue and red dots denote significantly downregulated or upregulated proteins, respectively. Larger dot size represents FDR < 0.1. (E) Differentially regulated phosphosites in No pCR tumors within the low-TIL group show enrichment of Rho GTPase pathway (FDR < 0.001). (F) ATM kinase activity is enriched in No pCR tumors within the low-TIL group. A t test was performed on PTM-SEA kinases NESs, including all kinases that were found enriched (FDR < 0.05) in at least three patients. The immune-depleted tumor subset we identified in our dataset was not exclusive to No pCR patients. Therefore, we next examined which features within the low-immune group were enriched specifically in treatment-resistant patients. We identified the guanine nucleotide exchange factor Kalirin as the one of the most upregulated proteins in the No pCR patient subset (>3-fold increase, p < 0.0003, FDR < 0.1) ([146]Figure 5D). Kalirin is known to regulate specific Rho GTPases and plays a vital role in neuronal plasticity.[147]^34 Interestingly, loss-of-function mutations in Kalirin were recently identified as a biomarker for positive immunotherapy response across 10 different cancer types.[148]^35 Functional Kalirin protected tumors from DNA damage through activated DNA-damage-repair mechanisms mediated by Rho GTPases, thus contributing to reduced neoantigens and poor immune clearance upon reactivation of the immune system.[149]^35 Accordingly, we found Rho GTPase signaling as the most differentially regulated pathway between No pCR and pCR samples within the low-immune subset in our phosphoproteomics data ([150]Figure 5E). We also identified ATM kinase, a central activator of DNA damage response, as significantly activated in No pCR patients compared with pCR patients in this subset ([151]Figure 5F). In contrast with previous findings, however, we identified a positive correlation between Kalirin and PD-L1 protein expression levels (p = 0.0014) in tumor biopsies[152]^35 ([153]Figure S6B). Together, our data demonstrate that non-responsive patients in the low-immune patient cohort have uniquely elevated levels of Kalirin and increased activation of associated downstream pathways. This finding suggests that higher Kalirin levels may participate in protecting tumors from immune cell clearance in No pCR patients, even after targeted therapy. Treatment resistance is multifaceted, and the combined resistance mechanisms are strongly associated with treatment success We aimed to assess whether the resistance mechanisms we identified constitute a protein signature associated with treatment outcome. We found that multiple mechanisms coincide in any given treatment-resistant tumor ([154]Figure 6A). Thus, rather than the alteration of a single protein or pathway, the combination of drug-evading mechanisms drives tumor fate. This highlights the multifaceted nature of drug resistance in vivo and illustrates the need to take a combination of these features into account to enable a more accurate prediction of treatment outcome. For example, a subset of treatment-resistant patients had both a pseudo-HER2^+ signature and showed ESR dependency. This group of patients might benefit from ESR-directed therapy rather than HER2-targeted therapy. Figure 6. [155]Figure 6 [156]Open in a new tab Multifaceted nature of treatment resistance in vivo (A) Overlap of resistance mechanisms was found among treatment-resistant tumors. (B) ROC curve of proteins representative of the main resistance mechanisms shows great sensitivity and specificity. To evaluate the proteome signatures we identified, we selected five proteins representative of the most prevailing resistance mechanisms: HER2, HER4, ER, IGF1R, and Kalirin. Accurate quantification of this panel of proteins by proteomics analysis was strongly associated with treatment response with high sensitivity and specificity, yielding a receiver operating characteristic (ROC) score of >0.95 ([157]Figure 6B). The improved area under the curve (AUC) for a combination of resistance mechanisms compared with the individual proteins ([158]Figure 6B) highlights the complex nature of therapy resistance in vivo and underscores the need for deep tumor profiling by system-wide quantitative analysis. Discussion In this study, we performed deep proteomics and phosphoproteomics profiling of 45 HER2^+ breast cancer tumors prior to the start of neo-adjuvant treatment with PTC-Ptz to identify molecular signatures predictive of treatment response. We show that enhanced IBC subtype classification, based on proteomics quantification of IBC biomarkers (HER2 and ER) combined with the assessment of associated biological features (such as phosphorylation abundance, expression levels of related proteins and activity of downstream signaling), improves treatment outcome prediction. We demonstrate that these observations are conserved in a second patient cohort of HER2^+ pre-treatment biopsies. Furthermore, we identify cellular mechanisms that precondition tumors to evade drug treatment: UPR-induced cellular dormancy, a metabolic switch toward OXPHOS and low levels of immune cell infiltration. These resistance mechanisms combined constitute a strong signature associated with treatment success within our patient cohort. Our study highlights the multifactorial nature of drug response in vivo and demonstrates the necessity of deep tumor profiling. We find that the current tumor classification of HER2 and ER can be improved by quantification of these IBC biomarkers combined with assessment of associated biological features. By categorizing patients as pseudo-HER2 or ER-driven using this combined approach, we were able to identify an enriched population of patients who were resistant to treatment. In contrast, we did not observe any predictive value in the expression of PR in IHC scores or proteome data (data not shown). This is in line with common clinical standards in which PR scoring is not utilized to drive treatment decisions.[159]^36 We validated our pseudo-HER2 and ER-driven findings in an independent cohort of neoadjuvant samples with matched treatment regimens in the I-SPY2 dataset. Extended analysis of the entire HER2^+ patient group within this large clinical trial further supports our observations that no pCR outcomes for patients are associated with low HER2^+ protein expression and low HER2^+ activity.[160]^37 The concept that better patient stratification can be achieved through enhanced IBC biomarker assessment is also in line with research that demonstrates the highly heterogeneous nature of tumors classified as HER2^+ by IHC/fluorescence in situ hybridization (FISH)[161]^9^,[162]^14 and research showing the subjective nature of ER classification by IHC scoring.[163]^38 It is important to note that tumor heterogeneity can contribute to discrepancies between measurements and cannot be ruled out as a factor influencing IHC and proteome differences in this study. Nevertheless, the shortcomings in the existing classification procedures hamper the effective deployment of subtype-specific therapeutic interventions. Improvements in patient stratification could directly affect clinical decision-making because suitable treatment strategies for the different IBC subtypes are readily available. Although our data and the data from Clark et al.[164]^37 indicate that improved treatment precision can be achieved through combining protein activation data with protein abundance measurements in neoadjuvant biopsies, it is still challenging to sufficiently stratify patients to a degree that warrants clinical implementation. In the current dataset, 30% of patients within the “ER-dependent” classification achieved pCR, indicating the applied treatment regimen was appropriate. To achieve stronger separation between pCR and No pCR patients using pre-treatment biomarkers, more in vivo data are needed to identify molecular signatures that correlate with patient response. Furthermore, quantification of relevant biomarkers must be translated from relative to absolute abundances for clinical implementation. This can be achieved by using multiplexed targeted proteomic analysis with spiked-in standard for absolute quantification. Our data suggest that overexpression of HER4 could compensate for HER2 inhibition, contributing to therapy resistance. In support of our findings, we observed a similar increase in HER4 expression in pre-treatment biopsies from patients who did not respond to PTT therapy in the I-SPY2 clinical trial. Nevertheless, the role of HER4 in breast cancer remains ambiguous; increased HER4 levels have been linked to a favorable disease course, yet have also been associated with poor outcomes and trastuzumab resistance.[165]^14^,[166]^39^,[167]^40^,[168]^41 Furthermore, our findings are inconsistent with previous findings in which low HER2, HER3, and HER4 levels were linked to treatment resistance against trastuzumab and pertuzumab.[169]^14 This suggests that the predictive value of HER4 expression is highly context dependent and requires further investigation. Our results also indicate a potential role of IGF1R signaling in therapy resistance. This finding is in agreement with previous research and supports efforts into the development of IGF1R inhibitors.[170]^16 The clinical benefits of these inhibitors have been very limited, however, likely because of the lack of patient stratification prior to study inclusion.[171]^42 Our data clearly indicate that only a subset of treatment-resistant tumors could benefit from IGF1R-targeted therapy ([172]Figures 2D and 2E). Furthermore, the efficacy of IGF1R inhibition is context dependent, as shown by the potential interplay between the ER and IGF1R.[173]^18^,[174]^19 This again highlights the need for precision oncology by system-wide molecular profiling to identify these therapeutic vulnerabilities and put them in their biological context on an individual patient basis. Another factor we identified to be associated with poor treatment response in our dataset is low levels of immune cell infiltration in the tumor biopsies. Consistent with previous findings in triple-negative BC patients,[175]^43 where pCR rates correlated with TIL abundance, we found that No pCR patients were significantly enriched in the low-TIL patient subset. We validated our findings in the I-SPY2 patient cohort, where we saw similar enrichment of no PCR patients within the low-TIL patient subset. Interestingly, within our low-TIL subset, we observed a significant increase in Kalirin protein expression and subsequent enrichment of Kalirin-associated signaling. Mutations in the KALRN gene have been identified as a biomarker for positive immunotherapy response, which is attributed to the loss of Kalirin function in activating DNA-damage repair.[176]^35 Accordingly, enhanced tumor mutational burden was associated with KALRN mutations in six independent patient cohorts.[177]^44 We hypothesize that the heightened Kalirin levels within the low-TIL No pCR population we observe in our dataset may contribute to treatment resistance through enhanced DNA-damage-repair mechanisms. This could enable improved resilience to chemotherapy treatment and continued immune escape by reduced neoantigen loads. Assessing the post-treatment mutational burden in these tumors would provide useful information to support this hypothesis and an important consideration for future study designs. Taken together, our results indicate that there is no singular mechanism that can be attributed to treatment resistance, but rather a panel of molecular features. From this panel we identified a treatment-resistant signature that could predict response based on the expression of five proteins (HER2, HER4, ESR, IGF1R, Kalirin). It is important to note that this signature is derived from a limited number of samples and requires further validation at the protein level in additional patient cohorts. Interestingly, we observed that in npCR tumors, the expression levels of this treatment-resistance signature were halfway between the pCR and No pCR tumors ([178]Figures 1A, 1G, and [179]S3B). In addition, we were unable to identify npCR-specific protein signatures. This indicates that the npCR tumors may harbor a mixed cell population, composed of cells that are treatment responsive and resistant. This group of patients likely benefits from the administered therapy but may require additional therapeutic intervention to exploit the therapeutic vulnerabilities of the resistant tumor cells. We think that this study highlights the feasibility and necessity of clinical proteomics in the study of drug resistance. The lack of predictive protein biomarkers for treatment outcome despite years of research is troubling and necessitates a different approach.[180]^45 A shift toward system-wide analysis of clinical samples can overcome the challenging translation of in vitro findings to patient samples and provide insights into the heavy crosstalk between signaling pathways in vivo. In this study, we show that the technological advancements made in recent years allow for great sampling depth of a standard needle biopsy. Because no additional surgery is needed, this approach could be applied to much larger cohorts of patients, paving the path to accurate patient stratification and precision oncology. Further research is needed to gain insights into in vivo tumor biology, especially in the field of phosphoproteomics. Our understanding of the biological relevance and clinical significance of phosphosites is lacking, hampering interpretation of phosphoproteomics data. To make advancements in this field, large clinical phosphoproteomics datasets are needed to correlate phosphorylation status with biological outcomes, such as our current study, where we identified a treatment-resistant signature primarily based on our proteome data, and further supported by the phosphoproteomics data. Limitations of the study We emphasize that the identified resistance mechanisms require further investigation, especially seeing the small sample size of this study. The aim of this study was to gain insights into resistance mechanisms observed in vivo. Some resistance mechanisms were identified in only a very limited number of patients (such as a switch toward OXPHOS). In addition, because all patients within the cohort received both chemotherapy and targeted therapy, we cannot delineate which treatment regimen correlates with our observed resistance signatures, or if our observations are associated with the combined effects from both PTC and TTZ/Ptz treatment together. We acknowledge that these results should be interpreted cautiously and require further study before they can be considered clinically actionable. STAR★Methods Key resources table REAGENT SOURCE IDENTIFIER Biological samples __________________________________________________________________ Human IBC tumor samples TRAIN2 clinical trial [181]https://clinicaltrials.gov/ct2/show/NCT01996267 __________________________________________________________________ Chemicals, peptides and recombinant proteins __________________________________________________________________ Sodium deoxycholate (SDC) Sigma Aldrich MFCD00064139 Tris(2-carboxyethyl) phosphine hydrochloride (TCEP) Sigma Aldrich MFCD00145469 Chloroacetamide (CAA) Sigma Aldrich MFCD0008027 TRIS Sigma Aldrich PhosSTOP Merck 4906837001 cOmplete, Mini, EDTA-free Protease Inhibitor Cocktail Merck 11836170001 TMT 11-plex reagent Thermo Scientific Cat# A34808 Trypsin Thermo Scientific Prod# 90057S LysC Wako 125–05061 Criterion XT Gels Bio-Rad 3450130 PVDF Membranes Bio-Rad 1620177 XT MOPS Bio-Rad 1610788 Pierce™ ECL Plus Western Blotting Substrate Thermo Scientific Cat#32132 __________________________________________________________________ Antibodies __________________________________________________________________ Estrogen Receptor α (D8H8) Rabbit mAb Cell Signaling #8644 HER2/ErbB2 (D8F12) XP® Rabbit mAb Cell Signaling #4290T Anti-Rabbit IgG, HRP-linked Ab Cell Signaling #7074 β-actin (14E5) Rabbit mAb Cell Signaling #4970S __________________________________________________________________ Critical commercial assays __________________________________________________________________ AssayMap Cartridge Rack, Fe(III)-NTA 5 mL Agilent Technologies Cat#G5496-60085 Kinetex 5u Evo C18 100A 150x2.1mm Phenomenex 00F-4633-AN Poroshell 120 EC-C18 2.7 micron Agilent Technologies AG699975-902 AssayMAP Cartridge Rack C18 5ul Agilent Technologies Cat# 5190-6532 Sep-Pak C18 1 cc Vac Cartridge Waters WAT023590 Bradford Protein Assay Bio-Rad 5000006 __________________________________________________________________ Deposited data __________________________________________________________________ Raw data PRIDE PXD034643 Proteomics search results PRIDE PXD034643 __________________________________________________________________ Software and algorithms __________________________________________________________________ Proteome Discoverer 2.2 Thermo Scientific OPTON-30812 GraphPad Prism 9.3.0 Graphpad Software Inc [182]https://www.graphpad.com/scientific-software/prism/ R v4.3.3 [183]https://www.r-project.org/ [184]https://www.r-project.org/ RStudio (v2022.02.1 461) RStudio, PBC [185]https://www.rstudio.com/ [186]Open in a new tab Resource availability Lead contact Further information and requests for resources and reagents should be directed to and will be fulfilled by the lead contact, Maarten Altelaar ([187]m.altelaar@uu.nl) Materials availability The study did not generate new unique reagents. Experimental model and subject details Forty-five patient biopsies of treatment-naive primary breast tumors were obtained from patients enrolled in the TRAIN-2 clinical trial ([188]NCT01996267)[189]^46 and were collected at the Netherlands Cancer Institute between 2013 and 2016. The use of biobank samples for this study was approved by the institutional review board (IRB) of the Netherlands Cancer Institute under number CFMPB672. All patients received neo-adjuvant treatment consisting of paclitaxel, carboplatin, TTZ, Ptz (PTC-Ptz). The proteomics study as presented here was performed on remaining frozen biopsies in the institute’s biobank. After careful consideration, we were able to include 45 high quality biopsies for proteomics studies. As the current study is primarily explorative in nature we did not perform a power calculation. Method details Biopsy preparation Patient pre-treatment biopsies were taken with a 14G needle and flash frozen in liquid nitrogen. Biopsies were sliced and approximately one-third of the material was allocated for this study (ca. 6–12 mg of tissue). Only biopsies with >60% tumor cells were selected for analysis. Patient information and tumor details can be found in SI [190]Table S1. Sample lysis and preparation Tissue lysis was performed in a 1% (w/v) sodium deoxycholate lysis buffer containing 10mM TCEP, 100mM TRIS, 40mM chloroacetamide, and protease inhibitor and phosphatase inhibitor tablet. Tissue slices were homogenized by multiple freeze-thaw cycles in combination with grinding by pestle in a 1.7mL sample tube in lysis buffer. Samples were then boiled for 5 min at 95°C and sonicated in a Bioruptor 300 (Diagenode) water bath for 30 min using 30 s cycles. Protein quantification was performed using Bradford Protein Assay (Bio-Rad) and 150ug of protein per sample was digested overnight with Lys-C (1:75) and trypsin (1:25) at 37°C. Samples were acidified and desalted using C18 cartridges on the AssayMap BRAVO Platform (Agilent Technologies). The TMT reference channel was made by pooling 25 μg of peptide from each sample. Samples were dried and resuspended in 50mM HEPES buffer, randomized, and then labeled with 10-plex TMT reagent (Thermo Scientific) in a 1:2 ratio (peptide: label) for 1.5 h at room temperature. TMT labeling reaction was quenched using a 5% hydroxylamine solution before samples were mixed in equal ratios to generate 5 complete TMT sets. Pooled samples were then desalted using Sep-Pac C18 cartridges (Waters), and fractionated on a high-pH reversed-phase C18 column (Kinetex 5u Evo C18 100A, 150 × 2.1mm, Phenomenex) coupled to an Agilent 1100 series HPLC over a 60 min gradient. Fractions were concatenated to 20 fractions for proteome analysis and further pooled to 10 fractions for phosphoproteome enrichment. Phosphoproteome samples were enrichened using Fe(III)-IMAC cartridges on the AssayMap BRAVO platform (Agilent Technologies) following the method described previously.[191]^47 LC-MS/MS analysis Fractionated TMT samples were analyzed by nanoLC-MS/MS on a Q Exactive HF-X mass spectrometer (Thermo Scientific) in-line with an Agilent 1290 HPLC system possessing a Reprosil pur C18 trap column (100 μm × 2 cm, 3 μm, Dr. Maisch) and a Poroshell 120 EC C18 analytical column (75 μm × 50 cm, 2.7 μm, Agilent Technologies). Samples were trapped for 5 min at a flow rate of 0.05 mL/min in 100% buffer A (0.1% FA) followed by elution with buffer B (0.1% FA, 80% ACN) at a flowrate of 300 nL/min over an LC gradient of 65 min (15%–45% B) for proteome fractions and a 95 min gradient (9%–35% B) for phosphoproteome fractions. MS settings were as follows: full MS scans (375–1500 m/z) were acquired at 60,000 resolution with an AGC target of 3e6 charges and max injection time of 20 msec. HCD MS2 spectra were generated for the top 12 precursors using 45,000 resolution, 1e5 AGC target, a max injection time of 80 msec, a fixed first mass of 120m/z, and a normalised collision energy of 32%. MS2 isolation windows were 0.7 Th for proteome samples and 1.2 Th for phosphoproteome samples. Western Blot analysis Tissue lysates containing a total of 40ug proteins were loaded on Criterion XT Gels (Bio-Rad). SDS-PAGE Electrophoresis was performed in XT-MOPS running buffer (Bio-Rad), and afterward, proteins were transferred into Immun-Blot PVDF Membranes (Bio-Rad). After blocking with 5% skim milk at room temperature for 1 h, membranes were incubated with primary antibodies at 4°C overnight. The primary antibodies against Estrogen receptor α, Her-2 and β-Actin (Cell Signaling Technology) were diluted 1:1000 with tris-buffered saline containing 1% Tween 20 (TBST). Incubation with secondary horseradish peroxidase-conjugated anti-rabbit antibody (Cell Signaling Technology) in 1:3000 dilution was performed at RT for 1h. The chemiluminescence detection was performed with Pierce ECL Plus Western Blotting Substrate (Thermo Fisher Scientific) and detected with Amersham Imager 600 (GE Healthcare, UK). Quantification and statistical analysis Preprocessing of datasets Raw data files were processed with Proteome Discover 2.2 (Thermo Scientific) using a Sequest HT search against the Swissprot human database. Results were filtered using a 1% FDR cut off at the protein and peptide level. TMT fragment ions were quantified using summed abundances with PSM filters requiring an S/N ≥ 10 and an isolation interference cut off of 35% or 50% (proteome and phosphoproteome respectively). Normalised protein and peptide abundances were extracted from PD2.2 and further scaled and analyzed using Rstudio. Data was normalised by the pooled reference, Log2 transformed and normalised by median subtraction. Phosphoproteome data was filtered to include only phosphopeptides with a class I phosphosite localization (ptmRS score >0.75). Phosphopeptides containing identical phosphorylation site localizations with different methionine oxidation sates or peptide missed cleavages were summed together to generate one quantitative value per unique phosphosite. Phosphosites quantified in peptides with different phosphorylation multiplicity states (i.e., doubly or singly phosphorylated) were not combined together and left as separate quantified values. Statistics To compare three means, a one-way ANOVA test was used (using the aov-function, followed by TukeyHSD-function in Rstudio). To compare two means, a two-sample T-test was used (using the t.test-function in Rstudio). A p value below 0.05 was regarded as statistically significant. Correlation analysis was performed by Pearson correlation using cor-function in RStudio. PTM-SEA analysis was performed using PTMsigDB.[192]^48 Pathway enrichment analysis was performed using GSEA[193]^49^,[194]^50 or Metascape using a custom background of all detected proteins in our MS analysis.[195]^51 Heatmaps were generated using pheatmap-function in Rstudio. Data was z-scored and Euclidean distance was used for clustering. ROC analysis was performed using pROC-function in RStudio, comparing No pCR and pCR patients only (near pCR patients were removed from the dataset prior to analysis). ROC curves for individual proteins were generated using roc-function, ROC curve of the combination of proteins was generated using multiclass.roc-function. Acknowledgments