The impact of tumor epithelial and microenvironmental heterogeneity on treatment responses in HER2+ breast cancer

Despite the availability of multiple human epidermal growth factor receptor 2–targeted (HER2-targeted) treatments, therapeutic resistance in HER2+ breast cancer remains a clinical challenge. Intratumor heterogeneity for HER2 and resistance-conferring mutations in the PIK3CA gene (encoding PI3K catalytic subunit α) have been investigated in response and resistance to HER2-targeting agents, while the role of divergent cellular phenotypes and tumor epithelial-stromal cell interactions is less well understood. Here, we assessed the effect of intratumor cellular genetic heterogeneity for ERBB2 (encoding HER2) copy number and PIK3CA mutation on different types of neoadjuvant HER2-targeting therapies and clinical outcome in HER2+ breast cancer. We found that the frequency of cells lacking HER2 was a better predictor of response to HER2-targeted treatment than intratumor heterogeneity. We also compared the efficacy of different therapies in the same tumor using patient-derived xenograft models of heterogeneous HER2+ breast cancer and single-cell approaches. Stromal determinants were better predictors of response than tumor epithelial cells, and we identified alveolar epithelial and fibroblastic reticular cells as well as lymphatic vessel endothelial hyaluronan receptor 1–positive (Lyve1+) macrophages as putative drivers of therapeutic resistance. Our results demonstrate that both preexisting and acquired resistance to HER2-targeting agents involve multiple mechanisms including the tumor microenvironment. Furthermore, our data suggest that intratumor heterogeneity for HER2 should be incorporated into treatment design.


INTRODUCTION
Amplification and overexpression of ERBB2 encoding for the human epidermal growth factor receptor 2 (HER2) distinguishes a subtype of breast cancers that accounts for approximately one fifth of all invasive breast cancer cases (1).Inhibition of HER2 was one of the first examples for targeted cancer therapy based on the development and use of the anti-HER2 antibody trastuzumab (2).Over the past two decades, the combination of trastuzumab with chemotherapy became a standard of care for patients with HER2 + breast cancer.Even though this targeted approach significantly improves the disease-free and overall survival of patients with HER2 + breast cancer, virtually all patients with advanced HER2 + disease will eventually develop resistance and progressive disease.Thus, to further improve treatment efficacy, numerous other HER2-targeting agents have been developed and evaluated in the clinic including various HER2 antibodies and small molecule inhibitors (SMIs) of the HER2 kinase (3).Trastuzumab and pertuzumab are two FDA-approved monoclonal antibodies that bind to the extracellular domain of HER2 and inhibit its activity while also activating the anti-tumor immune response via antibody-dependent cellular cytotoxicity (3).HER2-targeting antibodies were also used to engineer antibody-drug conjugates (ADCs) such as trastuzumab emtansine (T-DM1) (4).Upon binding to HER2, the T-DM1 is internalized into lysosomes, where it is degraded, releasing its microtubule-inhibitor payload (DM1) directly into the HER2+ cancer cell.SMIs inhibiting the tyrosine kinase activity of HER2 receptor complexes, such as lapatinib, were also shown to have some activity in a subset of patients (3).Several potential mechanisms of resistance to HER2-targeted therapy have been identified from preclinical and clinical studies.These include genetic alterations such as mutations in PIK3CA and ERBB2 leading to constitutive activation of downstream signaling pathways (5).
Upregulation of multiple other pathways, such as MET or SRC-FAK signaling, can also promote protumorigenic signaling during HER2 inhibition (6).However, despite accumulating knowledge in this area, the actual molecular changes driving resistance in human cancers have not been definitively demonstrated, and accurate predictions of the likelihood of resistance based on diagnostic biopsy profiles are not yet feasible.
A major obstacle to the effective treatment of HER2 + breast cancers is intratumor heterogeneity (ITH) for HER2 itself (7).The latest ASCO CAP guidelines support reporting of HER2 status as positive for tumors if at least 10% of the cancer cells stain positive for HER2 by immunohistochemistry (8).Thus, within each HER2 + tumor, there may be many cancer cells that lack HER2 and are genetically and functionally different from their HER2 + counterparts.We previously described a novel method, STAR-FISH, to assess cellular genetic heterogeneity for ERBB2 copy number and the PIK3CA H1047R hotspot mutation in a small cohort of HER2+ tumors subjected to neoadjuvant chemotherapy (9).Using this approach, we showed that patients with a significant increase in spatial cellular genetic heterogeneity after neoadjuvant treatment had shorter recurrence-free survival compared to patients with no change.Furthermore, two patients who received a combination of neoadjuvant chemotherapy and trastuzumab showed a decrease in spatial cellular genetic heterogeneity, yet still had poor outcome, implying the presence of pre-existing subpopulations resistant to HER2-targeted therapies.
However, how ITH impacts and is impacted by different types of HER2-targeted therapies remains unclear.In this study, we explored the effects of intratumor cellular heterogeneity for ERBB2 copy number and PIK3CA H1047R on the response to different types of HER2-targeting therapeutic strategies and changes in this heterogeneity during treatment in human breast tumors and patient-derived xenografts (PDX).

