Characterization of breast cancer cells and breast cancer cell-derived exosomes

We analyzed BCCs and BCDEs from the MDA-MB-453 (Her2), MCF7 (LA), ZR-75 (LB), and MDA-MB-231 (TN) cell lines to identify potential diagnostic/prognostic biomarkers that may help to stratify patients and guide therapy choice (Fig. 1A describes the overall workflow). We isolated BCDEs by ultracentrifugation of BCC supernatants cultured in exosome-depleted media. We confirmed the receptor profile of BCCs by Western blotting and evaluated the presence of the CD63, CD9, and CD81 tetraspanins (TSPANs) (given their function in exosome-biogenesis and release) and other exosome-related proteins (Hsc70 and Alix) in BCCs and BCDE extracts. The absence of clathrin and calnexin from BCDE samples indicated the lack of contamination with extracellular vesicle populations derived from the plasma membrane or endoplasmic reticulum (Fig. 1B). We determined exosome size and concentration via nanoparticle tracking analysis (NTA) (Fig. 1C), finding that BCDEs possessed sizes of 110–125 nm and that BCCs exhibited specific exosome production/secretion profiles, with LA and LB BCCs giving rise to the most significant number of exosomes. We confirmed the significantly elevated level of BCDEs isolated from LA and LB BCCs found by NTA and the presence of the well-described cup-shaped BCDEs by transmission electron microscopy (TEM) (26–27) (Fig. 1D-E).

Differential protein expression analysis in BCCs and BCDEs

We preprocessed and normalized proteomic data from BCCs and BCDEs to ensure comparability; overall, we quantified 1668 (BCCs) and 389 (BCDEs) UniProtKB/SwissProt proteins. The exploratory data analysis by clustering and principal component analysis (PCA) on both data sets discounted abnormal sample behavior (Supplementary Figure S1). We compared BCCs and BCDEs from each subtype against the remaining three subtypes combined: TN vs. (Her2 + LA + LB), Her2 vs. (TN + LA + LB), LA vs. (TN + Her2 + LB) and LB vs. (TN + Her2 + LA), which will be referred to as TN, Her2, LA, and LB, respectively, for simplicity. Table S1 reports the differentially-expressed proteins (DEPs) in each BCC and BCDE comparison. Overall, BCCs displayed a more significant number of dysregulated proteins, although BCDEs exhibit more significant levels of alterations. We observed the most severe dysregulation (in terms of number) in the TN (265 DEPs) and Her2 comparisons (121 DEPs) in BCCs and the TN (109 DEPs) and Her2 (118 DEPs) comparisons in BCDEs.

We represented significant DEPs from each subtype comparison for BCCs and BCDEs and considered their direction of change (upregulation/downregulation) as shown by Venn diagrams (Fig. 1F). In general, we observed a very low level of intersections among subtypes, with the most common intersection being between LA and the remaining subtypes. Our findings also suggested that most DEPs displayed subtype-specificity in BCCs and BCDEs; the more significant number of specific DEPs makes this especially true for the TN subtype.

We also explored proteomic correlation levels between BCCs and BCDEs through the statistics of the 239 commonly measured proteins. This analysis demonstrated significant positive correlations between parental BCCs and BCDE samples and significant negative correlations between non-parental BCCs and BCDE samples, except for LA and LB (Fig. 1G; Tables S2 and 3). This result suggests that BCDE protein cargo reflects the proteomic profile of their cell of origin and provides evidence for exosomes as a potential source of non-invasive breast cancer biomarkers.

Fig. 1

Characterization of BCCs and BCDEs and differential protein expression analysis. (A) Characterization workflow. (B) Western blot of BCC extracts confirming characteristic receptor expression and BCC/BCDE extracts for exosome marker expression (Clathrin/Calnexin as quality control for contamination with other extracellular vesicle types). (C) The concentration and distribution of BCDEs were analyzed with the NanoSight NS300 (vertical axis, x107 particles/ml) and by size (horizontal axis, in nanometers [nm]). (D) Quantitative analysis of BCDEs. (E) Representative TEM images of BCDEs. Scale bar = 200 nm. (F) Intersection of significantly dysregulated proteins by subtype in BCCs and BCDEs. Venn diagrams demonstrate the number of commonly significantly upregulated (“UP”) and downregulated (“DOWN”) proteins in intersections for subtype comparisons. (G) Statistical correlations between subtype comparisons in BCCs (columns) and BCDEs (rows). Correlations among parental BCCs and BCDEs are positive, while non-parental relations present near-zero or negative values. Significant correlations are depicted with a red star. n = 3 in each experiment; error bars represent mean ± SEM, *p < 0.01, **p < 0.001, ***p < 0.0001; one-way ANOVA with post-Bonferroni’s Multiple Comparison Test

