Tumor-derived exosomes RNA expression profiling identifies the prognosis, immune characteristics, and treatment in HR+/HER2-breast cancer

Tissue collection and public data sets

From February to March 2022, a cohort of consecutive patients with early-stage breast cancer was enrolled at the Tangshan People’s Hospital. These patients were required to meet all the following inclusion criteria:1) age between 55 and 65 years with no prior history of breast cancer or other malignancy and without any past chemotherapy treatment, 2) confirmed diagnosis of breast cancer through needle biopsy histology, and 3) no history of medical intervention for breast cancer prior to surgical resection. The following exclusion criteria were applied: inflammatory breast cancer, hypertension (blood pressure > 140/90 mmHg), hyperlipidemia (triglyceride > 1.7 mmol/L, total cholesterol > 6.0 mmol/L), diabetes, and tumors with a diameter larger than 5 cm.

Carcinoma and healthy tissue samples were obtained from surgical specimens immediately after surgery. Fresh tissues were used to prepare single-cell suspensions and subsequent primary cultures, which were then frozen and stored at -80° C for later analysis.

In The Cancer Genome Atlas (TCGA) Breast Cancer (BRCA) project, masked copy number segments, RNASeq expression (STAR – Counts), and clinical data were downloaded using the TGCAbiolinks R package [8]. The TCGA RNA sequencing FPKM data extracted from “STAR – Counts” files were transformed into log2 (FPKM + 1). Survival data and breast cancer subtypes were compiled from clinical data and a total of 532 patients with HR+/HER2-breast cancer with a follow-up time of >30 days were selected (Supplementary Table 4).

The Gene Expression Omnibus (GEO) dataset GSE25066 [9], which includes 292 patients with HR+/HER2- breast cancer with distant relapse-free survival (DRFS) data (Supplementary Table 5), was used for validation.

Dataset GSE145325 [10] contains RNA-seq data of patients with ER+ breast cancer treated with letrozole, and SRR files (SRP249306) were downloaded from Sequence Read Archive (SRA) stores and converted to FASTQ format using the SRAtoolkit.

Cell lines and cell culture

The human breast epithelial cell line MCF10A was procured from the Peking Union Medical College Cell Resource Center (PUMCCRC) in Beijing, China. The MCF-7 human breast cancer cell line was obtained from the Shanghai Cell Bank of the Chinese Academy of Sciences (CAS). Cells were cultured according to the cell instructions respectively.

Exosomes isolation

Exosomes present in the supernatants of breast cancer cell cultures were purified using differential ultracentrifugation. Briefly, the cells were grown in cell culture medium supplemented with 10% Fetal Bovine Serum (FBS)-exosomes depleted (SBI, USA) for 48 h. Thereafter, the cell culture supernatant was collected and subjected to centrifugation at 2× 10³ g for 20 min at 4° C to eliminate cells. A subsequent centrifugation at 1× 10⁴ g for 30 min at 4° C was performed to remove cellular debris. The resulting supernatant was sieved using 0.2 μm filters (Millipore, USA) and subsequently ultra-centrifuged at 1× 10⁵ g for 1h at 4° C. The pellets were resuspended in phosphate-buffered saline (PBS) and once again subjected to ultra-centrifugation at 1× 10⁵ g for 1h at 4° C.

Purified exosomes transmission electron microscopy

A volume of 10 μl of exosomes was administered onto copper transmission electron microscopy grids, measuring 3.05 mm with 200 mesh, and left to rest for 5 min. The grids were rinsed with PBS and coated with 2% uranyl acetate for 3 min. Exosome images were captured using a transmission electron microscope (Tecnai G2 Spirit Biotwin; FEI Company, USA).

RNA isolation and qRT-PCR

Cells total RNA was extracted using Trizol Reagent (Invitrogen) according to the manufacturer’ protocol. Subsequently, 1 μg of RNA was subjected to reverse transcription using the PrimeScript™ RT Master Mix kit (Takara, China). The SYBR premix Ex TaqTM II kit (Takara, China) was used to detect the expression levels of the specified genes, and the results were analyzed using the Stratagene Mx 3000P software (Agilent Technologies, USA). The 2^–ΔΔCt method was used to calculate the relative expression levels of mRNA. Primer sequences used in this analysis are listed in Supplementary Table 1. All qRT-PCR experiments were repeated in triplicate, and each group had three technical replicates. Statistical analyses were performed using GraphPad Prism 8.0. The outcomes of the experiments were expressed as the mean ± SD, and a two-tailed Student’s t-test was employed to compute the p-value. Statistical significance was considered when p-values were <0.05.

Microarray analysis and data processing