Cellular genetic heterogeneity and response to neoadjuvant HER2-targeting therapies
To assess the associations between cellular heterogeneity for ERBB2 copy number and mutant PIK3CA H1047R as well as response to treatment and changes in heterogeneity due to therapy, we performed STAR-FISH on cases from two uniformly treated cohorts of HER2 + breast cancer patients.
One cohort is from Norway (NOR) and consists of 30 cases treated with neoadjuvant chemotherapy and trastuzumab, 10 of which have both pre-and post-treatment samples; the other cohort (T-DM1) contains 16 cases treated with neoadjuvant T-DM1 and pertuzumab (7) (Supplemental Figure 1A,B, Supplemental Table 1).To assess spatial heterogeneity, we analyzed both pre-and post-treatment biopsies in the NOR cohort, whereas in the T-DM1 cohort only pretreatment samples were analyzed, with two spatially-separate biopsies per case.We utilized STAR-FISH, a combination of mutationspecific in situ PCR and fluorescence in situ hybridization which allows for detection of point mutation and copy number alterations at the single cell level in an intact archival sample (9) (Figure 1A).The signals from each fluorescent channel were quantified in individual nuclei and cells were assigned to one of five genotypes reflecting ERBB2 amplification and PIK3CA H1047R mutation status (Amp, Mut_Amp, Mut, WT_Amp, and WT; see Methods for details).The genotypes were then projected as genotype topology maps (Figure 1A).The number of cells quantified in each biopsied varied, thus, we compared relative genotype frequencies across samples and patients rather than absolute cell numbers.
First, we investigated whether the frequency of cells with ERBB2 amplification (Amp) and PIK3CA H1047R (Mut) changes during neoadjuvant treatment.We observed a significant enrichment for cells with the Mut_Amp genotype in the post-treatment as compared to the pre-treatment samples in the NOR cohort (Figure 1B), both in the entire cohort (p = 0.01, mean frequency change 0.0723) and also when individual cases were analyzed separately (Figure 1C,D and Supplemental Figure 2A,B, Supplemental Table 2).The fraction of Amp cells did not change significantly (p = 0.173, mean frequency change 0.0291), while the frequency of wild-type (WT) cells was significantly lower after treatment (p = 0.0497, mean frequency change -0.129).The spatial distribution of cells with different genotypes, based on recorded coordinates for each individual nucleus analyzed and measured by kmeans clustering (see Methods), was not significantly altered during treatment (Supplemental Figure 2C).
We then investigated whether the frequency of cells with a certain genotype in pre-treatment samples was associated with response to neoadjuvant therapy defined at the time of surgical excision as complete pathologic response (pCR) or no-pCR.We performed hierarchical clustering of patients using pre-treatment samples only.Due to the variability between different regions within the same sample (Figure 1C, Supplemental Figure 2A), we used the overall genotype frequency for each patient, rather than using each sampled region individually.We identified four major clusters (groups), with group 4 enriched in no-PCR cases.(Figure 1E).While the genotype frequencies differed significantly between the 4 groups, we detected no clear association between group 4 and the other groups (Figure 1E,F), suggesting there is no discernable difference between pCR and no-pCR groups based on pre-treatment biopsies.Both groups of samples had a similar spatial distribution of cells with different genotypes, except for WT-Amp cells, which were more dispersed in tumors with pCR (Supplemental Figure 2D).Thus, although we did observe an enrichment of Mut-Amp cells in the residual tumors, higher frequency of these cells in pre-treatment samples did not significantly impact the response to neoadjuvant therapy.
We have also explored potential associations between genotype heterogeneity and clinicopathologic information.Hormone receptor status of all tested samples did not correlate with changes in frequency of cells with distinct genotypes, except for cases with high HER2 expression having higher frequency of Mut_Amp cells, compared to cases with medium HER2 expression (p=0.02;Supplemental Figure 2E).There was also no association between the overall heterogeneity of the pretreatment tumor samples from both responders and non-responders and the presence of distant metastases (p-0.27Supplemental Figure 2F).
Next, we investigated whether the overall heterogeneity of the pretreatment tumor samples or a change in the extent of tumor heterogeneity pre-and post-treatment reflects a risk for breast cancerspecific death.While in many cases, there was a significant change in tumor heterogeneity in pre-vs.post-treatment samples, there was no significant association between the change in tumor heterogeneity and survival, a finding potentially caused by small sample size (Supplemental Figure 3A,B).There was also no significant association between the overall heterogeneity of the pretreatment tumor samples from both pCR and no-pCR cases and long-term overall survival, although patients with higher heterogeneity pretreatment appeared to have better survival (Supplemental Figure 3C).This observation was unexpected because higher diversity is generally associated with poor outcomes.
Upon further investigation, we found that higher diversity tumors have significantly more cells with ERBB2 amplification (p=0.0079) and significantly fewer cells without the ERBB2 amplification or PIK3CA mutation (p=0.0003;Supplemental Figure 3D).This suggests that the lower diversity tumors in this patient cohort consist primarily of HER2 negative cells that do not respond to HER2-targeted therapy.The higher diversity tumors, on the other hand, have larger populations of ERBB2 amplified cells and thus respond better to HER2-targeted therapy.Moreover, we found that patients with a higher level of ERBB2 amplification have better survival outcomes, although this difference was not significant likely due to the small sample size (p=0.8;Supplemental Figure 3E).
Similarly, in the T-DM1 cohort, the frequency of cells with distinct genotypes was not significantly different between patients who achieved pCR and those with no-pCR (p = 0.65 for Amp, p = 0.51 for Mut, p = 0.38 for Mut_Amp, p = 1 for WT, p = 0.57 for WT_Amp; Figure 1G and Supplemental Figure 4A,B).Similar to the NOR cohort, untreated samples from patients with pCR and no-PCR in the T-DM1 cohort had similar spatial distributions of cells with distinct genotypes and the overall pretreatment heterogeneity did not stratify patients by response to neoadjuvant therapy (p = 0.56; Supplemental Figure 4C,D).We also did not find significant associations between pretreatment genotype frequency and hormone receptor and HER2 status (Supplemental Figure 4E).Since this T-DM1 trial is still ongoing, we were not able to analyze associations with long-term outcome.

Cellular phenotypic heterogeneity and response to neoadjuvant HER2-targeting therapies
To test whether phenotypic heterogeneity could be a better predictor of response to neoadjuvant therapy than metrics based on genetic features, we performed cyclic immunofluorescence (CycIF) for 22 protein markers on consecutive slides of pretreatment samples from T-DM1 cohort (Figure 2A).This allowed us to classify individual cells within each tissue sample as tumor, stroma or immune subtype based on marker combinations (Supplemental Table 3).Tumor cells were classified according to their expression of HER2 and estrogen receptor (ER) into four categories (Figure 2B, C).
The frequency of these cell types in each case correlated with response to neoadjuvant therapy, with pCR cases having higher frequency of HER2 + cells and no-pCR cases having higher proportion of ER - HER2 -cells (Chi-squared test p = 0.022, Figure 2B-D).Moreover, pCR was also associated with increased fraction of CD8 + T cells among tumor infiltrating immune cells (p = 0.016, Figure 2D).Next, we asked if the proximity between distinct cell types could be predictive of response to neoadjuvant therapy (Figure 2E).We found that tumors with pCR are characterized by more frequent contacts between tumor cells and tumor cells of luminal HER2 + subtype (p = 0.004).Tumors with no-PCR had higher frequency of tumor cells neighboring granzyme B positive (GZM + ) macrophages (p = 0.023).
The immune cells in tumors with no-pCR also had fewer contacts with vimentin + stromal cells (nontumor, non-immune, non-endothelial cell gate, p = 0.048), potentially representing mesenchymal and more migratory phenotype, than those in pCR cases.Moreover, immune cells in no-pCR cases had significantly (p = 0.016) fewer neighboring FoxP3 + CD8 + regulatory T cells compared to tumor tissues with pCR.Thus, our results suggest that resistance to combined T-DM1 and pertuzumab treatment might be associated with impaired immune targeting and killing of tumor cells.
Using the single-cell CycIF quantitative data, we tested whether ITH for expression of HER2 and ER protein or the diversity of immune cells within a tumor could stratify the samples according to response to neoadjuvant treatment.Neither of the heterogeneity metrics showed a significant association with pCR (Figure 2F).There was also no significant correlation between HER2 expression by CycIF and ERBB2 amplification measured by STAR-FISH (Supplemental Figure 5A), and we did not find significant trends in CycIF data that was associated with the frequency of PIK3CA mutant cells (Supplemental Figure 5B).Since the T-DM1 trial is ongoing, it remains to be seen if cellular ITH for ERBB2 and PIK3CA H1047R or the phenotypic and immune microenvironment diversity is associated with long-term clinical outcome.