Definition of a BCC protein signature with subtype discriminatory potential

We defined a 33-protein signature of BCCs by selecting each comparison’s top ten DEPs and removing duplicates (Table S2). Figure 2A depicts the expression level of said proteins, highlighting protein-specific profiles. Of note, 24 proteins presented a pattern of specific up- or down-regulation in one of the different subtypes (bold in Table S2). Proteins specifically overexpressed in (i) Her2 BCCs mainly associate with metabolism and energy production (e.g., KRT17 and CKB), (ii) LA and LB BCCs play roles in tumor growth and angiogenesis (e.g., GFRA1 and CTSD) and (iii) TN BCCs function in cytoskeletal remodeling, migration, epithelial-mesenchymal transition, and chemoresistance. Of note, specific signature proteins do not present subtype specificity; instead, they display a gradual change across subtypes (e.g., SLC25A22, PDLIM1, CAPN2, FASN, and REEP6), while others displayed specificity to both luminal subtypes (e.g., KRT8, KRT18, HSPB1, and EPHX1).

The 33-protein signature can stratify cancer subtypes, as seen in the sample classification after the unsupervised hierarchical clustering performed (Fig. 2B) and the inter-subtype separation in the PCA plot (Fig. 2C). Both approaches demonstrated the significant difference between TN and the remaining subtypes and the relative similarity between the LA and LB subtypes.

Using STRING, we explored the BCC signature protein-protein interaction (PPI) network and functional profile [28]. The resulting PPI network displayed significantly high levels of connections (p = 6.88e-15, Fig. 2D) and presented a clear structure with a prominent central cluster comprising two-thirds of the 33 proteins and multiple isolated groups of one to three proteins. Functional enrichment revealed twenty-one biological processes (BP), six molecular functions (MF), and thirty-three cellular component (CC) enriched GO terms. These included BP GO functions Cytoskeleton organization (GO: 0007010), Cornification (GO:0070268), Oxidation-reduction process (GO:0055114), Cellular modified amino acid metabolic process (GO:0006575), and CC GO function Invadopodium (GO:0071437), which we highlighted given the interesting pattern of upregulation/downregulation of their related proteins as this provides subtype-specific functional detail.

Her2-specific downregulated BCC proteins related to functions such as Cytoskeleton organization, Cornification, and Invadopodium and upregulated proteins related to Oxidation-reduction and Cellular modified amino acid metabolic processes. LA-specific downregulated BCC proteins related to Cytoskeleton organization, Cellular modified amino acid metabolic process, and Invadopodium, and upregulated proteins related to Cornification, while LB upregulated proteins related to Cornification and downregulated proteins related to Cytoskeleton organization and Invadopodium. Finally, TN-specific upregulated BCC proteins related to Cytoskeleton organization and Invadopodium, and downregulated proteins related to Cornification and Cellular modified amino acid metabolic process. Interestingly, keratins (KRTs) represent currently applied cancer diagnosis biomarkers: TPS (KRT18), TPACYK (KRT8/18), and CYFRA 21 − 1 (KRT19) [28,29,30,31]; here, they displayed common behavior, independent of the remaining proteins related to Cytoskeleton organization and opposite to KRT7 (Fig. 2E).

In summary, a 33-protein signature of BCCs supported the clear discrimination of breast cancer subtypes and, more importantly, provided subtype-specific functional detail that could explain the differential prognosis of each breast cancer subtype.

Fig. 2

A 33-protein Signature of BCCs Discriminates Subtypes and Provides Subtype-specific Functional Detail. (A) Barplots reporting the expression levels of the 33 BCC signature proteins (Her2, LA, LB, and TN representing MDA-MB-453, MCF7, ZR-75, and MDA-MB-231). Boxes highlight upregulated (red) and downregulated (green) proteins for each BCC. (B) Heatmap of the expression levels of the 33 BCC signature proteins, with non-supervised clustering of the samples colored by the cell-of-origin. (C) The first and second components of a PCA of samples considering only the protein signature, colored by the cell-of-origin, with ellipse plots with a confidence level of 0.9. (D) PPI analysis of the BCC protein signature. Proteins colored by their association to GO functions Cytoskeleton organization (GO: 0007010), Cornification (GO:0070268), Oxidation-reduction process (GO:0055114), Cellular modified amino acid metabolic process (GO:0006575) and Invadopodium (GO:0071437). (E) Upregulated and downregulated proteins are depicted in red and blue for the specific Her2, LA, LB, and TN networks. For further information on the network, please visit https://version-11-5.string-db.org/cgi/network?networkId=bT6IRAwsY9iC

