Single-cell sequencing analysis reveals the relationship between tumor microenvironment cells and oxidative stress in breast cancer bone metastases

Data collection and processing

BM-related single cell RNA sequencing (scRNA-seq) data were obtained from Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/), dataset GSE190772 based on the GPL24676 platform, including BM of two bilateral bone metastases collected in a patient initially diagnosed with ER+ primary breast cancer. In addition, dataset GSE131007 based on the GPL20301 and GPL24676 platforms, including three primary tumor tissues from mouse xenografts of human BC patient origin and one BM tissue. Among them, cells from murine cells or human-murine doublets were excluded. Total 3 BM tissues and 3 primary samples, and the primary samples were the control samples in this study.

Construction of single cell atlas

Single-cell data were merged using the IntegrateData function [13] of the Seurat package [14] in R language, and cell clustering analysis was performed according to default parameters, filtering for cells with top and bottom 1% gene count and >10% mitochondrial content. The clustering results were downscaled and visualized [15] based on a uniform manifold approximation and projection (UMAP) for dimension reduction technique and projected onto a two-dimensional image defined as a single-cell atlas. In addition, cell types were annotated according to cell markers known from previous studies [16].

Differential gene expression analysis

The differentially expressed genes (DEGs) in each cluster between single cells of primary tumor tissue and BM tissues were identified using the “FindAllMarkers” function, and differences with a adjusted P-values < 0.05 and |log fold change (logFC)| > 0.5 were considered significant.

Functional enrichment and gene enrichment analysis

To further explore the biological processes and pathways involved with genes that showed dysregulated expression in different cell clusters, Gene Ontology (GO) terminology and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were performed based on the expression of marker genes. The R package clusterProfiler [17] for enrichment analysis regarding biological processes (BPs) of GO and KEGG signaling pathways, and P<0.05 were considered significant.

Cellular scoring of oxidative stress-related gene sets

The AddModuleScore function [18] in the Seurat package was used to score oxidative stress-related pathways.

Pseudo-time analysis

The differentiation developmental trajectory of dysregulated cells in primary tumor tissues and BM tissues was reconstructed using the Monocle 3 package [19] in R language and visualized by UMAP. Subsequently, the cells were sorted according to their progression through the developmental program.

Gene regulatory network (GRN) analysis

In addition, using the Python module tool pySCENIC [20], this study comprehensively reconstructed the transcription factor-centered gene regulatory network to further explore the regulatory mechanisms of dysregulated cells.

The workflow started with describing the input single-cell expression level profile matrix, and then using a regression method for each target (GRNBoost2) to infer co-expression modules. The results allowed us to determine which indirect targets were trimmed based on the discovery of cis-regulatory patterns (cisTarget). Subsequently, AUcell was used to quantify the activity of those regulators by enriching and scoring the regulator target genes to obtain a regulon activity score (RAS). The single-cell data were further downscaled using the RAS matrix and a regulon specificity score (RSS) was calculated based on the Jensen-Shannon divergence (JS scatter) and used to identify regulators specific for certain cell populations. The most specific and significant regulons were mapped to single cell cluster profiles and validated using massively parallel sample sequencing (SEEK database). Finally, a connection specificity index (CSI) matrix was calculated, and the regulators were hierarchically clustered according to CSI to define regulator modules that could be used to identify relationships between regulator modules and regulators. Those relationships were then visualized using the R package ComplexHeatmap.

Cellular communication

Signal transduction emphasizes the manner and outcome of signal reception and the signal conversion after reception, with ligand-receptor binding being one of the main forms of signal transduction between neighboring cells. In this study, high confidence ligand-receptor interactions between subpopulations of cells were identified by the R language package iTALK. It preferentially identifies genes that are highly or differentially expressed in cell clusters that will be matched by a ligand-receptor database to discover important intercellular communication events.

Data analysis and statistics

Comparisons between the two groups were made using Student’s t test and correlation coefficients were calculated using Spearman analysis. P<0.05 was considered significant. Regarding the code used in this study we have uploaded it as Supplementary File 1.

Data availability statement

Data used in this study were obtained from Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo).

Global single-cell landscape of bone metastases from breast cancer