Variable responses to single agent HER2-targeted therapies in PDX models
Our STAR-FISH data showed that cells with PIK3CA H1047R are often present in untreated HER2 + breast tumors as minor subpopulations that increase in frequency after neoadjuvant chemotherapy with or without a HER2-targeted agent (Figure 1B), implying that mutant PIK3CA might play a role in therapeutic resistance (9).In order to better understand the effects of mutant PIK3CA on response and resistance to different types of therapies, the same tumor should be treated with different agents, which is not possible in the clinic.Thus, we tested the response of two different HER2 + PIK3CA-mutant breast cancer patient-derived xenografts (PDX) to paclitaxel and four different HER2-targeting agents (lapatinib, T-DM1, trastuzumab, and pertuzumab).PDX1 was derived from untreated HER2 + breast cancer, carrying E365K mutation in C2 PI3K-type domain of PIK3CA gene, while PDX2, with PIK3CA H1047R hotspot mutation in kinase domain, was derived from a tumor pre-treated with several rounds of trastuzumab with vinorelbine, another tubulin-binding chemotherapeutic drug (see Methods for details).Both patients received trastuzumab after tumor resection; however, their disease course differed dramatically, resulting in stable remission (PDX1) versus metastatic disease and death (PDX2).
To generate a treatment cohort, we performed mammary fat pad injections in NOG mice (n=30 for each PDX, two injections, both inguinal fat pads in each animal).Once all tumors reached over 0.5 cm diameter, the animals were randomized to six treatment groups (n=5 animals per group,) and treated for three weeks (Figure 3A and Supplemental Figure 6A).None of the treatments induced a complete tumor regression of either PDX model in this time frame.The most significant decrease in tumor volume was observed with T-DM1 treatment in both PDX models (p = 0.0001, mean fold change 2.3 for PDX1 and 1.3 for PDX2; Figure 3B,C, Supplemental Figure 6B,C, Supplemental Table 4).The most pronounced difference between the two different PDXs was their response to paclitaxel where only PDX1 showed a decrease in size during treatment (p = 0.0001).PDX1 was also significantly more sensitive to T-DM1 (Figure 3B,C).T-DM1 is a conjugate of trastuzumab with emtansine, which exerts its cytotoxic function by binding to tubulin (10).Paclitaxel is also a tubulin-binding agent, and thus the higher sensitivity of PDX1 to both T-DM1 and paclitaxel might be due to the anti-tubulin effects of these agents.PDX2 was significantly more sensitive to pertuzumab than PDX1 (p = 0.0032, mean fold change 1.2 for PDX1 and 1.5 for PDX2; Figure 3B,C) and both tumor models were refractory to trastuzumab and lapatinib treatment (p = 0.11 and 0.24 for PDX1, p = 0.1 and 0.35 for PDX2, respectively).
To assess treatment-induced changes and explore mechanisms of resistance, we analyzed the cellular and molecular profiles of the residual tumors.The histology of T-DM1-and paclitaxel-treated samples showed some differences in the PDX1 model, with T-DM1 resulting in more stroma admixed within the tumor (Figure 3D and Supplemental Figure 7A), but none of the other treatments affected tumor cellularity.Cell death and proliferation rates measured by cleaved caspase-3 and phosphohistone H3 immunofluorescence were essentially the same in all treatment groups as in untreated controls (Figure 3D, Supplemental Figure 7A-C, Supplemental Table 5).These results suggest that after three weeks of treatment, all residual tumor cells might be therapy-resistant and that PDX1 and PDX2 differ in their initial responses to certain drugs.We reasoned that the residual tumors can be classified into two groups with distinct resistance mechanisms: 1) pre-existing resistance to treatments that had no effect on tumor growth (groups treated with lapatinib, trastuzumab, and pertuzumab for PDX1 and paclitaxel, lapatinib and trastuzumab for PDX2) and 2) adaptive resistance related to regression in tumor size due to elimination of treatment-sensitive cells.Given the distinct disease history of our PDXs, our observations suggest that strong selective and adaptive mechanisms specific for each patients' tumor may determine long-term outcomes.