RNA quantity and quality were measured using a NanoDrop ND-1000. RNA integrity was assessed via standard denaturing agarose gel electrophoresis or an Agilent 2100 Bioanalyzer.

Sample labeling and array hybridization were performed according to the Agilent One-Color Microarray-based Gene Expression Analysis protocol (Agilent Technologies). The Agilent Feature Extraction software (version 11.0.1.1) was used to analyze the acquired array images. Long non-coding (lnc) RNAs and mRNAs that were flagged at least in three out of six samples as Present or Marginal (“All Targets Value”) were chosen for further data analysis.

Microarray probes re-annotation

The Human reference sequence (GRCh38.d1.vd1) and annotation file (GDC.h38 GENCODE v36 GTF) were downloaded from The Genomic Data Commons (GDC) database (https://gdc.cancer.gov/about-data/gdc-data-processing/gdc-reference-files).

We mapped the probe sequences to the human genome (GRCh38.d1.vd1.fa.tar.gz) using the SeqMap software (http://www-personal.umich.edu/~jianghui/seqmap/). Using an annotation file (gencode.v22.annotation.gtf.gz), we reannotated the probes for the two chips and removed those corresponding to multiple genes.

Identification of prognostic genes and construction of a risk model

Univariate Cox regression analysis was performed to explore the correlation between genes and OS in TCGA dataset. DEGs with p <0.05 were considered candidate genes.

Random forest analysis was used to reduce the scope of the gene screening. Finally, the coefficient for each gene was obtained through multivariate Cox regression, and the exosome score (ES) = sum of coefficient × expression level of the gene was calculated.

Immune-related features and tumor microenvironment data

Thorsson conducted an extensive immunogenic study of more than 10,000 tumors containing 33 different types of cancer using data collected from TCGA [11]. We used Thorsson’s outcomes to investigate the association between ES and immune-related features. Similarly, Tamborero et al. [12] provided data on immune infiltration patterns, which we used to investigate the association between the ES and tumor immune microenvironment. All data were numerical and transformed into standardized z-scores for statistical analysis and presentation.

Data availability

Microarray data is available in the Gene Expression Omnibus: GSE207304 Source code for model construction is available at https://github.com/QZhou-Ch/BC-Exosomes.

Statistical analysis

Most statistical analyses were performed using R (version 4.2.0) with default arguments unless mentioned otherwise. The Shapiro–Wilk test was used to test data normality. The Wilcox test was performed to verify the statistical significance between two groups, whereas the Kruskal–Wallis test was applied to test for multiple groups.

Differential expression analysis was performed using the R package DESeq2 [13]. Survival and survminer packages were used for survival analysis, and the time-dependent receiver operating characteristic (ROC) and area under the curve (AUC) were determined using the R package survival ROC. P-values were two-sided and adjusted for multiple testing using the Benjamini–Hochberg False Discovery Rate (FDR), and statistical significance was set at p < 0.05 or FDR < 0.1.

Differential gene expression in exosomes of luminal breast cancer cells and healthy breast cells

Three pairs of breast cancer and healthy breast tissues were collected at the Tangshan People’s Hospital to establish primary cultures (Supplementary Table 2). Subsequently, the supernatant of the cell culture medium was collected, and exosomes were extracted via differential ultracentrifugation (Figure 2A). Total RNA was extracted from exosomes and used to determine gene expression patterns using microarrays. Background correction was performed with the backgroundCorrect function using the “normexp” method, and normalization was performed with the normalizeBetweenArrays function using the “normexp” method. When a gene contained multiple probes, the maximum value was considered as the gene expression value. Differential expression analysis was performed using the R Bioconductor package, limma. DEGs were selected using a fold-change cut-off of >1.5 or < -1.5 and a p-value < 0.05. In total, 817 DEGs were identified in cancer vs. healthy exosomes, including 315 upregulated tumor-derived exosomes genes (UTEGs) and 502 downregulated tumor-derived exosomes genes (Figure 2B, 2C).

Figure 2. Differentially expressed genes from exosomes. (A) Exosomes were observed using transmission electron microscopy. Scale bar, 200 nm. (B) Heatmap of Different Gene Expression in exosomes. (C) Volcano plot of Different Gene Expression in exosomes. (D) Consensus clustering matrix for k = 2. (E) Consensus clustering CDF for k = 2–10. (F) Kaplan–Meier survival curves of OS according to two clusters in TCGA cohort. CDF, cumulative distribution function; OS, overall survival; TCGA, The Cancer Genome Atlas.

UTEGs cluster analysis

We used the 315 UTEGs in patients with HR+/HER2- BRCA from TCGA to perform a consensus cluster analysis (maxK = 10, reps = 1000, clusterAlg = “km,” distance=“euclidean”) and the proportion of ambiguous clustering (PAC) score showed an optimal K of two (Figure 2D, 2E). The 532 patients with breast cancer were divided into two categories according to the optimal K value: ExoA and ExoB, with 244 and 288 cases, respectively. Kaplan–Meier survival curves (Figure 2F) revealed that patients in the ExoA group had poorer survival than those in the ExoB group (HR, 0.286; 95% CI, 0.136-0.600; log-rank p = 0.33 e-03). Therefore, we believe that this new type of cluster classification will be beneficial for accurate clinical diagnosis and treatment.

Development and evaluation of ES in TCGA cohort

To meet the cross-platform availability of the model, we hope that it can be applied to the Affymetrix Human Genome U133A Array platform because it is the most widely used chip with the largest amount of data available. After Univariate Cox analysis of the TCGA dataset, five UTAGs genes were eventually included in the candidate model. Feature importance was determined via random forest survival analysis (Figure 3A, 3B) with default parameters (ntree = 5000), and we used four genes to construct the final model: PDPK1, WSB2, PIR, and MTHFD2 (details in Supplementary Table 3).

Figure 3. Construction and validation of exosome scores. (A) The relationship between error rate and the number of trees. (B) Random Forest feature importance ranking for the five predictive features.
(C) Kaplan–Meier survival curves of OS according to ES groups in TCGA cohort. (D) ROC curves at 3-, 5-, and 7 years of OS according to ES groups in the TCGA cohort. (E) Forest plots of OS. (F) Kaplan–Meier survival curves of DRFS according to ES groups in the GEO validation set. (G) ROC curves at 3- and 5 years of DRFS according to ES groups in the GEO validation set. ES, exosome score; OS, overall survival; TCGA, The Cancer Genome Atlas; ROC, receiver operating characteristic; AUC, area under the curve; DRFS, distant relapse-free survival; AIC, Akaike information criterion; GEO, Gene Expression Omnibus.

The patients in the TCGA set were split into two groups, high- and low-score groups, according to the median ES. The survival curves (Kaplan–Meier estimates) revealed that the OS in the low-score group was significantly higher than that in the high-score group (HR, 0.170; 95% CI, 0.081–0.359; log-rank p = 0.380e-04) (Figure 3C). The AUC of ES for OS were 0.754 at 3 years, 0.739 at 5 years, and 0.783 at 7 years (Figure 3D).

Moreover, to evaluate the independent prognostic effect of ES, we used a multivariate COX regression model to adjust for other factors, including age and TNM stage. Results indicate that ES remained a significant and independent prognostic indicator in TCGA cohort (Figure 3E).

Prognostic power evaluation of ES in the GEO cohort

To validate the robustness of the risk signature, the ES for each patient with HR+/HER2- BRCA in the GEO cohort was also evaluated. The median ES was used to divide patients into high- and low-score groups. Kaplan–Meier survival curves revealed that DRFS was significantly prolonged in the low-scoring group (HR, 0.4594778; 95% CI, 0.241–0.877; log-rank p = 0.020) (Figure 3F). The AUC of IS were 0.634 and 0.683 at three and five years, respectively (Figure 3G).

Different immune-related features between high- and low score patients in TCGA cohort

First, we examined the distribution of ES in five immune types in patients with breast cancer. A substantial difference was observed between the five immune subtypes of the Kruskal–Wallis test (p = 0.016), as shown in Figure 4A. Furthermore, patients with type C3 (Lymphocyte Depleted) had a lower ES than type C2 (IFN-γ Dominant) and type C4 (Lymphocyte Depleted) subtypes.

Figure 4. Exploration of the exosome score predictive role in the TCGA cohort. (A) Raincloud Plot shows the comparison of ES between the different immune subtypes. (B) Lollipop plot showing the comparison of immune-related features between the low- and high-score groups. The length of the stick represents the difference between the medians of the features in the high- and low-score groups. (C) Circos Plot displays the CNVs distribution in ES-relevant groups, and red dots represent the difference between the incidence of CNVs in the high- and low-score groups. (D) Violin plot showing the distribution of the two groups in different immune cells, NA p > 0.05, * p < 0.05, ** p < 0.01, *** p < 0.001. ES, exosome score; TCGA, The Cancer Genome Atlas; CNVs, copy number variations; aDC, activated dendritic cells; CD8+ T, CD8 T cells; iDC, immature dendritic cells; NKbright, NK CD56 bright cells; NKdim, NK CD56 dim cells; Tcm, central memory T cells; Tem, effector memory T cells; Tfh, follicular helper T cells; Tgd, gamma delta T cells; Th, helper T cells; Treg, regulatory T cells.

In the following analysis, we explored the differences in the composition of the tumor immune infiltrate, somatic variation, immunogenicity, and genomic state between the ES-based subtypes of the TCGA dataset (Figure 4B).

The high-score group had higher proliferation, altered fraction, silent and non-silent mutation rates, as well as a reduced stromal fraction than the low-score group.

Potential factors that presented tumor somatic or germline mutations, including homologous recombination defects, intratumor heterogeneity, number of segments, aneuploidy score, non-silent silent mutation rates, were compared between the high- and low-score groups. The median values for the variables of the high-score group were substantially higher than those of the low-score group.

Finally, given the differences in tumor immune infiltrates between the low- and high-score groups, we analyzed the immune microenvironments in both groups. The Wilcox test revealed that high-scoring patients had significantly higher abundances of central memory T cells (Tcm) and T helper cells (Th), but lower abundances of B cells, immature dendritic cells (iDC), mast cells, neutrophils, NK CD56 bright cells (NKbright), effector memory T cells (Tem), and gamma delta T cells (Tgd) (Figure 4D) than the low-score group.

CNV incidence between low- and high-score patients in TCGA cohort

We found significant differences between the two groups in tumor somatic or germline mutations, which were both calculated from CNV files; therefore, we performed a statistical analysis of CNVs events in both groups of patients.

The ‘Masked Copy Number Segment File’ and ‘SNP6 GRCH38 Remapped Probeset File for CNV Analysis,’ were used for GISTIC2 analysis [14] on the GenePatter website [15]. The CNV values were further thresholder using a noise cutoff of 0.3. A chi-square test was performed to confirm the difference in focal CNVs events between the high- and low-score groups.

Patients in the high-score group had a higher incidence of CNVs events in 6861 genes than those in the low-score group (FDR < 0.1) (Figure 4C) (Supplementary Table 6).

ES can reflect the treatment sensitivity

We further identified associations between the model and the treatment effects of chemotherapeutic agents and endocrine therapy. The Genomics of Drug Sensitivity in Cancer (GDSC) database was adopted to evaluate the correlation between the ES and drug sensitivity values (IC50) of molecules in 50 breast cell lines (Tissue = “breast”). Consequently, vincristine, paclitaxel, and docetaxel, which target microtubules, appeared to be significantly and positively correlated with ES (Figure 5A–5C). Other common chemotherapeutic drugs, such as doxorubicin and cisplatin; and endocrine therapy drugs, such as tamoxifen and fulvestrant, were not significantly associated with ES (p>0.05, Figure 5D–5G). Tamoxifen, fulvestrant or anastrozole, and exemestane, mechanism lies in reducing the interaction of estrogen with cancer cells. Therefore, we evaluated the relationship between ES and endocrine therapy in another clinical study of letrozole treatment.

Figure 5. The ES predicting the sensitivity to chemosensitivity and endocrine therapy. Scatter plots for associations between drug sensitivity values (IC50) and (A) vincristine, (B) paclitaxel.
(C) docetaxel, (D) doxorubicin, (E) cisplatin, (F) tamoxifen, and (G) fulvestrant; (H) Raincloud Plot showing the comparison of ES between endocrine therapy non-responder and responder groups. ES, exosome score; IC50, half maximal inhibitory concentration.

Because the study did not provide the exact reference genome or annotation file, we re-aligned the reads using the fastp-Hisat2-featurecounts pipeline. Based on the counts data (Supplementary Table 7), we calculated FPKM and ES in GSE145325 dataset. The Wilcoxon test revealed that ES in the non-responder group was significantly higher than that in the responder group (Figure 5H).

Validation of the expression of four tumor-derived exosome genes

The mRNA expression levels of PDPK1, WSB2, PIR, and MTHFD2 were evaluated in exosomes extracted from the MCF-7 and MCF10A cell lines. The results suggest that the mRNA expression of the four genes was higher in MCF-7 exosomes than in MCF10A exosomes (Figure 6A–6D).

Figure 6. The expression of risk model genes. (A) MTHFD2 mRNA expression levels in exosomes from MCF10A and MCF-7 cell line; (B) PDPK1 mRNA expression levels in exosomes from MCF10A and MCF-7 cell line; (C) PIR mRNA expression levels in exosomes from MCF10A and MCF-7 cell line; (D) WSB2 mRNA expression levels in exosomes from MCF10A and MCF-7 cell line, * p < 0.05, ** p < 0.01, *** p < 0.001. MTHFD2, methylenetetrahydrofolate dehydrogenase 2; PDPK1, phosphoinositide-dependent protein kinase 1; PIR, pirin; WSB2, WD repeat, and SOCS box containing 2.