Definition of a BCDE protein signature with subtype discriminatory potential

We next defined a BCDE protein signature in the same manner by selecting the top ten DEPs in each subtype comparison and removing duplicates, which also comprised 33 proteins (Fig. 3A). Table S3 summarizes the upregulated and downregulated proteins in BCDE samples, highlighting their subtype specificity. In general, we observed a general downregulation of BCDE signature proteins in TN (25 proteins under the median) and Her2 (20 proteins under the median) and an upregulation in LA (22 proteins above the median) and LB (19 proteins above the median). Similar to the BBC signature, the BCDE signature could stratify cancer subtypes, as evidenced in sample classification after unsupervised hierarchical clustering (Fig. 3B) and inter-subtype separation in the PCA plot (Fig. 3C). Importantly, these findings suggest that the protein cargo of BCDEs partly reflects their cell of origin.

As before, we explored the BCDE signature PPI network and functional profile using STRING [28]. The resulting PPI network displayed significantly elevated levels of connections (p < 1.0e-16, Supplementary Figure S2) and possessed a similar structure to the BCC network, with a prominent central cluster comprising two-thirds of the total proteins and multiple isolated clusters of one to two proteins. Functional enrichment of signature proteins with STRING revealed 14 BCC-enriched GO terms in BCDEs; however, PPI analysis failed to reveal associations between protein clusters and biological functions in this case. Supplementary Figs. S2B-E report the different patterns of dysregulation in Her2, the less differentiated patterns in LA and LB, and the general downregulation of signature proteins in TN, respectively.

We next explored the efficacy of our protein signatures to stratify breast cancer subtypes through The Cancer Genome Atlas (TCGA) and The Cancer Proteome Atlas (TCPA). PCAs revealed that both BCC and BCDE protein signatures discriminated TN from the LA and LB subtypes in the TCGA dataset but had more difficulty separating Her2 and could not discriminate between LA and LB (Fig. 3D). Interestingly, the BCDE protein signature represented the most powerful means of separating subtypes, especially for the TN subtype, providing evidence for the potential of exosomes as a source of subtype biomarkers in liquid biopsies. In the TCPA dataset, BCC and BCDE signatures could discriminate the TN subtype from the remaining subtypes (Fig. 3D). While further studies will be needed to corroborate our findings, these preliminary results offer an interesting starting point when considering liquid biopsies as a tool for TN breast cancer diagnosis.

Fig. 3

Validation of a 33-protein signature in BCDEs. (A) Barplots reporting the expression level of the 33 BCDE signature proteins (Her2, LA, LB, and TN representing MDA-MB-453, MCF7, ZR-75, and MDA-MB-231). Boxes highlight upregulated (red) and downregulated (green) proteins for each BCDE. (B) Heatmap of expression levels of the 33 BCDE signature proteins, with non-supervised clustering of the samples colored by the cell-of-origin. (C) The first and second components of a PCA of samples considering only the protein signature, colored by the cell-of-origin, with ellipse plots with a confidence level of 0.9. (D) PCA plots of samples in the TCGA (left) and TCPA (right) datasets for the BCC and BCDE signatures proteins. For both validation steps, Her2, LA (LumA), LB (LumB), and TN (Basal) represent MDA-MB-453, MCF7, ZR-75, and MDA-MB-231

Pathway analysis reveals functional differences in BCC subtypes

We next performed pathway inference analysis (PIA) with the Hipathia package [32,33,34,35] in BCCs to explore dysregulated pathways and subpathways for each subtype from the Kyoto Encyclopedia of Genes and Genomes (KEGG) signaling pathways (34–35). Hipathia divides KEGG pathways into subpathways, which are defined as all the nodes in a path ending in a particular final node (also called the effector protein). A subpathway exists for each final node within a KEGG signaling pathway; thus, Hipathia infers the activation level for each subpathway separately as an independent biological entity. After analyzing all available subpathways, we discovered 234 (Her2), 38 (LA), 26 (LB), and 146 (TN) altered subpathways in BCCs (Table S4).

A summary of the number of altered subpathways within each pathway provided potentially interesting differences between BCC subtypes. Figure 4A depicts the top altered pathways, which support the definition of two pattern types: (i) pathways altered mainly in one subtype (e.g., cAMP, HIF-1, and NF-kB signaling pathways in Her2; AGE-RAGE signaling pathway in diabetic complications in LB; or Chemokine signaling pathway in TN) and (ii) pathways altered at the same time in different subtypes (e.g., Toll-like receptor signaling pathway, Rap1 signaling pathway, or Proteoglycans in cancer).