Genetic mechanisms of resistance to HER2-targeted therapies
One mechanism by which cancer cells can evade targeted therapy is the lack or loss of the target in tumor cells (11).HER2 + tumors often display cellular or spatial heterogeneity for ERBB2 gene amplification and overexpression (12).Thus, HER2-targeting therapies may select for cancer cells with fewer ERBB2 copies or lower levels of HER2 protein expression.Indeed, our cycIF data showed that patients with higher fraction of HER2-negative tumor cells prior to treatment are less likely to achieve pCR (Figure 2B,D).To explore if this observation could explain the relative treatment resistance of our PDX models, we performed fluorescence in situ hybridization (FISH) to assess ERBB2 copy number and immunofluorescence for HER2, phosphorylated EGFR, and estrogen receptor (ER), which have been associated with treatment responses in patients (12).We did not observe any noticeable differences in the expression of these proteins due to treatment in either PDX model, and both PDXs retained high levels of ERBB2 amplification in all treatment groups (Figure 3D, Supplemental Figure 7A-E).These results suggest that treatment resistance in these PDX models cannot be explained by the lack of drug target (i.e., HER2) expression.
Next, we performed exome sequencing of untreated and treated PDX1 and PDX2 samples to identify pre-existing or acquired mutations that could explain the differences in response to various treatments.Both PDXs were PIK3CA-mutant, and exome-seq confirmed the presence of the PIK3CA E365K mutation in PDX1 and of the H1047R mutation in PDX2 (Supplemental Figure 8A).The allelic frequency of these PIK3CA mutations was close to 1 in all tumors tested, except in PDX1 treated with T-DM1, where a decrease to 0.75 was observed in both sequenced samples.Thus, changes in the frequency of PIK3CA mutant cells could not explain the observed differences in response to HER2targeted treatments.
To investigate whether the residual tumors acquired new mutations after treatment, we compared SNVs of treated tumors to the untreated controls.We set out to identify mutations that would be present in both replicates (two independent tumors) within a treatment group compared to the untreated samples.In PDX1, we identified 10 non-synonymous mutations that were present in both sequenced tumors from each treatment group (Supplemental Table 6).Only 3 out of these 10 mutations were pathogenic according to the COSMIC database, namely ANKRD17 p.S785N, TNRC6A p.M1320T, and USP34 p.F2756L (Supplemental Figure 8B).Except for the TNRC6A mutation, no other mutations were common between PDX1 and PDX2.TNRC6A is a component of the miRNA processing machinery also involved in RNA splicing and is frequently mutated in colorectal and gastric tumors with microsatellite instability (13).TNRC6A mutations were not present in all samples; however, almost every treatment group had a single tumor sample with this alteration, suggesting it may represent a more general resistance mechanism.ANKRD17 is a mask protein that regulates the nuclear import of the YAP proteins (14) and is also a substrate of CDK2, which plays a role in DNA replication during G/S cell cycle transition (15).USP34 was shown to inhibit stemness and epithelialto-mesenchymal transition (EMT) in breast epithelial cells (16).In PDX1, the variant allele frequencies of the ANKRD17 and USP34 mutations were higher in T-DM1-treated samples compared to all other treatments and untreated tumors, while in PDX2 the same mutations were found in only one of two sequenced samples and had lower frequencies compared to PDX1 (Supplemental Figure 8B).
We did not detect any new non-synonymous mutations in T-DM1-treated PDX2 samples, suggesting that either epigenetic or microenvironmental factors may be responsible for acquisition of resistance in this model or that the tumor already acquired resistance in the patient.Residual tumors after pertuzumab treatment harbored only one new pathogenic mutation, ARID4B p.V600A.ARID4B is a subunit of the histone deacetylase-dependent SIN3 repressor complex and promotes mammary tumorigenesis and metastasis (17).This observation again points to epigenetic remodeling as a potential mechanism for acquired resistance in PDX2.Interestingly, trastuzumab-and lapatinib-treated PDX2 tumors that did not show differences in tumor growth did acquire novel mutations, many of which have not previously been reported in the COSMIC database (Supplemental Table 7).These findings suggest that different HER2-targeting therapies exert distinct evolutionary pressures on the tumors.

Intrinsic resistance to HER2-targeted therapy is linked to activation of signaling pathways
Since the differential mutational patterns and acquisition of or selection for new mutations did not fully explain the observed treatment resistance, we next performed transcriptional profiling of the PDX tumors to address the intrinsic resistance to paclitaxel and HER2-targeting agents in our two PDX models.Genes highly expressed in untreated PDX1 compared to untreated PDX2 were significantly enriched in synaptogenesis, including neuroligin, neuregulin 3, and synaptotagmin genes, Hedgehog and WNT signaling, characterized by overexpression of GLI-1, Frizzled, Hedgehog, TCF/LEF, WNT7B, as well as the ESR2 pathway, and androgen receptor and ERBB-family signaling (Figure 4A, Supplemental Table 7).These pathways are known to be involved in resistance to HER2 signaling inhibition (6).Genes highly expressed in PDX2 compared to PDX1 were enriched in cell cycle, cytoskeleton, and DNA damage checkpoint-related functions, including tubulin, myosin II, Chk2, Cyclin B, and Aurora A.
Next, we compared the gene expression profiles of tumors from all treatment groups with corresponding untreated tumors separately for PDX1 and PDX2.In PDX1 the expression patterns of every treatment were enriched in different pathways, implying distinct mechanisms of resistance (Figure 4B, Supplemental Table 7).Only leukocyte chemotaxis, interferon signaling, and kallikreinkinin system were shared in two treatment groups.High expression of kallikreins has previously been linked to tamoxifen and paclitaxel resistance in breast and ovarian cancer, respectively (18,19).
Leukocyte chemo-attractants, such as CCL5/RANTES and CXCL16, may contribute to breast cancer progression and drug resistance by recruiting monocytes and macrophages (20)(21)(22)(23).In PDX1 tumors treated with paclitaxel, where tumor growth was significantly affected, the number of differentially expressed genes was too low for pathway analysis.This finding suggests that the residual tumor represents the resistant fraction that rebounded to the original tumor state.Alternatively, it is also possible that the bulk sequencing does not fully capture more fine-grained differences within the residual tumors.We observed clearer differences in gene expression of PDX2 tumors where protein folding and cell cycle-related pathways were enriched in multiple different treatment groups (Figure 4B, right; Supplemental Table 7).Increased expression of the components of the protein folding machinery, especially chaperones such as HSP90, may promote the evolution of new heritable traits and have been previously implicated in emergence of resistance to targeted therapy in estrogen receptor-positive breast cancer (24).Thus, the overall upregulation of genes involved in protein folding and genes responsible for cell cycle progression in several treatment groups of PDX2 tumors points to the universal nature of these resistance mechanisms.
In our analysis of genetic variants that could contribute to drug resistance, we found a mutation in a splicing factor, TNRC6A, present in at least one tumor from all treatment groups, but not in the untreated controls.Therefore, we performed an alternative splicing analysis of our bulk RNA sequencing data to assess which phenotypic traits could be attributed to changes in spliced isoforms.
In PDX1, we found a fusion of GJC3-AZGP1 genes and PPP1R1B-STARD3 genes present only in T-DM1 treated samples, but no other splicing events (Supplemental Table 8).Treated PDX2 tumors had multiple splice isoform alterations (Supplemental Table 8).Two genes, CD46 and GRB7, were differentially spliced in all PDX2 residual tumors compared to untreated controls (Figure 4C,D).
Alternative splicing of CD46, a transmembrane glycoprotein involved in innate and adaptive immune response, has been shown to produce various isoforms, with different O-glycosylation and variable cytoplasmic tails (25).Exclusion of exons 7-8 and exon 13 of CD46 was shown to be non-random in activated and memory/effector T cells (26) and it may play a role in immune evasion.GRB7 is an adapter protein involved in HER2 signaling, and its splice variants have been implicated in ovarian cancer progression (27).These findings suggest that altered splicing may affect immune response regulation as well as boost HER2 downstream signaling to promote drug resistance (28).