We try to draw the global single cell atlas of BM through scRNA-seq technology, and further explore the potential ecological panorama of BM, in order to find the potential therapeutic target of BM. The flow of this study was shown in Figure 1A. After standardized data processing and quality control, a total of 16,409 high-quality single-cell transcriptional profiles were captured and clustered to generate 43 cell clusters, and differential gene expression analysis revealed a wide range of gene expression dysregulation in different kinds of cell types in BM compared to controls. The cell clusters were further identified into eight cell types (Figure 1B), including BC cells, osteoblasts (OC), mesenchymal stem cells (MSCs), B cells, CD8⁺ T cells, CD4⁺ T cells, fibroblasts and endothelial cells (En). The markers positively expressed by the cells are consistent with recent published gene signatures such as scRNA-seq and laboratory studies (Figure 1C). In addition, BC cells were reduced in the BM microenvironment compared to controls, whereas MSCs, Fibroblasts and other cells were significantly more abundant in the BM (Figure 1D). In addition, we assessed the differences in single-cell oxidative stress levels between control and BM samples and found that GO_CELL_DEATH_IN_RESPONSE_TO_OXIDATIVE_STRESS and GO_RESPONSE_TO_OXIDATIVE_STRESS had higher scores in bone metastasis samples (Figure 1E). In summary, we initially constructed a global landscape of the dynamic single-cell ecology of the BM microenvironment by single-cell histology, and we found significant concomitant dysregulation of gene expression between different cell types and explored the altered cellular ecology of BM patients, in which BC, MSCs and fibroblasts may play an important role in BM. In addition to this, oxidative stress may play a key role in BM.

Figure 1. Global single-cell landscape of patients with bone metastases from breast cancer. (A) Flow chart underlying this study. Primary tumors were used as controls in the follow-up study. (B) Single-cell atlas mapping cell types. (C) Cell marker genes for annotation. (D) Differences in cell abundance between control and breast cancer bone metastasis patients. (E) Comparison of differences in single-cell oxidative stress levels between control and bone metastasis samples.

Landscape of BC cell subpopulations in breast cancer bone metastases

When cancer cells metastasize to the bone, these cells enter some part of the bone through the blood or lymphatic system and become metastatic cancer cells. When these cancer cells deposit in bone, they release substances that form osteoclasts and osteoblasts. This may account for the decrease in BC cells in Figure 1D. Based on the cellular ecological atlas at single-cell resolution, we explored the subpopulations of BC cells in depth and identified a total of 10 subpopulations of BC cells (Figure 2A) and found that these subpopulations were heterogeneous among different subgroups (Figure 2B). These BC cell subpopulations all significantly expressed their marker genes, with BC_MUC1, BC_SCGB2A2, BC_FN1, BC_BGN, and BC_PEG10 subpopulations being significantly more abundant in BM (Figure 2C, 2D). By enrichment analysis, we found that these subpopulations were significantly enriched in oxidative stress response and BPs associated with oxidative stress and cell death (Figure 2E). In addition, Oxidative phosphorylation, HIF-1 signaling pathway, TGF-beta signaling pathway, Wnt signaling pathway, Apoptosis-associated oxidative stress and cell death-related pathways, MAPK signaling pathway, p53 signaling pathway, PI3K-Akt signaling pathway Breast cancer, and JAK-STAT signaling pathway of cancer-related signaling pathways were also significantly enriched (Figure 2F). Meanwhile, in control and BM samples, BC cells were significantly different from GO_CELL_DEATH_IN_RESPONSE_TO_OXIDATIVE_STRESS, GO_CELLULAR_OXIDANT_DETOXIFICATION and GO_RESPONSE_TO_OXIDATIVE_ STRESS oxidative stress pathway scores were significantly different (Figure 2G and Supplementary Figure 1). Oxidative stress-related pathway scores were significantly higher in BM samples than in controls. These results further suggest that BC cells in BM are susceptible to oxidative damage and exhibit high levels of oxidative stress, which plays a key role in apoptosis.

Figure 2. Breast cancer cell subpopulations in patients with bone metastases from breast cancer. (A) Single-cell atlas showing breast cancer cell subpopulations. (B) Single-cell atlas showing breast cancer cell subpopulations in control and breast cancer bone metastasis patients. (C) Differential abundance of breast cancer cell subpopulations in control and breast cancer bone metastasis patients. (D) Marker genes specifically and highly expressed in subpopulations of breast cancer cells. (E, F) Biological processes (E) and signaling pathways (F) enriched in breast cancer cell subpopulations. (G) Comparison of the differences in oxidative stress levels in breast cancer cells between control and bone metastasis samples.

Clonal evolution of BC cells in bone metastases from breast cancer