The subpathways associated with the second group also presented with subtype-specific dysregulation patterns for BCCs. For instance, the Rap1 signaling pathway (Fig. 4B), which promotes tumor invasion and metastasis in breast cancer [36], includes seven dysregulated subpathways in at least one subtype (identified in Fig. 4C). The subroute ending in the THBS1 effector protein (the THBS1 subroute, in blue), which relates to defective angiogenesis in the pathway (Fig. 4C), displayed overactivation in Her2 and TN and underactivation in LA and LB. The ITGA2B subroute (lilac), which functions in cell adhesion and migration, displayed a similar pattern. Cell adhesion and migration also associated with the ACTB, PRKC1/PARD6A/PARD3, and ITGAL subpathways (pink, red, and orange, respectively), which displayed deactivation in Her2. Additionally, the ITGAL subroute displayed deactivation in LA and activation in TN. The MAPK1 subroute (yellow), related to proliferation and survival, displayed overactivation in Her2, while the RAC1 subroute (dark lilac), related to SMG crosstalk, displayed overactivation in TN. We observed similar differential patterns to the Rap1 signaling pathway in additional pathways, including the Proteoglycans in cancer, Pathways in cancer, Hepatitis B, and Toll-like receptor signaling pathways (Fig. 4D-G, with subpathways of Fig. 4D and E identified in Supplementary Figs. S3 and 4). Our results indicate that although the Rap1 signaling pathway becomes systematically altered in breast cancer, each subtype possesses specific alterations in the pathway (activating or deactivating specific subpathways) and causes different potential functional outcomes. This knowledge may support the development of tailored treatments that respond to each patient’s specific requirements.

Fig. 4

Pathway analysis reveals altered functional outcomes according to BCC subtype. (A) Number of subpathways altered in the top altered KEGG pathways in each BCC (Her2, LA, LB, and TN representing MDA-MB-453, MCF7, ZR-75, and MDA-MB-231). Framed pathways are further represented in B and D-G. (B) The dotplot represents the altered subpathways’ statistical value within the Rap1 signaling KEGG pathway. (C) Representative example, detail of the Rap1 signaling KEGG pathway in each BCC. Specific parts of the altered subpathways are colored following the column Subpathway in (B). Common components are colored in light gray. For further, visit https://www.genome.jp/pathway/hsa04015. (D-G) Dot plots representing the statistical value of altered subpathways within the (D) Proteoglycans in Cancer, (E) Pathways in cancer, (F) Hepatitis, and (G) Toll-like receptor signaling pathway KEGG pathways in each BCC subtype

Functional impact of BCDE cargo in receptor cell pathways

We hypothesized that the functional characterization of BCDEs should be interpreted as the functional change that BCDEs induce in receptor cells; therefore, we analyzed altered pathways in BCDEs with Hipathia analogously to BCCs. We interpreted the PIA in BCDEs as the impact of BCDE protein cargo on signaling pathways in target cells, finding 49 (Her2), 47 (LA), 2 (LB), and 133 (TN) subpathways altered in BCDEs (Table S4). Summarizing the number of altered subpathways within each KEGG pathway, we found considerable differences among BCDEs (Fig. 5A), with specific significant alterations to Ovarian steroidogenesis, the Glucagon signaling pathway, and Adrenergic signaling in cardiomyocytes in TN BCDEs, Prostate cancer in Her2 BCDEs, and Retrograde endocannabinoid signaling in LA BCDEs. We also discovered KEGG pathways commonly altered in all BCDEs, such as Adherens Junction (Fig. 5B), whose altered subpathways are characterized in Fig. 5C. We also found the ErbB signaling pathway, Proteoglycans in cancer, and Pathways in cancer KEGG pathways altered in Her2 and TN BCDEs (Fig. 5D and F).

Interestingly, the Her2 and TN subtypes present with overexpressed pathways whose effector proteins relate to migration and invasion processes. For instance, the Adherens junction subpathway with SNAI2 effector protein (plays an essential role in promoting invasion and metastasis) and the subpathway from Proteoglycans in cancer ending in HOXD10 (related to cell migration and invasion) become upregulated in Her2 (Supplementary Fig. S5). Furthermore, the activity of specific proteins (PXN, PTK2, PRKCA, NUDT16L1, HSPB2, and ACTB) through the Proteoglycans in cancer pathway (related to migration and invasion) became increased in TN. The LA and LB subtypes present with the overactivation of subpathways concerning the effector protein CTNND1 (orange and red in Fig. 5B), which functions in angiogenesis and metastasis repression.