Cancer cell subpopulations associated with acquired resistance
Since bulk RNA sequencing may obscure the presence of minor pre-existing resistant clones and does not allow for the accurate separation of tumor epithelial and stromal cells, we performed single cell RNA sequencing (scRNA-seq).Two tumors were profiled from each treatment group, with an average of 4,521 cells per sample (Supplemental Table 9).Based on our cluster tree analysis, we used the cell clusters identified using a resolution of 0.3 (Supplemental Figure 9).Clustering of cells of human origin from each untreated tumor revealed the presence of five (PDX1) or six (PDX2) distinct clusters (Supplemental Figure 10, 11A,B).To identify the changes exerted by the different treatments on the distinct cancer cell subpopulations, we analyzed the frequencies of each cell cluster in each treatment group using chi-squared tests.In PDX1, T-DM1 treatment, which reduced tumor growth, resulted in a significant decrease in the number of cluster 2 cells, characterized by high VEGFA expression (p= 3.69e-08; Supplemental Figure 10C, Supplemental Table 10).In PDX2, the main changes associated with treatments decreasing tumor growth rates were a higher frequency of cells in cluster 1 in T-DM1-treated tumors (p=3.1e-05) and cluster 1 and cluster 4 decreased cellularity in tumors after pertuzumab treatment (p = 0.0011 for cluster 1, p = 2.61e-07 for cluster 4; Supplemental Figure 10D, Supplemental Table 10).
In search for common mechanisms of acquired resistance in the two PDX models, we performed clustering followed by UMAP dimensionality reduction of the combined data from both PDXs (Figure 5A, Supplemental Table 10).The contribution of cells from the two PDXs to each cluster was only significantly different for cluster 5 (p=3.15e-26), as this cluster mostly contained cells from PDX1 (Figure 5B).The distribution of cells in clusters 0, 3, 4, and 5 was significantly influenced by treatment (Figure 5C).T-DM1 treatment, which affected the growth of both PDXs, resulted in a higher frequency of cells in cluster 3 (p=0.004)and cluster 5 (p=0.0003)compared to other treatment groups.Cluster 3 cells were also enriched in trastuzumab-treated samples (p= 0.004), yet pertuzumab-treated samples had fewer cells in this cluster (p=0.014)compared to other treatment groups.Two HER2 antibodybased treatments, trastuzumab and pertuzumab, increased the frequency of cells in cluster 0 (ptrastuzumab=1.08e-05,ppertuzumab= 1.16e-15) and decreased the frequency of cells in cluster 4 (ptrastuzumab= 1.62e-11, ppertuzumab= 1.88e-07) compared to untreated control.
To address the functional differences among clusters with differential treatment response, we performed process network analysis of genes most abundantly expressed in these clusters.We found that trastuzumab and pertuzumab treatment resistance is linked to survival of subpopulations of cancer cells able to regulate angiogenesis and EMT (cluster 0), while leading to elimination of cells with active translation, protein folding, and ER stress response linked to apoptosis (cluster 4; Figure 5D, Supplemental Table 10).Subpopulation of cells characterized by expression of DNA damage repair Cell cycle signatures can be a major confounding factor in cluster identification in single cell expression profiling experiments.Therefore, we investigated the relationship between the cell cycle and clustering and treatment effects.Cells in different phases of the cell cycle were indeed separated in the UMAP space.However, the only treatment-associated change was a decrease in S-phase cell numbers after pertuzumab treatment (p= 0.012), when both PDXs were analyzed together (Supplemental Figure 11A-C).Moreover, regressing out the cell cycle genes did not alter the clustering, suggesting that the biological processes found in our pathway analyses drive the cell cluster identity regardless of cell cycle involvement (Supplemental Figure 12).Overall, our scRNA-seq profiling of the cancer cells from residual tumors identified several subpopulations of cells differentially affected by the treatments.Yet, none of the treatments resulted in the appearance of a new population or eradication of a specific cellular subtype after three weeks of treatment, suggesting that the acquisition of resistance is a complex process where a balance between different subpopulations might influence outcome.

Stroma involvement in HER2-targeted treatment resistance
Next, we explored differences in stromal composition between the two PDX and all treatment groups.
We found that on average, 57.34% of the single cells extracted from the tumors were of murine origin (Supplemental Table 9).This set of cells enabled us to characterize the changes in stromal cell populations associated with differential treatment response.Since the clustering of cells based on scRNA-seq is driven by cell type-specific programs, we expected to identify similar cell types associated with the tumor microenvironment of the two PDX models studied.However, combined analysis of the stroma from PDX1 and PDX2 revealed that the contribution of each PDX to the 10 identified clusters is mutually exclusive (Figure 6A,B and Supplemental Figure 13).This observation and the differential responses to treatment prompted us to analyze the stroma of the two PDXs independently (Figure 6C,D,G,H).
The analysis of cell cycle signature distribution across the stromal cells revealed a heterogeneous pattern, with cells in different cell cycle stages scattered over the clusters (Supplemental Figure 14).Thus, cell identity and/or expression programs are dominant over cell cycle programs in tumor stromal cells.Unlike in tumor cells, where different treatments did not affect the cell cycle, in the case of stroma cells from both PDXs, the frequencies of cells in G1 and S phase change upon treatment.T-DM1, which was initially effective for both PDXs, depleted S phase (p= 4.88e-09) and increased G1 phase (p= 4.12e-07) stromal cell frequency, suggesting a higher degree of response within T-DM1 affected stroma.
To investigate mechanisms of acquired resistance in the stroma, we focused on clusters in which cell distribution was altered by paclitaxel and T-DM1, the treatments that affected the growth of PDX1 (Figure 6D, Supplemental Table 11).The stroma of the residual tumors after paclitaxel and T-DM1 treatment had higher frequencies of cells in cluster 8 (ppaclitaxel =1.79e-07, pTDM1= 1.61e-08).ImmGen analysis and an additional literature search for top overexpressed markers identified these cells as alveolar epithelial cells, expressing high levels of Wfdc18 (29) (Figure 6E, Supplemental Figure 13).