We then depicted the differentiated developmental trajectories of tumor cell subpopulations. BC_MUC1, BC_SCGB2A2, and BC_BGN subpopulations were at the end of the developmental trajectory of tumor cells, while BC_SERPIINA6 subpopulation was at the beginning of the developmental trajectory and gradually differentiated into different other subpopulations as BM progressed (Figure 3A, 3B). Markers were clustered into four modules, where BC_MUC1 and BC_SERPIINA6 were regulated by the same transcription factor HOXB2 and BC_BGN was regulated by the TCF4 transcription factor (Figure 3C), and the expression of these transcription factors was mapped in a single cell atlas (Figure 3D). In conclusion, we further clarified the differentiation and development trajectory of tumor cell subsets in BM and the transcriptional regulation targets of different subsets.

Figure 3. Clonal evolution of breast cancer cell subpopulations in patients with bone metastases from breast cancer. (A, B) Pseudo-time values (A) and developmental trajectories (B) of breast cancer cell subpopulations, with pie charts representing the proportion of control and breast cancer bone metastasis patients in breast cancer cell subpopulations. (C) Co-expression modules of transcription factors in breast cancer cell subpopulations of patients with breast cancer bone metastases. Left: Identification of regulator modules based on the regulator’s linkage specificity index matrix. Middle: representative transcription factors and their binding patterns in the modules. Right panel: cell subpopulations in which transcription factors are located. (D) Single-cell atlas showing transcription factors regulating breast cancer cell subpopulations.

Landscape of fibroblast subpopulations in bone metastases from breast cancer

The aforementioned of others’ analysis revealed the presence of a large number of fibroblasts in bone metastasis sites, and fibroblasts are becoming important cellular players in bone metastasis [21]. Eight fibroblast subpopulations were further identified by subpopulation analysis (Figure 4A), and these fibroblast subpopulations were essentially highly enriched in BM (Figure 4B), with differentially expressed genes for each subpopulation displayed in Supplementary Table 1. These fibroblast subpopulations all expressed different specific markers, and studies of their variable levels revealed significant abundance of Fibroblasts_CLDN10 and Fibroblasts_S100P in controls, and Fibroblasts_IBSP, Fibroblasts_TAGLN, Fibroblasts_ASPN, and Fibroblasts_OLFML2B, Fibroblasts_KRT19 and Fibroblasts_MGST1 were significantly abundant in BM (Figure 4C, 4D). Further, the BPs and KEGG signaling pathways involved in these fibroblast subpopulations were explored, and fibroblast subpopulations were found to be significantly involved in the biological processes of ossification and bone remodeling (Figure 4E), in addition to ECM-receptor interaction, Cytokine-cytokine receptor interaction, Apoptosis, Focal adhesion, TNF signaling pathway, and TGF-beta signaling pathway were significantly enriched (Figure 4F). Taken together, our results suggest that fibroblasts in BM may regulate the ecological niche formation of bone metastases and play an important role for tumor cell inoculation into the bone marrow and growth.

Figure 4. Fibroblast subpopulations in patients with bone metastases from breast cancer. (A) Single-cell atlas showing fibroblast subpopulations. (B) Single-cell atlas showing fibroblast subpopulations in control and breast cancer bone metastasis patients. (C) Differences in abundance of fibroblast subpopulations in control and breast cancer bone metastasis patients. (D) Marker genes specifically and highly expressed in subpopulations of fibroblast subpopulations. (E, F). Biological processes (E) and signaling pathways (F) that enrich fibroblast subpopulations.

Clonal evolution of fibroblasts in bone metastases from breast cancer

The pseudo-time differentiation trajectory showed that the Fibroblasts_S100P and Fibroblasts_CLDN10 subpopulations were in an early position of differentiation development, while Fibroblasts_IBSP, Fibroblasts_TAGLN, Fibroblasts_ASPN, Fibroblasts_ OLFML2B, Fibroblasts_KRT19, and Fibroblasts_MGST1 were at the end stage and all were highly enriched in BM (Figure 5A, 5B). Subsequent GRN analysis showed that fibroblast subpopulation genes were organized into three modules (Figure 5C) and that different fibroblast subpopulations were regulated by different TFs to guide cell fate selection, respectively (Figure 5D). These results reflect the developmental trajectory of fibroblast subpopulations in BM and the transcriptional regulatory targets of different subpopulations.

Figure 5. Clonal evolution of fibroblast subpopulations in patients with bone metastases from breast cancer. (A, B) Pseudo-time values (A) and developmental trajectories (B) of fibroblast subpopulations, with pie charts representing the proportion of fibroblast subpopulations in control and breast cancer bone metastasis patients. (C) Co-expression modules of transcription factors in fibroblast subpopulations of patients with breast cancer bone metastases. Left: Identification of regulator modules based on the regulator’s linkage specificity index matrix. Middle: representative transcription factors and their binding patterns in the modules. Right panel: cell subpopulations in which transcription factors are located. (D) Single-cell atlas showing transcription factors regulating fibroblast subpopulations.