Fig. 5

Pathway analysis reveals altered functional outcomes according to BCDE subtypes. (A) The number of subpathways in the top ten altered KEGG pathways in each BCDE (Her2, LA, LB, and TN representing MDA-MB-453, MCF7, ZR-75, and MDA-MB-231). Framed pathways are further represented in B and D-E. (B) The dotplot represents the altered subpathways’ statistical value within the Adherens junction KEGG pathway. (C) Representative example, detail of KEGG pathway Adherens junction. Specific parts of the altered subpathways are colored following the column Subroute in B. All nodes in subpathway CTNND1** (red) are included in subroute LEF1 CTNNB1 (orange) except for pl20ctn, representing node CTNND1**. For further details, visit https://www.genome.jp/pathway/hsa04520. (D-F) Dot plots representing the statistical value of the altered subpathways within the (D) Proteoglycans in Cancer, (E) ErbB signaling pathway, and (F) Pathways in Cancer KEGG pathways in each analyzed BCDE

BCDEs differentially influence breast cancer progression as a reflection of the cell of origin

We next conducted two experiments to corroborate our hypothesis that the BCDE protein cargo functionally impacts target cells. As seen in Fig. 5, the functional results suggested that TN and Her2 BCDEs impact migration and invasion, while LA and LB BCDEs mainly impact the vascular endothelium and angiogenesis events during cancer progression; thus, we performed tube formation and scratch wound healing assays in the presence of BCDEs from specific subtypes to evaluate their potential to alter pathways according to their cargo.

Phase-contrast microscopy images from tube formation assays using human umbilical vein endothelial cells (HUVECs) revealed that LA and LB BCDEs contributed to more complete and complex nets with a more significant number of tube-like structures than TN and Her2 BCDEs or PBS (as vehicle control treatment) (Fig. 6A). The quantification of tubes, loops, and branching points (as representative measures of angiogenesis) demonstrated that LB BCDEs significantly increased the ability of HUVECs to form these structures compared to PBS (Fig. 6B-D). LA BCDEs significantly increased loop numbers compared to PBS (Fig. 6C) while LB BCDEs significantly increased loop numbers compared to TN BCDEs (Fig. 6C) and tube-like structure formation and branching points compared to TN and Her2 BCDEs (Fig. 6B and D).

In the scratch wound healing assay, we exposed fibroblasts to BCDEs and monitored scratch closure over 24 h as a proxy for migration and invasion. In this case, TN (but not Her2, LA, or LB) BCDEs significantly increased fibroblast migration compared to PBS and LA BCDEs (Fig. 6E).

In summary and in agreement with the previous pathways subroute analysis (Fig. 5), the tube formation assay results suggest that LA and LB BCDEs have a more significant impact on angiogenesis than TN and Her2 BCDEs, with LB BCDEs having the most significant impact. The scratch wound healing assay results suggest that TN BCDEs have a more pronounced invasive and metastatic potential. Together, these results support the proposed interpretation of the functional results of the PIA analysis as the different functional impact of BCDE subtypes in recipient cells according to their cell of origin.

Fig. 6

Differential contribution of BCDEs to migration and angiogenesis according to subtype. (A) Bright-field images of tube formation assays. HUVECs were seeded (4.0 × 104) on 48-well plates coated with Matrigel and incubated with PBS or BCDEs (5 µg). (B-D) Quantitative analysis representative of (B) tubes, (C) loops, and (D) branching points formed by HUVECs cultured on Matrigel-coated plates after a 12 h incubation with BCDEs. (E) Quantitative analysis of wound healing assays using fibroblasts (expressed as % fibroblast migration). Fibroblasts were treated with 5 µg BCDEs and monitored for 24 h; fibroblasts were treated with PBS as a control. Bars represent mean ± SEM, **p < 0.01, ***p < 0.001; one-way ANOVA with the post, “Bonferroni’s Multiple Comparison Test of three different experiments, n = 3 wells in each experiment

ExoCellBC web tool

We have made all data and results available within the ExoCellBC web tool (https://bioinfo.cipf.es/metafun-ExoCellBC), which is freely accessible to any user and allows the confirmation of the results described in this manuscript and the exploration of other results of interest. This easy-to-use resource is divided into different sections: [1] summary of analysis results in each phase. Then, for each of the datasets (BCC and BCDE), the detailed results of the [2] exploratory analysis, [3] differential expression, and [4] pathway analysis, as well as [5] the protein-protein interaction analysis results. The user can interact with the web tool through graphics and tables and search for specific information for a protein or pathway.