Functionally, these mammary alveolar epithelial cells seem to exert immune regulation based on
MetaCore process network analysis of the genes overexpressed in cluster 8 (Figure 6F, Supplemental Table 11).Thus, the development of resistance to paclitaxel and T-DM1 in this xenograft model might be regulated by a small population of alveolar epithelial cells potentially via the induction of anti-tumor immune responses.
Stromal cells within PDX2 largely classify into similar cell types identified in PDX1 (Figure 6G, Supplemental Table 11), with several additional clusters emerging.The distribution of cells within the clusters is more variable with treatments compared to PDX1 (Figure 6H), and three distinct macrophage clusters were differentially affected by treatments in PDX2.Macrophages in cluster 0 are most similar to the adipose tissue macrophages profile in ImmGen database (Supplemental Figure 12, Supplemental Table 11), and macrophages in cluster 0, 1, and 4 demonstrate distinct pathway activation (Figure 6I).The most distinct macrophage cluster, cluster 4, expresses high levels of Lyve1 (Figure 6J), a marker of specialized perivascular macrophages involved in the maintenance of arterial homeostasis (30).Lyve1 + macrophages are among the clusters depleted in residual tumors treated with T-DM1 (p= 8.81e-08) and pertuzumab (p=4.35e-07), the two most effective treatments for PDX2, compared to all other treatments and untreated controls.(Figure 6H).The stroma of these tumors also has fewer proliferating progenitors (pTDM1=0.028,ppertuzumab= 0.00476; cluster 7) and fewer cells in two clusters identified as fibroblastic reticular cells (FRC; pTDM1,5= 5.41e-08, pTDM1,8= 2.94e-05, ppertuzumab,8= 0.0151; clusters 5 and 8).FRC, normally found in lymph nodes, may exert pleiotropic immunomodulatory functions (31).In PDX2, the two different FRC populations take part in cell-matrix interactions and regulation of angiogenesis or antigen presentation, respectively (Figure 6J-M).
Altogether, our results show that specific subsets of cells in the tumor microenvironment facilitate acquired resistance to treatment.

DISCUSSION
Being able to predict which treatment would benefit a particular case of HER2 + breast cancer could impact a significant population of cancer patients, sparing them multiple rounds of treatment and decreasing the risk of relapse.Following our previous findings (9), we here set out to test whether HER2-targeted therapy could exert strong selective pressures on genetically distinct subpopulations of cancer cells.Our single-cell in situ analyses show that the frequency of cells with resistance-conferring PIK3CA mutation is higher in untreated samples of patients who did not respond to trastuzumab/chemotherapy combination, and in those patients, the subpopulations with PIK3CA mutation and ERBB2 amplification expand after treatment.These frequency changes were significant, yet our small sample size did not allow for assigning a threshold that could be used to predict the response to neoadjuvant HER2-targeting trastuzumab/chemotherapy or T-DM1/pertuzumab combinations.However, in our cyclic immunofluorescence analysis of samples from T-DM1 cohort the frequency of HER2-expressing tumor cells correlated with response to treatment.Simultaneous single cell mutation detection and accurate protein expression in the same cell in situ would be needed to integrate genetic and phenotypic heterogeneity measurements in intact tumors.It is possible that the lack of correlation between changes in genetic heterogeneity and response in our NOR cohort was due to non-uniform chemotherapy treatment.The only significant variable associated with worse outcome was the decrease in Mut_Amp cell clustering, found in only one case analyzed.This was the only tumor in which PIK3CA mut cells were more dispersed after treatment, despite having increased frequencies.This result is similar to our previous finding and suggest PIK3CA mutant cancer cells are more migratory .Our results point to the importance of comparing samples from the same tumor before and after therapy and suggest that spatial organization of resistance-conferring cells may be an important variable in predicting long-term survival.Indeed, treatment response in T-DM1 cohort was associated with spatial features, as tumors with no pathological complete response had tumor cells intermixed with activated macrophages, suggesting their resistance to immune ques.Temporal analysis of changes exerted by treatment on heterogeneous tumors will be necessary to fully capture the intricacies of tumor evolution and understand its directionality.
Lack of response to HER2-targeted treatment is a significant clinical challenge.Intratumor heterogeneity increases the chances that some subpopulations of cancer cells within a tumor already possess the resistance traits allowing them to survive therapy.Heterogeneous tumors are also more likely to adapt to the selective pressure posed by treatment and acquire new genetic and epigenetic features conferring resistance.In this study, we modeled the resistance of HER2 + PIK3CA mut breast cancer using two patient-derived xenografts.As expected, treatment with trastuzumab and lapatinib had little effect on the PDX growth.However, the pre-existing resistance in the two PDX models was different for paclitaxel and pertuzumab.PDX1 resistance to pertuzumab might be explained by the high activity of estrogen receptor pathways and WNT and Hedgehog signaling in the untreated tumor.
Similarly, paclitaxel resistance was associated with elevated expression of genes involved in cell cycle progression.The PDXs acquired resistance to T-DM1 and either paclitaxel or pertuzumab; however, the changes in the tumor cells themselves were not as pronounced as the effects of treatments on tumor stroma.
Single cell analysis of cancer cells from both PDXs revealed increased frequency of a subpopulation able to activate an immune response, stress response, cell motility, and cell cycle after T-DM1 treatment.Moreover, the acquired resistance-related changes were prevalent in the tumor stroma.In PDX1, the resistance to paclitaxel and T-DM1 was linked to higher than expected numbers of alveolar epithelial cells.These murine cells most likely constitute the luminal fraction lining the mammary ducts, as they express luminal cytokeratins.Our results showed that major pathways activated in these cells were related to regulation of immune response.This work is the first report of crosstalk between breast cancer cells and normal mammary epithelium in the development of resistance and immunomodulation.Since our study was performed in an immunocompromised animal model, it is possible that other immune cells, including B and T cells, would also participate in this crosstalk.Further studies with on-treatment biopsies or longitudinal analysis of the tumor microenvironment during treatment will be needed to address the causal relationship between these cell types in acquired resistance.
Although the same types of stromal cells were identified in PDX2, the acquired resistance to T-DM1 and pertuzumab affected frequencies of three distinct stroma cell populations.Depletion of Lyve1 + macrophages was previously shown to increase arterial stiffness and collagen deposition around blood vessels leading to arterial dysfunction (30).This may lead to inadequate blood flow, limiting further drug distribution to the tumor and creating a pro-tumorigenic hypoxic microenvironment.Fibroblastic reticular cells in lymphoid organs control migration and positioning of immune cells and support homeostasis and immune activation (32).However, in our scRNA-seq data these cells, present within the tumor microenvironment and displaying a pro-angiogenic phenotype, were depleted in T-DM1 and pertuzumab-resistant tumors.Thus, acquisition of resistance in PDX2 seems to rely on limiting the delivery of the drugs to the tumor controlled by stroma cells and upregulation of hypoxia and stress response pathways in the tumor cells.
Our study provides a detailed insight into the pre-existing and acquired resistance of HER2 + PDX models to HER2-targeted therapy.The molecular and cellular profiles associated with resistance identified in this study merit thorough comparison with transcriptome and exome profiling from clinical trials of HER2-targeted agents.As more of these techniques are incorporated into clinical practice and larger datasets become available, further studies will be possible to test whether a combination of vascular normalization and immunotherapy could prevent the onset of resistance in HER2 + breast cancer.