The MSC subpopulations landscape of breast cancer bone metastases

MSCs can differentiate into a variety of cell types, including osteogenic osteoblasts, chondrocytes, and adipocytes [22]. In addition, MSCs are a key component of tumor cell homing and adhesion to bone metastasis ecological niches [2]. Exploring the MSC subpopulation landscape by single cell resolution identified ten MSC subpopulations (Figure 6A), all of which were essentially significantly abundant in BM (Figure 6B), with differentially expressed genes for each subpopulation displayed in Supplementary Table 2. Further exploring the expression of specific markers for these subpopulations and differences in the abundance of subpopulations, the MSC_MARCKSL1 subpopulation was significantly more abundant in the control group, while all other subpopulations were significantly present in the BM group (Figure 6C, 6D). Enrichment analysis revealed that MSC subpopulations were involved in the regulation of BMP signaling pathway, tumor necrosis factor production, tumor necrosis factor superfamily cytokine production, and some skeletal growth and development and oxidative stress-related BPs (Figure 6E), Focal adhesion, Cell adhesion molecules, ECM-receptor interaction, which are KEGG signaling pathways associated with tumor metastasis, were enriched, in addition, PI3K-Akt signaling pathway, TGF-beta signaling pathway, Hippo signaling pathway and NF-kappa B signaling pathway, which are cancer-related KEGG signaling pathways, were also significantly enriched (Figure 6F). GO_RESPONSE_TO_OXIDATIVE_ STRESS scores were higher in the BM samples (Figure 6G and Supplementary Figure 2). In conclusion, by identifying MSC subpopulations and uncovering the functional roles of their subpopulations, we found that MSCs and oxidative stress may play an important role in BM.

Figure 6. MSC subpopulations in patients with bone metastases from breast cancer. (A) Single-cell atlas showing MSC subpopulations. (B) Single-cell atlas showing MSC subpopulations in control and breast cancer bone metastasis patients. (C) Differential abundance of MSC subpopulations in control and breast cancer bone metastasis patients. (D) Marker genes specifically and highly expressed in subpopulations of MSC subpopulations. (E, F) Biological processes (E) and signaling pathways (F) enriched in MSC subpopulations. (G) Comparison of the differences in oxidative stress levels in MSC between control and bone metastasis samples.

Clonal evolution of MSCs with bone metastases from breast cancer

Further exploring the differentiation trajectory of the MSC subpopulations, the MSC_MARCKSL1 subpopulation was located at an early position in the development of differentiation and had a high ability to differentiate further toward subpopulations specifically present in the BM (Figure 7A, 7B). Confirming the accuracy of the results in Figure 6C, the MSC_MARCKSL1 subpopulation was significantly more abundant in the primary lesions. By GRN, we explored the TFs regulating MSC subpopulations that are regulated by different TFs, including HOXB3, PPARG, SPI1, XBP1, ELF5, and ZBTB7B, respectively (Figure 7C, 7D). Taken together, we can know the differentiated developmental trajectory of MSC subpopulations and explored their transcriptional regulation.

Figure 7. Clonal evolution of MSC subpopulations in patients with bone metastases from breast cancer. (A, B) Pseudo-time values (A) and developmental trajectories (B) of MSC subpopulations, pie charts representing the proportion of control and breast cancer bone metastasis patients in MSC subpopulations. (C) Co-expression modules of transcription factors in MSC subpopulations of patients with breast cancer bone metastases. Left: Identification of regulator modules based on the regulator’s linkage specificity index matrix. Middle: representative transcription factors and their binding patterns in the modules. Right panel: cellular subpopulations in which transcription factors are located. (D) Single-cell atlas showing transcription factors regulating MSC subpopulations.

Intercellular communication in breast cancer bone metastases

Since we have successfully outlined cellular-level alterations in BM and characterized the functional and transcriptional profiles of different cellular subpopulations, we used a public ligand-receptor database to infer intercellular communication during BM. By comparing cell identity-specific genes with ligand receptors, we classified hypothetical ligand-receptor pairs for different cell populations in control and BM samples. For controls, we found a stronger interaction between fibroblasts and BC cells (Figure 8A). Notably, in BM, BC cells showed the most interactions with other cell types (Figure 8B and Supplementary Table 3), and we observed strong ligand-receptor pairs in BC cell subpopulations.

Figure 8. Intercellular communication in breast cancer bone metastases. (A) Intercellular communication of cell subpopulations in control samples. (B) Intercellular communication of cell subpopulations in bone metastasis samples.