Prolonged treatment with Doxorubicin leads to various alterations in endothelial cells that are associated with a different response pattern to the therapy

We utilised the MMTV-PyMT (Mouse Mammary Tumour virus-polyoma middle tumour-antigen) [22,23,24] transgenic mouse model of breast cancer to examine the impact of tumour ECs on Dox sensitivity and development of chemoresistance. We implanted cancer cells isolated from MMTV-PyMT-tumours into wildtype (MMTV-PyMT-negative) FVB/N females to ensure consistent tumour growth (Fig. 1a). To evaluate the effectiveness of prolonged Dox treatment (4 mg/kg dose, one time per week, over three weeks), we compared daily measurements of tumour size (Fig. 1b; Fig. S1a). In contrast, to the placebo control group (Pla), the Dox-treated group exhibited two distinct tumour growth response patterns. Some tumours showed frequent increases and sharp decreases in tumour size (fold change) throughout the treatment schedule, which we classified as Responders (Rsp) to prolonged Dox treatment. However, others showed less variation in tumour size fold change with fewer growth reductions, and we categorised these as Partial Responders (Prt) to prolonged Dox treatment. Our analysis further revealed that the Rsp group had a higher percentage of days with size reduction compared to the Prt group (Fig. S1b). Additionally, while both groups displayed similar tumour growth rates at the beginning of treatment, the growth rate for the Prt group increased significantly in the later stages of treatment suggesting a partial response to Dox. In contrast, the Rsp group did not acquire this increase in growth rate at later stages of treatment suggesting a more robust response to Dox therapy (Fig. 1b; Fig. S1a-b). Similar tumour growth patterns were observed in additional independent experiments (Fig. S1c-d).

Fig. 1