Human tissue samples
The Norwegian cohort was comprised of 30 cases of HER2 + breast cancer, with matched untreated biopsy and post-neoadjuvant treatment sample.The Response Evaluation Criteria In Solid Tumors (RECIST) (33) was used to score the effect of the neoadjuvant treatment, with pathological complete response (pCR) defined as no invasive tumor cells in the primary tumor region or lymph nodes after neoadjuvant treatment.Non-pCR was defined as the presence of residual invasive tumor cells in primary tumor region or lymph nodes.In this cohort, 11 patients had a complete response, 13 had a partial response 1, 2 had partial response 2 and 2 had stable disease after neoadjuvant combination of chemotherapy and trastuzumab.The details of adjuvant treatment are shown in Supplemental Figure 1.The T-DM1 cohort was comprised of untreated FFPE tumor sections from twenty patients with HER2 + breast cancer from the NCT02326974 clinical trial.For each case two different blocks with spatially distinct diagnostic biopsies were used.Pathological response to neoadjuvant treatment was reported using the Residual Cancer Burden calculator (34).Since this trial is ongoing, long term responses are not yet available.
Patient-derived xenografts: PDX1 and PDX2 were derived from human HER2 + breast cancer by Dr.

STAR-FISH
Specific-to-allele PCR-FISH for PIK3CA H1047R mutation was performed as described previously .
The staining was imaged using a Leica SP5 scanning confocal microscope.Z-stacks from at least three areas of each sample were collected and maximum projections were used to quantify signals in each individual nucleus (ImageJ macro, Supplemental Methods).Based on this quantification, each nucleus was assigned a genotype and downstream analysis was performed with R v3.5.1 (https://www.Rproject.org).For data analysis details see Supplemental Methods.
PDXs were expanded in NOG mice. 2 million cells in 50% Matrigel (BD Biosciences) in DMEM media (Gibco), total volume 50µl, were injected per mammary fat pad in a total volume 50 µl.In our experiments, we injected 30 animals with either PDX1 or PDX2, two contralateral fat pads per each animal.Tumor diameter was monitored weekly.When all tumors reached over 0.5 cm diameter, the mice were randomized into 6 treatment groups.Paclitaxel (i.p.) T-DM1, trastuzumab and pertuzumab (i.v) were delivered once weekly, while lapatinib was dosed by gavage 5 days a week (Supplemental Figure 4).Tumor diameter was measured on treatment weekly and the experiment was terminated after 3 weeks of treatment.At the endpoint, tumors were extracted, weight and divided into 3 parts.
One was immediately processed into single cell suspension (as described before (35)), one was snap frozen and one was fixed in buffered formalin and processed into a paraffin block.

Immunofluorescence staining
Formalin-fixed paraffin embedded tissues were baked overnight at 65°C, deparaffinized and subject to antigen retrieval in citrate buffer (pH 6; Invitrogen) or Target Retrieval solution (pH 9; Dako) for 20 min in a steamer.After blocking with 10% goat serum in PBST (0.05% Tween-20 in PBS), the samples were incubated with primary antibody in 5% goat serum in PBST overnight at 4ºC and washed three times for 5 min in PBST.Fluorescent dye-conjugated secondary antibodies incubation was held for 45 min at room temperature.The details regarding antibodies used are reported in Supplemental Table 12.

Cyclic immunofluorescence
Formalin-fixed paraffin-embedded (FFPE) tissues were mounted on adhesive slides and baked overnight at 55 °C and an additional 30 minutes at 65 °C.Tissues were deparaffinized with xylene and rehydrated with graded ethanol baths.Two step antigen retrieval was performed in the Decloaking Chamber (Biocare Medical) using the following settings: set point 1 (SP1), 125 °C, 30 seconds; SP2: 90 °C, 30 seconds; SP limit: 10 °C.Slides were further incubated in hot Target Retrieval Solution, pH 9 (Agilent, S236784-2) for 15 minutes.Slides were then washed in diH2O and once for 5 minutes in 1x phosphate buffered saline (PBS), pH 7.4.Sections were blocked in 10% normal goat serum (NGS, Vector S-1000), 1% bovine serum albumin (BSA, Sigma A7906) in PBS for 30 minutes at 20 °C in a humid chamber, followed by PBS washes.Primary antibodies (Supplemental Table 12) were diluted in 5% NGS, 1% BSA in 1x PBS and applied overnight at 4° C in a humid chamber, covered with plastic coverslips (Bio-Rad, SLF0601).Following overnight incubation, tissues were washed 3 x 10 min in 1x PBS.Coverslips were mounted in Slowfade Gold plus DAPI mounting media (Life Technologies, S36938).For details of image acquisition and analysis see Supplemental Methods.For cell type determination and composition analysis, single cell mean intensity was used and threshold was assessed to generate a binary data set.Binary calls for each marker were combined into a gating strategy to define final cell types (Supplemental Table 3).For spatial analysis, nuclear centroids were used to calculate distances between cells.For each cell type A-cell type B pair, the number of cell type B within 75 micrometers of each single cell of type A was counted, and the mean number of cell type B proximal to cell type A was calculated for each case.For heterogeneity analysis, Shannon entropy of each case was calculated based on single cell phenotype frequency of tumor receptor status expression (HER2+, ER+, HER2+/ER+, HER2-/ER-) or immune cell types (B cell, CD4 + T cell, CD8 + T cell and CD68 + Macrophage).

Exome sequencing and analysis
Two independent tumors from each treatment group were used.DNA was extracted with Qiagen Blood and Tissue kit and shipped to Novogene for library construction and 100x coverage exome sequencing.
Mutec1 (36) was used to call single nucleotide variants (SNVs), and Strelka (37) and Mutec2 (38) were used to call insertions and deletions.This analysis was done in two ways: (1) a panel of normal samples from TCGA were used in place of matched controls to estimate absolute mutation levels in untreated and treated samples, and (2) the untreated samples were used as the matched controls for variant calling in the treatment samples to estimate relative changes compared to the untreated samples.
Oncotator (39) and Variant Effect Predictor (VEP) were used to annotate SNVs and INDELs.Variants appearing in non-cancer variant annotations (e.g., dbSNP and 1,000 Genomes) were filtered out from downstream analysis.Finally, variants were filtered using Orientation Bias Filters , the MAF Panel of Normals Filter and the BLAT Realignment Filter (40) The GATK3 CNV (38) pipeline was used to call copy number variants (CNVs), where the untreated samples were used in place of matched controls to call CNVs in the treated samples.R v3.5.1 (https://www.R-project.org) and Bioconductor v3.8 (41) was used for all downstream analysis.The list of CNVs was filtered to only include those with at least 100 reads.IRanges v2.16.0 (42) was used to identify CNVs with overlapping regions between the two replicates from each treatment group.CNVs with no overlaps in the other replicate were filtered out.

RNA-seq and analysis
For bulk RNA sequencing, RNA was extracted from a snap frozen tissue with Qiagen Blood and Tissue kit.Libraries were prepared using Illumina kit and sequenced according to standard protocol.
Demultiplexing and alignment of data was performed as previously described (43).To remove mouse contamination, reads that mapped uniquely to the human genome were kept for further analysis.A subset of these also mapped uniquely to the mouse genome, for which a threshold was determined for the minimum difference in alignment scores between human and mouse.The thresholds were determined, for each sample separately, by repeatedly taking random samples of the reads that mapped uniquely to both genomes, specified by varying alignment score differences, and statistically comparing sample alignment scores between human and mouse genomes using a two-tailed t-test.
The minimum difference by which samples confidently displayed p-values ≤ 0.01 defined the threshold, beyond which reads were kept.Read counts for individual genes were generated using the htseq-count script of the HTSeq framework (version 0.6.1p1)(44) using modified parameters (--stranded no) and the hg19/mm10 refGene annotation file available at the UCSC Genome .The DESeq2 (45) [R package version 1.16.1 was used to generate differential expression gene lists, with a fold-change value ≥2 and padj value ≤0.05.Pathway enrichment analysis was performed using MetaCore from Clarivate Analytics.Alternative splicing analysis was performed using LeafCutter v0.2.9 (46).GENCODE v19 (47) (GRCh37) was used to annotate exons and splice junctions.

Single-cell RNA-seq and data analysis
Single cell suspension from two independent tumors were used from each xenograft and treatment group.Cell and library preparation for single cell RNAseq was performed according to 10xGenomics Chromium v2 protocol, targeting 2,000 cells per sample.Cell Ranger v3.0.2 (48) was used to preprocess the single cell RNAseq output.`cellranger mkfastq` was used to demultiplex the base call files into FASTQ files; `cellranger count` was used to align and filter the FASTQ files individually, aligning to both the human (hg19) and mouse (mm10) genomes, and to generate feature-barcode matrices; `cellranger aggr` was used to normalize the separate runs to the same sequencing depth and aggregate the samples into one combined feature-barcode matrix.Altogether this resulted in separate feature-barcode matrices for stroma and tumor cells for each PDX model and replicate.R v3.5.1 (https://www.R-project.org) and Seurat v3 (49) were used for all downstream analysis.The data was filtered to include cells with at least 500 expressed genes and with a maximum of 25% of expressed genes coming from mitochondrial genes.The number of cells before and after filtering are shown in Supplemental Table 9. Cell cycle scores were assigned to each cell based on their expression of G2/M and S phase markers using `CellCycleScoring`, using the marker gene list `cc.genes.updated.2019`as provided by Seurat.Data was normalized using `SCTransform`.The percent of mitochondrial RNA per cell was regressed out during normalization and default values were used for all other parameters.Samples were integrated using `SelectIntegrationFeatures`, `PrepSCTransform`, `FindIntegrationAnchors`, and `IntegrateData`.The number of features, `nfeatures`, was set to 3000 and the default values were used for all other parameters.The resulting normalized data was used for all downstream analyses.`RunPCA` was used for PCA dimensionality reduction and `RunUMAP` was used for UMAP dimensionality reduction.The first 10 principal components were used for both tSNE and UMAP.`FindNeighbors` was used to construct a shared nearest neighbor (SNN) graph, with k = 20.`FindClusters` was used to identify clusters based on the constructed SNN graph.A range of resolution values from 0.1 to 1 were tested to determine the appropriate number of clusters.Clustree v0.4.3 (50) was used to visualize the resulting clustering tree and investigate how clusters at different resolutions are related to each other.`FindAllMarkers` was used to find the differentially expressed genes for each of the clusters.The resulting gene lists were filtered to include only genes with adjusted p-value less than 0.05.Pathway enrichment analysis was performed using MetaCore from Clarivate Analytics.The associations between treatment and cluster, cell cycle phase and cluster, and cell cycle phase and treatment were visualized using the `mosaic` function in vcd v1.4.4 (51).The size of each square in the mosaic plot represents the number of cells that fall into that respective category.The colors represent the Pearson residuals for each category in the contingency table, where a large value suggests greater frequency and a small value suggests smaller frequency than would be expected under the null hypothesis of independence.These Pearson residuals were used to compute separate p-values for each category.The p-value shown on the plot is computed using a chi-squared test of independence and represents the significance of the overall association between the rows and columns.

genes (cluster 3 )
Figure 6C), the scRNA-seq data suggest an overrepresentation of cycling cells in T-DM1-resistant

Figure 3 .Figure 4 .
Figure 3. Single-agent HER2-tageted treatment effects on two HER2 + PIK3CA mut patient-derived xenografts.(A) Schematic of experimental design.(B) Treatment response for PDX1 and PDX2.Waterfall plots show percent change in diameter for each tumor; n=10 tumors and 5 animals per treatment group.(C) Tumor weight at the experimental endpoint.P values in (B,C) indicate statistical significance based on unpaired two-tailed Student's t-tests.(D) Representative images of the histology (hematoxylin-eosin staining; upper panels) and immunofluorescence staining for apoptosis marker (cleaved caspase-3), cell proliferation marker (phospho-histone H3) and treatment targets (HER2 and phospho-EGFR proteins).Scale bars, 100µm.Staining was repeated twice with similar results.

Figure 5 . 16 CFigure 6 .
Figure 5. Single-cell RNAseq analysis of two PDXs identifies clusters associated with differential drug response.(A) Combined analysis of tumor cells from n=2 samples per each PDX and each condition (total 24 samples, average cell number per sample = 2355).UMAP plots colored by cluster (top panel), PDX (middle panel) and treatment (bottom panel).(B) Cell distribution among clusters based on PDX from which they were derived.Red color indicates lower than expected frequency, blue -higher than expected.P value of  !test is shown.(C) Cell distribution in different clusters based on treatment.Red color indicates lower than expected frequency, blue -higher than expected.P value of  !test is shown.(D) MetaCore Pathway enrichment analysis for genes enriched in clusters differentially affected by distinct treatments (Cluster 0, 3, 4, and 5).Circle size corresponds to number of genes found per pathway.