Prolonged Doxorubicin treatment promotes distinct response patterns in the MMTV-PyMT breast cancer model. (a) Cartoon showing the prolonged Doxorubicin treatment model (created with BioRender (BS294BTN5B). Wildtype (FVB/N) females were implanted with 10^5 cancer cells isolated from treatment naïve MMTV-PyMT (FVB/N) tumour bearing mice. When tumours reached approximately 100 mm³ in size, they were treated weekly, intraperitoneally, with 4 mg/kg of Doxorubicin (Dox) for three weeks, while size was monitored throughout the whole experiment. (b) Plots depict tumour size related to the previous measure size (as fold change). Placebo (Pla, n = 3 mice), Responders (Rsp, n = 6 mice), and Partial Responders (Prt, n = 6 mice). (c) Principal Component Analysis graph showing tumour endothelial cells (ECs) and tumour cancer cells (Non-ECs) bulk-RNA sequencing samples; prepared from Pla, Rsp and Prt to Dox treatment mice (n = 20 cell preparations in total; n = 3 for EC Pla, n = 2 for Non-EC Pla, n = 4 for EC Rsp, n = 5 for Non-EC Rsp, n = 3 for EC Prt, and n = 3 for Non-EC Prt). Bar plot shows top 20 pathway terms significantly deregulated in (d) ECs and (e) in Non-ECs, Prt versus Rsp comparisons. Colour gradient of bars represents adjusted p-values (Fisher’s exact test with Benjamini-Hochberg post hoc test for multiple hypotheses). (f) Heatmap shows differential expression of various secreted molecules (secretome) from Non-ECs and ECs isolated from tumours from Prt and Rsp to Dox-treated mice, Prt versus Rsp comparisons (DE as fold change, from normalized counts transformed using voom)

To elucidate the factors influencing these differential responses to Dox, we examined transcriptomic changes in tumour ECs, which serve as barriers to drug delivery and are implicated in therapeutic resistance [7, 13, 19]. ECs were isolated from tumour cell suspensions, prepared from tumours that were harvested one week after last treatment, using magnetic antibody cell separation (MACS) combined with CD31 and ICAM2 antibodies specific for ECs (Fig. 1a-b). The unbound fraction of cells was enriched in epithelial cancer cells (Non-ECs), as reported previously [25].

Total mRNA from tumour ECs and Non-ECs preparations was subjected to bulk-RNA sequencing (Fig. 1a). Analysis revealed distinct enrichment of endothelial markers in tumour ECs, while Non-ECs exhibited increased epithelial markers (Fig. S1e). Principal Component Analysis (PCA) showed a clear differentiation between tumour EC and Non-EC samples, with Dox-treated tumour ECs clustering separately, suggesting distinct regulatory mechanisms (Fig. 1c). To gain insight into the distinct changes promoted by prolonged Dox in tumour ECs, we performed differential expression gene analysis by comparing Prt and Rsp, in both Non-ECs and tumour EC samples (Fig. S2a). Distinct regulatory patterns in tumour ECs were observed between the Prt and Rsp groups. The proportions of upregulated and downregulated unique genes differed compared with Non-EC (Fig. S2b). Compared with placebo, the EC Prt group showed less deregulated unique genes than the Rsp group (Fig. S2c), thus highlighting distinct regulatory responses between the two Dox-treatment groups.

EnrichR analysis [26,27,28] showed significant and distinct pathway changes in ECs and Non-ECs between Prt and Rsp to Dox groups (Fig. 1d, e). Tumour ECs had alterations in NF-kB transcription, secretome production, and cytokine signalling pathways (Fig. 1d). Compared with placebo, Rsp ECs mainly showed differences in replication and cell division pathways (Fig. S1f), whereas Prt ECs showed changes in NF-kB signalling, the secretome, and extracellular matrix pathways (Fig. S1g). Furthermore, differential expression of various secreted molecules was identified between Prt and Rsp tumour ECs (Fig. 1f). Collectively, our data suggests that Dox treatment reshapes EC transcriptional states in a manner that associates with divergent tumour responses to Dox therapy.

Identification of distinct endothelial cell populations in prolonged Doxorubicin treated breast tumours

To investigate the effects of prolonged Dox treatment on tumour EC dynamics at single-cell resolution, we carried out single-cell RNA sequencing (sc-RNAseq). For this, tumour EC enrichment was achieved by combining MACS and fluorescence-activated cell sorting (FACS) for CD31 + CD45- viable cells (Fig. 2a). Transcriptomic analysis excluded cells lacking CD31/Pecam1 and Cdh5 and verified the absence of markers for other cell types, including lymphatic (Lyve1), epithelial (Cdh1, Ktr5), fibroblast (Acta2), pericyte (Pdgfrb), and immune (Ptprc) cells (Fig. S3a-b). Unsupervised cluster analysis identified twelve distinct tumour EC populations (Fig. 2b-c), represented across all samples (Fig. 2d). The typical arterial-capillary-vein axis was assessed via specific markers for arterial (Gja4, Hey1), capillary (Crxc4, Cd36), and vein (Vwf, Ackr1), revealing populations with arterial characteristics in clusters 5 and 0, capillary markers in clusters 6 and 11, and venous markers in clusters 1 and 10, with cluster 10 classified as postcapillary veins (PCV) due to the additional expression of Selp and Vcam1 (Fig. 2e).

Next, we sought to identify characteristic markers that signify the angiogenic properties of tumour ECs as described previously [17, 29,30,31,32](Fig. 2f). Cxcr4, Trp53i11, and Pdgfb were notably enriched and detected in cluster 7 (Fig. 2f). Additionally, pathways related to migration, chemotaxis, and extracellular matrix remodelling were significantly represented in this cluster (Fig. S4a), further confirming the identity of these cells as “Tip” ECs [31, 32]. Interestingly, some markers associated with this population were also found in another cluster, which we referred as cluster 2, albeit at much lower levels. Literature review indicates that this cluster displays genes typical of a specific angiogenic population known as “Breach” ECs [32]. Furthermore, cluster 4 expresses markers indicative of a precursory state of these “Breach” ECs and we therefore classified it as “Pre-Breach” (P-Br) ECs (Fig. 2f). These cells have only been identified in tumours and are believed to assist “Tip” cells in spreading due to their podosome-like characteristics [32]. Consequently, pathways related to migration, chemotaxis, and extracellular matrix remodelling were also significantly represented in this cluster, similar to what was observed in cluster 7 (Fig. S4a).

We then examined proliferative cells by assessing Ki67 expression, alongside with cyclins and Cdks expression, revealing that cluster 9 comprises “Proliferative” ECs (Fig. 2f). Notably, Cdk4 was also expressed in cluster 3, indicating these cells may be in the cell cycle stage. This cluster is characterised by an enrichment of immature markers, suggesting it represents cells that are transitioning into or out of the cell cycle, thus we have designated them as “Immature” ECs (Fig. 2f). Pathway analysis further demonstrated that cluster 9 exhibits a significant upregulation of cell cycle and mitotic pathways (Fig. S4a).

Gene expression related to nucleotide metabolism, oxidative phosphorylation, fatty acid oxidation, and glycolysis across clusters (Fig. 2g) were identified, revealing that both the “Proliferative” (cluster 9) and “Immature” (cluster 3) cells are metabolically active, exhibiting high levels of glycolytic and oxidative phosphorylation pathways to fulfil their energy demands. “Tip” cells (cluster 7) are characterised by the expression of glycolytic enzymes that facilitate migration, corroborating previous findings [33]. Fatty acid oxidation is essential for the migration and survival of active, capillary proliferating EC (“stalk” cells) during angiogenesis. It also allows mature capillary EC (“phalanx” cells) to maintain a non-proliferative, quiescent state and support vascular integrity and adaptation to metabolic stress [34,35,36]. This metabolism was found to be restricted to capillary clusters 6 and 11. Interestingly, we discovered a subpopulation of ECs that has not been described in the literature. This cluster (Cluster 8, “Oxidative-active”) exhibited similarities with “Immature” and “Proliferative” cells, characterised by increased energetic metabolic pathways and a lack of specific angiogenic markers, while also expressing the highest levels of Dox clearance markers (Akr1a1, Cbr4). Furthermore, analysis of Hallmark pathways indicated possible involvement in hypoxic regions and of glycolysis pathways (Fig. S4b), reported features of an aggressive tumour cell phenotype [37].

Dox treatment and limited nutrient availability create stress in ECs that they must overcome for survival. Our analysis identified three clusters enriched in genes linked to EC stress response and activation: Hsp1a1, Plk2, Junb, Sox4, and Cyr61 (Fig. 2h). These include cluster 0, cluster 6, and cluster 10. Cluster 6 consists of active capillary ECs, while cluster 10 is characterised by postcapillary veins. Cluster 0, which shares markers with arterial cells (Stress-active-Art) (Fig. 2e), shows limited metabolic activity (Fig. 2g). Notably, genes associated with the NF-κB pathway, such as RELA and NFKBIA, are enriched in these clusters (Fig. 2h), with Hallmark pathways analyses revealing significant enrichment for the TNF/NF-κB pathway (Fig. S5). Overall, these findings highlight the presence of a highly heterogenous EC landscape in prolonged Dox treated tumours.

Fig. 2

Fig. 2. Distinct endothelial cell populations detected by sc-RNAseq in Doxorubicin treated MMTV-PyMT tumour bearing mice. (a) Experimental design for the single-cell RNA sequencing (sc-RNA-seq) approach. Created with BioRender (MI294B6DD9). (b) UMAP visualises single endothelial cells (EC) isolated from MMTV-PyMT tumours from Placebo and Doxorubicin (Dox)-treated mice, after filtering low quality and all non-vascular endothelial cells (n = 3 pooled tumour samples per condition). Numbers in parenthesis represent cell counts. (c) Unsupervised clustering of distinct single cell EC populations. Numbers in parenthesis represent cell counts. (d) Percentage of cells in each EC cluster. Bubble plots for (e) arterial, capillary and vein (f) angiogenic, (g) drugs and cell metabolism, and (h). stress response and NF-kB pathway markers expressed by each EC cluster. Bubble size represents the percentage of cells in the cluster expressing a particular marker, and bubble colour represents marker expression levels (generated with ShinyCell app[91] based on normalized expression values)

NF-kB pathway activation influences angiocrine signals in tumour endothelial cells from prolonged Doxorubicin treated mice

We conducted a comparative analysis of the abundance of cellular subpopulations in prolonged Rsp and Prt Dox-treated samples (Fig. 3a-b). Notably, angiogenic cell clusters accounted for 50–60% of total cells in both Dox-response samples, reflecting the high angiogenic characteristics of breast tumours (Fig. 3b-c). Multiple studies support that cancer therapy, including chemotherapy can promote angiogenesis, tumour regrowth, and progression after therapy [38,39,40,41,42]. We observed a decrease in Rsp ECs compared with partial Prt ECs, attributed to a higher proportion of proliferative cells (cluster 9) and oxidative-active cells (cluster 8) in the Prt group. Additionally, there was a significant reduction in the ratio of larger vessels versus small vessels, with a higher capillary (clusters 6 and 11) and postcapillary venule (PCV, cluster 10) proportions (Fig. 3b-c). Hallmark terms pathway analysis revealed that cell cycle regulation, interferon responses, Myc activation, hypoxia response, and NF-kB signalling were among the top pathways differentially regulated between Prt and Rsp across all tumour EC clusters (Fig. 3d). Given that our bulk-RNAseq showed a deregulation of the NF-kB pathway (from Fig. 1d) and differences in the expression of angiocrine signals between Prt and Rsp ECs (from Fig. 1f), we aimed to investigate the expression of angiocrine molecules to determine which clusters of ECs produce them.

Certain angiocrine molecules, such as Cxcl12, Il6st, Tgfb1, Tnfsf10, and Tnfsf12 were widely expressed, although their levels varied across different clusters (Fig. 3e). Notably, activated capillaries and postcapillary vein ECs (clusters 6 and 10) produced the highest levels of angiocrine molecules and were primary sources of various chemokines and cytokines involved in immunomodulation, including Il6, Il6st, Cxcl9, Cxcl10, Tnf, Csf1, Csf3, Cxcl12, Cxcl1, and Cx3cl1. They also presented the greatest number of significant pathway terms related to cytokines and interleukins in the gene and Hallmark terms analysis (Fig. 3e; Fig. S6; Fig. S7). Comparison between the two response groups to Dox, showed higher expression of angiocrine molecules in Prt ECs (Fig. 3f). Moreover, the increased expression of these genes in Prt ECs corresponds with pathways for stress and NF-κB (Fig. 3g), suggesting that resistance to prolonged Dox treatment may be linked to higher angiocrine expression molecules in clusters 6 and 10 due to NF-kB pathway activation.

To investigate the initial changes in angiocrine signals, we next analysed the acute Dox treatment (hereinafter referred to as Shrt), which reflects the first week of the prolonged treatment, where tumour growth rates were similar between the Rsp and Prt groups (Fig. 1b). All previously described EC populations were represented across samples (Fig. S8a-b) with about half identified as angiogenic (Fig. S8c). After Shrt Dox treatment, we observed a reduction in immature and proliferative ECs (clusters 3 and 9) and a slight increase in oxidative ECs (cluster 8) (Fig. S8c). Myc targets showed significant differences in the Shrt EC clusters, possibly accounting for the reduction in proliferative ECs (cluster 9) (Fig. S8d). In contrast to prolonged Dox treatment, the NF-kB pathway did not show significant changes, although clusters 6 and 10 show variation in some stress and NF-kB pathway gene expression (Fig. S8e). Characteristic angiocrine molecules in these clusters also exhibited varied regulation levels (Fig. S8f).

Fig. 3

Angiocrine factors alterations in specific endothelial cells populations associates with response to prolonged Doxorubicin treatment. a. UMAP visualises single endothelial cells (ECs) isolated from MMTV-PyMT tumours from Responders (Rsp), and Partial Responders (Prt) prolonged Doxorubicin (Dox)-treated mice. (n = 3 pooled tumour samples per condition). Numbers in parenthesis represent cell counts. b. Cell cluster abundance in percentage, in each Dox treatment group. c. Abundance of clusters (in percentage) in each indicated EC cluster category, in each Dox-treated group. d. Significantly deregulated Hallmark pathway terms in Prt versus Rsp ECs. Colour gradient scale represents adjusted p-values (Fisher exact’s p-value). e. Bubble plot shows representative angiocrine molecules expression within each cluster in Rsp versus Prt ECs. f. Angiocrine molecules, and g. Stress response and NF-kB pathway molecules from angiocrine expression within clusters 6 and 10, in Rsp versus Prt ECs. Bubble size represents the percentage of cells in the cluster expressing a particular marker, and bubble colour represents marker expression levels. Plots generated with ShinyCell app[91] from normalized gene expression data.

Doxorubicin-driven angiocrine signals influence the response of breast cancer cells to treatment in vitro

To study the impact of angiocrine signals from ECs in response to Dox on tumour cancer cells behaviour in vitro, we used Dox at 0.125µM dose. This dose was chosen based on that (1) isolated tumour cancer cells from treatment-naïve MMTV-PyMT tumours (Non-EC fractions) showed approximately a 50% reduction in survival rate when assessed by crystal violet assays (Fig. S9a), (2) previous work from our laboratory [7] shown that 0.125µM Dox treatment of endothelial cell cultures generates cytokine enriched conditioned media, thus a relevant feature for our current study; and (3) is similar to that used in the in the literature [7, 43, 44]. Isolated ECs from treatment-naïve MMTV-PyMT tumours (EC fractions), were treated with 0.125µM Dox for 24 h and then subjected to cytokine arrays analysis. Dox treatment of ECs resulted in the upregulation of several angiocrine molecules, including Ccl5, Ccl2, Cxcl10, Cxcl16, Cxcl1, and Selp (Fig. 4a, Fig. S9b-d). Notably, we found that NF-kB is a major regulator of these significantly regulated molecules (Fig. 4b) and that bulk-RNAseq results further support the increased expression of these molecules (Fig. S9e). Moreover, these angiocrine molecules were predominantly upregulated in clusters 6 and 10 in our sc-RNAseq analysis, in the Prt ECs when compared to Rsp ECs (Fig. 4c).

Tumour spheroid cultures provide a three-dimensional model that closely mimics the architecture and microenvironment of in vivo tumours, making them valuable for evaluating anticancer drug responses [45,46,47]. To analyse the effect of Dox-driven angiocrine signals on cancer cells, treatment-naïve MMTV-PyMT tumour-derived cancer cells, were used to form spheroids, and subsequently were treated with conditioned media from ECs pre-treated with Dox and subsequently washed (Fig. 4d-e). In comparison to conditioned media (CM) from PBS control (Veh)-treated ECs, CM from Dox-treated ECs protected cancer cells from DNA damage and death, indicated by a reduction in p-H2AX and active Caspase-3 levels respectively (Fig. 4d-e; Fig. S10a-b). Further, crystal violet survival assays performed in MMTV-PyMT tumour-derived cancer cell in two-dimensional cultures, subjected to same CM and Dox treatments above described, shown that CM from Dox-treated ECs protected Dox-treated-, but not Veh-treated-, cancer cells from death (Fig. S11a).

Fig. 4

Doxorubicin-driven endothelial angiocrine signals promote molecular changes in tumour cancer cells related with development of resistance to Dox-treatment. a. Cytokine arrays analysis of endothelial cells (ECs) isolated from MMTV-PyMT tumours and treated in culture with PBS vehicle (Veh) or 0.125µM Doxorubicin (Dox) for 24 h. Selected cytokines and quantification bar plot are shown. n = 4 replicates; NO = not statistically significant; YES = p < 0.05 (two-sided Student’s t-test). (b). Bar plot of Trrust transcriptional factors using the significantly deregulated cytokines from the arrays in (a). Colour gradient of bars represents adjusted p-values (Fisher exact’s p-value). (c) Bubble plot shows expression of significant deregulated genes from arrays in (a) assessed in the angiocrine clusters 6 and 10 from the sc-RNAseq analysis. Genes were cluster according to their expression across the samples. Colour gradient of bars represents adjusted p-values (Fisher exact’s p-value). Representative images from MMTV-PyMT tumour-derived cancer cell spheroids pre-treated with conditioned media from ECs treated in culture with PBS vehicle (Veh) or Doxorubicin (Dox) for 24 h, cancer cells were next treated with Veh or 0.125µM Dox for further 48 h, and last stained for (d) p-H2AX or (e) Cleaved Caspase-3 (CC-3). Scales, 100 μm. Quantification plots of p-H2AX area or CC-3 area per DAPI-positive area are shown. ****p < 0.0001 (two-sided Student’s t-test). (f). LS/MS phosphoproteomics experimental design for MMTV-PyMT tumours-derived cancer cells sample analysis. Created with BioRender (CW294BUSKT). (g). Kinase Significant Enrichment Assay (KSEA) analysis. FC, fold change. n = 4 cell lysates for each condition; *P < 0.05, **P < 0.01, ***P < 0.001 (two-sided Student’s t-test). (h) Kaplan-Meier plot shows disease progression/survival in high and low phospho-PRKACA expression human breast cancer patient groups from LinkedOmics database [48]. Cox proportional hazard ratios, with the rank P value is shown

Mass Spectrometry phosphoproteomic analysis was performed to decipher changes occurring in the tumour cells after the treatment with the EC CM, which may protect them from Dox damage, and promote the resistance to the treatment. Analysis showed that CM from EC treated with Dox promoted significant reduction of ERK2 and mTOR and, in contrast, PRKACA was found significantly upregulated (Fig. 4f-g). Taken together these changes in signalling pathways may lead to the functional changes observed in the tumour cancer cells as we observed significant regulated pathways related with DNA damage repair from the phosphoproteomic data (Fig. S11b). Furthermore, to test whether our results could extend to human disease, we interrogated human breast cancer patient phosphoproteomic datasets [49], we found that high phospho-PRKACA expression levels in breast cancer patients associates with poor survival outcome (Fig. 4h). Overall findings not only delineate the heterogeneity of EC responses but also demonstrate how alterations in NF-κB/PRKACA dependent angiocrine signalling contribute to the development of partial resistance to Dox.