Single-Cell Transcriptomic Analysis Identifies Senescent Osteocytes That Trigger Bone Destruction in Breast Cancer Metastasis

Bone metastases represent an advanced stage of breast cancer, marked by malignant cells that escape the primary breast tumor and colonize the bone tissue (1). Skeletal metastases affect the overall course of the disease and are a major cause of morbidity, severe pain, impaired motility, pathologic fractures, and mortality. Existing studies on skeletal metastasis have focused on how cancer cells benefit from the bone microenvironment by establishing a vicious cycle in which breast cancer cells reprogram osteoblasts (OB) or osteoclasts to cause aberrant bone destruction or bone formation, which in turn provides tumor cells with bone-derived growth factors that fuel tumor growth in the bone (2). Advances in the last decade have provided a deeper understanding of the complexity of the tumor microenvironment in bone and led to the identification of new cellular interactions between cancer and bone cells contributing to cancer progression beyond those originally described in the “vicious cycle” paradigm (1). Thus, understanding the interplay between cancer cells and other tumor microenvironment cells is critical for identifying new targets that interfere with the progression of tumors in bone.

Osteocytes (Ot), the primary resident cells of the bone tissue, have traditionally been recognized for their role in maintaining bone health and homeostasis (3) but have been overlooked in bone cancer research for decades. Recent evidence supports the idea that Ots play a key role in the complex interplay between cancer and the skeletal system. Preclinical and clinical studies in multiple myeloma and metastatic prostate cancer have demonstrated that Ots influence tumor growth, bone destruction, and even the efficacy of cancer therapies (4–8). Further, the knowledge acquired on the role of Ots and their derived factors in these cancers has provided the rationale for developing novel therapeutic interventions to treat cancer in the bone (5, 9–11). Mounting evidence suggests that Ots interact with breast cancer cells and influence their proliferation, migration, and invasion abilities (12–15). However, how cancer cells affect Ots in the tumor microenvironment is largely unknown.

In this study, we used single-cell RNA sequencing (scRNA-seq) and a combination of murine in vitro, ex vivo, and immunocompetent preclinical models to explore the impact of established breast cancer bone metastases on Ots. Further, we validated our findings with human breast cancer cells in immunodeficient models, human bone tissue, and samples from patients with established breast cancer bone metastases. Our study shows that established breast cancer bone metastases induce premature senescence in Ots. Moreover, our results reveal that senescent Ots acquire a pro-osteoclastogenic senescence-associated secretory phenotype (SASP). Furthermore, we provide evidence that depleting senescent cells using senolytics may represent an attractive adjuvant therapy to blunt bone loss in bones with established breast cancer metastasis.

Liberase was purchased from Roche. Dulbecco's phosphate-buffered saline, Hank’s Balanced Salt Solution, 4′,6-diamidino-2-phenylindole (DAPI), and TRIzol were purchased from Thermo Fisher Scientific. Dasatinib (cat. #D-3307) was ordered from LC Laboratories and quercetin (cat. #Q4951) from Sigma-Aldrich. Doxorubicin was ordered from Tocris Bioscience (cat. #2252). Annexin V/propidium iodide Apoptosis Kit was obtained from BD Biosciences. DMEM, α-Minimum Essential Media (αMEM), FBS, bovine calf serum, antibiotics (penicillin/streptomycin), and TRIZol were purchased from Life Technologies. Trypan Blue was purchased from Sigma-Aldrich. Normocin, plasmocin, and puromycin (cat. #58-58-2) were obtained from InvivoGen. D-Luciferin was obtained from PerkinElmer and coelenterazine from Nanolight Technology. Anti-GAPDH (cat. #2118S, RRID: AB_561053) and anti-P16 (cat. #ab211542, RRID: AB_2891084) were purchased from Cell Signaling Technologies and Abcam, respectively.

Murine EO771 (RRID: CVCL_GR23) and EO771-luciferase (EO771-luc)-expressing mammary cancer cells were provided by Dr. Karbassi (University of Arkansas for Medical Sciences) and Dr. Norian (University of Alabama at Birmingham, Birmingham, AL), respectively. Human MDA-MB-231 (RRID: CVCL_0062) breast adenocarcinoma cell line was purchased from ATCC. MDA-MB-231-luciferase cells for ex vivo experiments were generated by transducing cells with lentiviral particles carrying a nonsecreted Gaussia luciferase vector purchased from BPS Bioscience (cat. #79893-C). MDA-MB-231-luciferase cells for in vivo experiments were provided by Dr. Zhang (Baylor College of Medicine, Houston, TX). Human MCF10 epithelial cells (RRID: CVCL_0598) and MCF7 (RRID: CVCL_0031) breast cancer cells were provided by Dr. Steven Post and Dr. Katie R. Ryan (University of Arkansas for Medical Sciences, Little Rock, AR), respectively. All breast cancer cells were cultured in DMEM with 10% FBS, 1% penicillin and streptomycin, 0.2% normocin, 1.5 mg/mL sodium bicarbonate, and 2% HEPES buffer. Murine Ocy454 (RRID: CVCL_UW31) Ot-like cells were provided by Dr. Pajevic (Boston University, Boston, MA; ref. 16) and cultured in αMEM medium with 10% FBS, 1% penicillin and streptomycin, and 0.2% normocin in rat type I collagen–coated flasks. Prior to performing the experiments, Ocy454 cells were cultured at 37°C for 2 weeks. Murine MLOA-5 (RRID: CVCL_0P24) and MLOY-4 (RRID: CVCL_M098) Ot-like cells were obtained from Kerafast and cultured in 2.5% FBS and 2.5% BCS with 1% penicillin and streptomycin and 0.2% normocin. MLO-Y4-GFP cells have been described previously (17). Cell lines were routinely assessed for Mycoplasma and authenticated by morphology, gene expression profiles, and tumorigenic capacity. Cell culture studies were performed by (i) treating Ots with conditioned media (CM; 50%) from breast cancer cells or (ii) coculturing breast cancer and Ot cells in a cell-to-cell manner (1:1) for 24 hours. Breast cancer and osteocytic CM were prepared by culturing 2 × 10⁶ cells in 10 mL of culture medium for 48 hours. Ot-like cultures were treated with dasatinib (200 nmol/L) and quercetin (50 µmol/L) after 48 hours of incubation with the breast cancer CM.

Single-cell studies were performed in NuTRAP^−/+; DMP1-8kb-Cre^−/+ reporter mice were generated by crossing B6;129S6-Gt(ROSA)26Sortm2(CAG-NuTRAP)Evdr/J mice (NuTRAP; #029899; Jackson Laboratory; ref. 18) with DMP1-8kb-Cre mice (19). Seven-week-old NuTRAP^−/+; DMP1-8kb-Cre^−/+ female mice were inoculated intratibially with 10⁵ EO771-luc cells or PBS as a control and were sacrificed after 14 days. To assess the impact of established breast cancer bone metastases on cellular senescence, C57BL/6 female mice were injected intratibially with 10⁵ murine EO771-luc cells or saline, and 3 days later randomized by body weight to the following groups: (i) naïve mice orally receiving vehicle (10% EtOH, 30% PEG, 60% Phosal-50 PG); (ii) EO771-luc–bearing mice orally administered vehicle; or (iii) EO771-luc–bearing mice orally receiving a senolytic cocktail (DQ) of dasatinib (5 mg/kg) and quercetin (50 mg/kg) once a week. The C57BL/6 model injected with murine EO771 was used to avoid confounding factors such as those arising from the immune system or cross-species differences. To assess the impact of human breast cancer cells on cellular senescence, (i) 8-week-old SCID female mice were injected intracaudally (modeling late stages of cancer metastases) with 10⁵ human MDA-MB-231-luc cells or PBS as control and sacrificed after 4 weeks or (ii) 7-week-old NOD.Cg-Prkdc^scidIl2rg^tm1Wjl/SzJ (Jackson Laboratory, strain, #05557) mice were injected intratibially (modeling established breast cancer bone metastases) with 10⁵ human MDA-MB-231 breast cancer cells (left limb) or saline (right limb). The sample size for these studies was calculated based on previous studies (5, 9). Mice were housed in ventilated cages and maintained in a pathogen-free, accredited facility under a 12-hour light–dark cycle at constant temperature (23°C) and access to food and water ad libitum.

Cells were isolated from the tibias of NuTRAP^−/+; DMP1-8kb-Cre^−/+ mice after removing the muscle, periosteum, and epiphyses 2 weeks after saline/breast cancer cell injection. Next, we flushed out the bone marrow and performed serial Liberase/EDTA digestions, as described previously (20). Cells from the digestions were pooled, and GFP⁺ cells were sorted using FACS. Cells per condition were encapsulated using a Chromium Controller (10x Genomics), and libraries were constructed using a Chromium Single Cell 3′ Reagent Kit (10x Genomics) by the UAMS Genomics Core. Libraries were sequenced using an Illumina NovaSeq 600 machine to generate FASTQ files. Three independent samples were sequenced from the breast cancer group (Supplementary Fig. S1) and pooled for the final bioinformatics analysis.

The FASTQ files were preprocessed using Cell Ranger software version 6 (10x Genomics) to produce feature barcode matrices. The alignments were performed using the mouse reference genome mm10 and imported for further analysis into the R suite software environment using the Seurat package v4.2.0 (21–23). Quality control protocols were applied to remove outlier barcodes based on the depth, number of genes, and proportion of mitochondrial genes. The harmonization/integration of different samples was performed using the reciprocal principal component analysis method based on dimensional reduction using Uniform Manifold Approximation and Projection (UMAP; ref. 22). Subpopulation identification and clustering were performed using the Louvain algorithm with multilevel refinement (24). Gene-specific markers of individual clusters were identified using the FindMarkersAll function using the MAST algorithm for cell type identification (25). Functional pathway enrichment analysis was performed using PIANO software (26). Senescence and SASP gene ontology (GO) terms were created by compiling publicly available datasets of differentially regulated genes in senescent cells (27, 28). The senescence gene sets from refs. 27, 28 were used to calculate gene signature scores using PIANO (25). Z-score was used to calculate statistical differences. The polar figure was generated using publicly available datasets from Agoro and colleagues (29), Wang and colleagues (30), and Youlten and colleagues (31).

Ex vivo murine bone cultures were established using femurs from C57BL/6 female mice (EO771 cells) or NOD.Cg-Prkdcscid Il2rgtm1Wjl/SzJ (NSG; #005557; Jackson Laboratory) female mice (MDA-MB-231). We (i) injected 10⁵ EO771 breast cancer cells into femoral bones, (ii) treated femurs with 50% CM from breast cancer cells in the presence/absence of dasatinib (400 nmol/L) and quercetin (100 µmol/L), or (iii) injected them with 10⁵ MDA-MB-231 cells and treated them in the presence/absence of doxorubicin (750 nmol/L) with/without dasatinib (400 nmol/L) and quercetin (100 µmol/L) for up to 6 days. Treatments were added on day 1 and refreshed on day 3. Ex vivo human bone organ cultures were established using human cancellous bone fragments similar in size obtained from the femoral head discarded after hip arthroplasty. Bone samples were obtained from two females with no pathologies or medications that could affect bone mass or architecture. For these ex vivo cultures, we (i) plated 2 × 10⁵ MDA-MB-231 breast cancer cells on human bones or (ii) treated them with 50% MDA-MB-231 CM for 5 days. The CM was refreshed every 72 hours. Tumor bioluminescence was imaged 10 minutes after incubating bones with D-luciferin (150 μg/mL) or coelenterazine (100 μmol/L) using an IVIS Lumina XRMS system (PerkinElmer).

Cell proliferation, death, and apoptosis were estimated using the Trypan Blue exclusion method as previously described, by flow cytometry using the Annexin V Apoptosis Detection Kit following the manufacturer’s recommendations, or chromatin condensation and nuclear fragmentation of cells transfected with nuclear green fluorescent protein (4, 6). To block apoptosis, Ot-like cells were pretreated with 50 nmol/L of the caspase inhibitor DEVD (Sigma-Aldrich) 1 hour before the addition of breast cancer CM. DEVD was refreshed every 24 hours.

Ocy454 Ot-like cells were cultured with 50% EO771-CM for 9 days. The cells were washed twice with PBS, and fresh culture medium was added. After 48 hours, CM was harvested from control and senescent Ocy454 cells. CD11b⁺ mononuclear cells were cultured in αMEM containing 10% FBS and 10 ng/mL of M-CSF for 3 days and then treated with a suboptimal dose of RANK ligand (RANKL; 10 ng/mL) for 4 days in the presence/absence of 25% CM from control and senescent Ots. RANKL and CM were replenished every 2 days. Cells were stained for tartrate-resistant acid phosphatase (TRAP) using a leukocyte acid phosphatase kit (Sigma-Aldrich), and TRAP-positive mononuclear cells and multinuclear cells (≥3 nuclei/cell) were scored as previously described (4).

For in vitro studies, Ot-like cells were treated with 50% of CM from breast cancer cells for 2 to 9 days or etoposide (10 μmol/L) for 4 days. After incubation, cells were washed with PBS, fixed, and stained with senescence-associated beta-galactosidase (SA-β-Gal) staining solution (1 mg/mL X-Gal; cat. #X1220, Teknova, 40 mmol/L citric acid, pH 6.0, 5 mmol/L potassium ferrocyanide, 5 mmol/L potassium ferricyanide, 150 mmol/L NaCl, 2 mmol/L MgCl₂) at 37°C for 16 to 18 hours (32). For in vivo studies, bones were fixed in 4% paraformaldehyde, embedded in gelatine, and cut in 30-μm–thick sections. Senescence β-galactosidase staining was adapted from the manufacturer’s protocol (Cell Signaling, Senescence β-Galactosidase Staining Kit, #9860). Solutions were prepared according to the manufacturer’s instructions. Bone sections were fixed with 1× fixative solution (200 μL/bone) for 15 minutes at room temperature. After washing twice with PBS, 200 μL of β-galactosidase staining solution was added to the bone sections and incubated at 37°C overnight. The sections were washed twice with PBS, counterstained with hematoxylin, and mounted using Aqua/Polymount. The number of positive SA-β-Gal cells was imaged and quantified using a brightfield microscope (EVOS FL Auto) at 20×. Two independent investigators imaged and quantified three random areas from each well in a blinded fashion.

Total RNA was isolated from cells and bone tissues using TRIzol and converted to cDNA (Applied Biosciences) following the manufacturer’s instructions. Gene expression was quantified by quantitative real-time PCR (qPCR) using TaqMan assays from Applied Biosystems, following the manufacturer’s instructions. MLOY-4 GFP⁺ cells were cocultured with EO771 breast cancer cells at a 1:1 ratio for 24 hours and sorted using a FACS sorter (FACSAria III, BD Biosciences). Gene expression levels were calculated using the comparative threshold method and normalized to the housekeeping gene GAPDH (4, 6). Fold changes were calculated using control/vehicle conditions as the reference.

Cell lysates (20 μg) were boiled in the presence of SDS sample buffer (NuPAGE LDS sample buffer; Invitrogen) for 10 minutes and subjected to electrophoresis on 10% SDS-PAGE (Bio-Rad Laboratories). Proteins were transferred to PVDF membranes using a semidry blotter (Bio-Rad) and incubated in a blocking solution (5% nonfat dry milk in TBS containing 0.1% Tween-20) for 1 hour to reduce nonspecific binding. Immunoblots were performed using anti-GAPDH (1:1,000) and p16 antibodies (1:2,000) followed by goat anti-rabbit secondary antibodies, conjugated to horseradish peroxidase (1:5,000) in 5% milk (Santa Cruz Biotechnology). Western blots were developed using an enhanced chemiluminescence detection assay following the manufacturer’s directions (Bio-Rad). Protein bands were quantified using ImageJ.

RNA in situ hybridization was performed using the RNAscope 2.5 HD Reagent Kit-RED assay from Advanced Cell Diagnostics following the manufacturer’s instructions, as previously described (33). The following probes were incubated on paraffin-embedded tissue sections for 2 hours at 40°C: murine p16^Ink4a (cat. #411011) and positive/negative controls (cat. #313911; cat. #310043). The signal was then detected at room temperature for 10 minutes. Sections were counterstained with hematoxylin, dehydrated at 60°C for 20 minutes, and mounted using VectaMount permanent mounting medium (Vector Laboratories). The number of positive and negative Ots was quantified using a brightfield microscope at ×40 magnification. Analyses were performed in the cortical bone of an 800-μm region of the tibia, starting 200 μm below the growth plate, in a blinded fashion by two independent investigators.

Pericentromeric satellite heterochromatin undergoes decondensation and elongation in senescent cells, and the large-scale unraveling of pericentromeric satellite heterochromatin DNA, termed senescence-associated distension of satellites (SADS), is a well-established marker of cell senescence in vivo (34). We used non-decalcified tibiae from naive and breast cancer–bearing mice embedded in methylmethacrylate. Fluorescence in situ hybridization staining was performed as previously described (35). SADS were visualized in Ots (senescent Ots ≥4 SADS per Ot) using c-FISH [Cy3-labeled (F3002)], CENPB-specific (ATT​CGT​TGG​AAA​CGG​GA) peptide nucleic acid (PNA) probe (Panagene Inc.), and quantified by confocal microscopy (Zeiss LSM 880, 100× oil) in the cortical bone starting 200 μm below the growth plate. At least 50 nuclei were analyzed for each sample in a blinded fashion by two independent investigators.

Telomere dysfunction–associated foci (TAF) assays were performed on murine nondecalcified tibias embedded in methyl methacrylate using a protocol established by Farr and Khosla labs (36). TAFs were visualized in Ots (senescent Ots ≥3 TAFs per Ot) using a primary antibody for γH2AX (1:200; anti–γH2A.X rabbit monoclonal antibody, Cell Signaling Technology; 9718) and a Cy3-labeled telomere-specific (CCCTAA) peptide nucleic acid probe (TelC-Cy3, Panagene Inc.; F1002). The mean number of TAF per Ot was quantified by confocal microscopy (Zeiss LSM 880, 63× oil) in the cortical bone starting 200 μm below the growth plate. At least 35 nuclei were analyzed for each sample in a blinded fashion by two independent investigators.

MMP13 immunofluorescence staining was performed on decalcified paraffin–embedded bone tissue sections. Tissue sections were incubated with anti-MMP13 (1:50; Abcam, AB39012) overnight at room temperature, washed, incubated with Goat Anti-Rabbit IgG H&L—Alexa Fluor 594 (1:1,000; Abcam, AB150080) for 1 hour and with copper sulfate (cat. #209198, Sigma-Aldrich) for 10 minutes, and mounted with Prolong Gold Anti-Fade DAPI mounting medium (cat. #P36935, Invitrogen). MMP13 positive/negative Ots were imaged using a Zeiss Axio Imager.M2 image system at ×40 magnification. The percentage of MMP13-positive Ots was assessed using Fiji (37) in a blinded manner by two independent investigators.

Static and dynamic bone histomorphometric analyses were performed using the OsteoMeasure High-Resolution Digital Video System (OsteoMetrics) as previously described (5, 9). Analyses were conducted in the cancellous bone of an 800-μm region of the tibiae, starting 200 μm below the growth plate.

The bone resorption biomarker C-telopeptide of type 1 collagen (CTX; Immunodiagnostic Systems, cat. #AC-06F1) and the bone formation marker propeptide of type 1 collagen (P1NP; Immunodiagnostic Systems, cat. #AC-33F1) were analyzed in the serum from mice or in CM from bones cultured ex vivo, as previously described (5, 9).

Osteolytic lesions were imaged using a Faxitron X-ray radiography system (Hologic) as previously described (5, 9). MicroCT imaging was performed in live mice using a vivaCT 80 (Scanco Medical AG). Analyses were performed at the cancellous bone of the proximal tibia, in an area 20 µm below the growth plate, using 10 μm resolution.

Tumor growth was monitored weekly using the IVIS Lumina XRMS system. Mice were intraperitoneally injected with 150 mg/kg D-luciferin, and luminescence imaging was initiated 10 minutes after luciferin injection.

Archived diagnostic transiliac biopsies collected at the Pathological Biobank at Odense University Hospital from four consented female patients with established breast cancer metastatic bone disease were used for this analysis. The average patient age was 63 years (range, 51–73 years), and all patients received radiotherapy for their primary cancer and zoledronic acid at the time of dissemination to the bone (4–5 mg/year, one dose). None of the patients had received chemotherapy within the last 5 years before dissemination. Three-millimeter bone biopsies were fixed in 4% formalin for 24 hours, decalcified in 10% formic acid for 7 hours, and embedded in paraffin. One 3.5-μm-thick section from each sample was multiplex immunostained for cytokeratin 7 and 19 (CK7 and CK19), along with fluorescent in situ hybridization for SPP1 and CDKN2A expression. Briefly, sections were deparaffinized in a xylene and ethanol gradient and pretreated in Custom Reagent (Advanced Cell Diagnostics) for 20 minutes at 40°C. Sections were then hybridized overnight at 40°C with a channel 1 probe targeting the 6 to 1,500 nucleotide region of SPP1 mRNA (NM_001251829.1) and a channel 2 probe targeting nucleotides 95 to 1,206 of CDKN2A mRNA (NM_000077.4). Signal amplification was performed per the manufacturer’s recommendations and visualized with opal 690 and opal 620 dyes (Akoya Biosciences). Sections were then HRP blocked (Advanced Cell Diagnostics) for 15 minutes at 40°C, followed by 5% casein blocking buffer (20 minutes at room temperature) and incubated with an antibody cocktail against CK7 and CK19 (mouse IgG2a anti-CK19 clone A53-BA2.26, Sigma-Aldrich). Matched negative controls were conducted by omitting the target probe and primary antibody. Slides were counterstained with Hoechst and mounted in Prolong Gold mounting medium before whole-slide scanning using a VS200 Olympus slide scanner. Scanning was performed at ×40 magnification with 30 ms exposure in the DAPI channel (455 nm), 80 ms exposure in the Cy3 channel (565 nm), 300 ms exposure in the Cy5 channel (670 nm), and 70 ms in Texas Red (615 nm). The image visualization settings (brightness and contrast) were adjusted for image quantification and representative image acquisition. Artificial intelligence (AI)–assisted histology was performed using the IF + FISH v2.1.5 module of HALO (Indica Labs). Cell segmentation and classification were followed by 10-µm band proximity analysis of Ots within 500 μm of CK7/19⁺ cancer cells.

Data were analyzed using the GraphPad software (GraphPad Software Inc.). Differences in means were analyzed using a combination of unpaired t test and ANOVA, followed by pairwise multiple comparisons (Tukey). Values are reported as the mean ± SD. Statistical significance was set at P ≤ 0.05. The data analysis was performed in a blinded manner.

All animal procedures were performed in accordance with the guidelines issued by the Institutional Animal Care and Use Committee of the University of Arkansas for Medical Sciences (AUP protocol, #2022200000489) and by the General Administration of the Free State of Bavaria (Regierung von Oberbayern), Germany (ROB-55.2-2532.Vet_02-23-14). Institutional and national guidelines for the care and use of laboratory animals were followed in these studies. The collection and de-identification of human bone samples were coordinated by the UAMS Winthrop P. Rockefeller Cancer Institute Tissue Biorepository and Procurement Service (TBAPS) and approved by the UAMS Institutional Review Board (IRB protocol, #262940). Archived diagnostic transiliac biopsies collected at the Pathological Biobank at Odense University Hospital, from four consented female patients with established breast cancer metastatic bone disease were included in this analysis under approval from the National Committee on Health Research Ethics (S-20180057). All participants provided written informed consent before the study procedures with continuous consent ensured throughout participation.

The scRNA-seq data generated in this study were deposited in the NCBI SRA database under the Bioproject PRJNA1033671. The additional Ot scRNA-seq datasets analyzed in this study were obtained from Gene Expression Omnibus repository at GSE208152 (29) and GSE154719 (30), and from ArrayExpress (https://www.ebi.ac.uk/arrayexpress) under the following accession number: E-MTAB-5533 (31). All other raw data are available upon request from the corresponding author.

To determine the impact of established breast cancer bone metastasis on Ots, we crossed NuTRAP reporter mice (18), a Cre-inducible strain that allows the labeling and simultaneous isolation of cell type–specific nuclei (mCherry) and mRNA (GFP), with Dmp1-8kb-Cre mice to induce active Cre recombination in Ots (38). The mice were intratibially injected with EO771-luc breast cancer cells (Fig. 1A). GFP⁺ cells were isolated 2 weeks after injection, when the mice displayed active tumor growth and established bone disease, and scRNA-seq was performed in the isolated cell population (Fig. 1B and C). Transcriptomic profiling of the GFP⁺ cells identified three distinct clusters: (i) pre-OBs, (ii) OBs, and (iii) Ots (Fig. 1D; Supplementary Fig. S1), as defined by distinct gene expression patterns (Fig. 1E; Supplementary Fig. S2). The Ot fraction represented ∼10% of the total cells and remained unchanged between groups (Fig. 1F and G). In contrast, the OB fraction decreased by 20%, and the pre-OB fraction increased by 23% in the breast cancer versus the control group (Fig. 1F and G). In addition to these three identified clusters, our analysis detected an additional cell type (FiX) characterized by high expression of Acta2, Myh11, Tagln, Rgs5, and Igfbp7, genes expressed by cancer-associated fibroblasts (Supplementary Fig. S3; refs. 39, 40). Intriguingly, this population expanded in bones with breast cancer metastasis (Supplementary Fig. S3).

We first performed GO analysis of the genes differentially expressed in cells isolated from control and breast cancer–bearing bones to identify signaling pathways altered in osteoblastic cells by breast cancer bone metastasis. When combining the three cell populations (pre-OBs + OBs + Ots), we found enrichment in GO terms associated with cellular senescence, SASP, and inflammatory response (Fig. 1H). Furthermore, cells isolated from bones with established breast cancer metastasis exhibited a higher senescence score [constructed with previously reported transcriptional biomarkers of senescence (27, 28)], than those from control bones (Fig. 1I) and upregulation of the SASP-related genes Mmp13, Spp1, Serpine2, Timp2, Igfbp7, Igfbp5, and Vegfa (Fig. 1J).

Next, we looked at each individual population, focusing our analysis first on the Ot population. Our dataset contains one of the largest Ot populations sequenced to date (719 naïve vs. 340 breast cancer Ots). The polar plot comparison displayed in Fig. 2A highlights common genes identified in our dataset and in previous studies, including the traditional markers Dmp1, Phex, Dkk1, and Pdpn, as well as new markers such as Cd44, Sema5a, and Tgfbr2. We identified 850 genes that were differentially expressed in Ots from control and breast cancer–bearing mice (Fig. 2B and C; Supplementary Table S1). In addition, GO term analysis revealed the upregulation of gene sets related to cellular senescence in Ots from bone metastases, whereas GO terms associated with cellular proliferation and cell cycle regulation were downregulated (Fig. 2D). Furthermore, we found an increased prevalence of Ots with a higher senescence score and higher mRNA levels of SASP-related genes Spp1, Mmp13, Cstb, Serpine1, and Bmp2 (Fig. 2E–G). Similar observations were made in the OB cluster, which showed enrichment in senescence GO terms and senescence scores and upregulation of SASP-related factors in bones from mice with breast cancer tumors (Supplementary Fig. S4A–S4E). Although the pre-OB cluster in mice with breast cancer bone tumors also exhibited a higher senescence score and upregulation of SASP (Supplementary Fig. S5A–S5E), the small number of pre-OB cells in the naïve mice limited the robustness of this observation. Collectively, these results suggest that established breast cancer metastases generate a bone microenvironment that is conducive to premature cellular senescence in Ots and OBs.

To determine whether breast cancer cells or their products directly promote Ot senescence, we performed a series of in vitro studies using Ot-like cell lines and ex vivo cultures containing human and murine primary Ots. CM from murine EO771 or human MDA-MB-231 breast cancer cells markedly decreased the cell number but only marginally increased apoptosis in Ot-like cells (Supplementary Fig. S6A–S6E). Inhibition of apoptosis using the caspase3 inhibitor DEVD prevented apoptosis induced by EO771 cells but did not alter the number of live Ots. Together, these findings suggest that breast cancer cells induce proliferative arrest in Ot-like cells, a feature of senescence (Supplementary Fig. S6F). Further, Ot-like MLO cell lines cultured in direct contact with breast cancer cells (24 hours) or treated with CM from breast cancer cells (48 hours) also showed upregulation of senescence-related genes p16^Ink4a, p21^Cip1, Mmp13, and Il6 (Supplementary Fig. S7A–S7C). Because MLO Ot-like cells are immortalized, they are not an optimal model for studying senescence. Thus, we used Ocy454 cells (16), which are conditionally immortalized when cultured at 33°C, for subsequent in vitro studies. After 2 weeks at 37°C, Ocy454 Ot-like cells exhibited morphological features and a gene expression profile consistent with mature Ots (Supplementary Fig. S8). Ocy454 cells treated with CM from EO771 (Fig. 3A–D) or MDA-231 (Fig. 3E–H) breast cancer cells exhibited upregulation of the senescent markers p16^Ink4a (RNA and protein) and p21^Cip1 and the SASP-related factors Il6 and Mmp13 compared to control Ocy454 cells, high SA-β-Gal activity, and increased prevalence of SA-β-Gal⁺ cells. Remarkably, these features of senescence were evident after 48 hours of treatment with CM from breast cancer cells (Supplementary Fig. S9A–S9F). Comparable results were obtained with CM from MCF7 cells, an ER⁺ human breast cancer cell line, but not with noncancerous breast epithelial MCF10 cells (Supplementary Figs. S10A–S10D and S11A–S11C).

To study the response of primary Ots to breast cancer cells, we next used ex vivo bone cultures (Fig. 4A), a system that recapitulates the spatial dimension, cellular diversity, and molecular networks of the tumor niche in a controlled setting (41). Treatment of bones with EO771-CM for 2 or 5 days increased the prevalence of p16^Ink4a+ primary Ots by 30% versus controls (Fig. 4B and C) and upregulated the expression of p16^Ink4a, p21^Cip1, Mmp13, Spp1, Il6, and Mmp9 (Supplementary Fig. 12A). Comparable results were obtained with CM from human MDA-MB-231 breast cancer cells (Supplementary Fig. S12B). Treatment with the senolytic drugs DQ prevented or reduced the increased expression of p16^Ink4a, Mmp13, and Il6 (Fig. 4D), supporting a direct link between changes in gene expression and cellular senescence induced by breast cancer cells. Bones from mice with established human MDA-MB-231 breast cancer metastasis exhibited a higher prevalence of TAF⁺ senescent Ots compared to the saline-injected contralateral leg (Fig. 5A–C). Similarly, we detected senescence-associated β-Gal⁺ osteoblastic cells (Ots and OBs) in bones colonized by metastatic breast cancer cells after intracaudal inoculation with human MDA-MB-231 breast cancer cells (Supplementary Fig. 13A). This was accompanied with an increased expression of SASP factors Mmp13, IL6, and Rankl in bones bearing MDA-MB-231 tumors (Supplementary Fig. S13B). We also explored whether human breast cancer cells induce senescence in human Ots. MDA-231-CM upregulated P16^Ink4a, P21^Cip1, MMP13, SPP1, and IL6 gene expression compared to the controls in human bones containing primary Ots (Fig. 5D). Moreover, we allowed the breast cancer cells to infiltrate human bones cultured ex vivo. We detected active tumor growth after 2 days and a similar upregulation of senescence-related markers in bones infiltrated with breast cancer cells compared to control bones, except P16, which remained unchanged (Fig. 5E). Finally, we examined Ot senescence in bone biopsies from a small cohort of breast cancer patients with established bone metastasis. We detected P16^Ink4a+ and SPP1⁺ Ots (Fig. 5F), which were preferentially located close to the bone marrow areas infiltrated by breast cancer cells (Fig. 5G). The percentage of P16^Ink4a+SPP1⁺ Ots showed a nonsignificant positive correlation with the tumor burden in the bone marrow (Supplementary Fig. S14). These findings, together with our bioinformatics results, demonstrate the induction of premature cellular senescence and SASP development in Ots in established breast cancer bone metastases.

Because the accumulation of senescent Ots has been linked to bone loss seen with aging or radiation therapy (42, 43), we hypothesized that the accelerated cellular senescence induced by breast cancer cells in the bone contributes to bone destruction. To test this hypothesis, we treated mice with DQ, a cocktail of senolytic drugs known to deplete senescent bone cells (Fig. 6A; refs. 42, 43), 3 days after injecting cancer cells intratibially. DQ therapy did not affect tumor progression (Fig. 6B). Bones bearing breast cancer cells exhibited a higher prevalence of p16^Ink4a+, SADS⁺, and MMP13⁺ Ots (Fig. 6C–E). Treatment with DQ prevented the increase in senescent SADS⁺ Ots and attenuated the increase in MMP13⁺ Ots in mice with bone metastasis (Fig. 6C–E). In addition, mice with established breast cancer bone metastasis treated with DQ had fewer osteolytic lesions (Fig. 7A), higher bone mass, and improved bone microarchitecture than vehicle-treated mice with bone metastasis (Fig. 7B and C). At the end of the study, due to the aggressiveness of the model, we could not detect differences in CTX or bone formation rate between mice with bone metastasis receiving vehicle or DQ (Supplementary Fig. S15A–S15D), suggesting that the protective effects of DQ occurred during the initial stages of tumor progression. To assess this possibility, we developed an ex vivo model resembling in vivo conditions in which tibiae from female mice were injected with EO771-luc cells and treated with vehicle or DQ (Fig. 7D). Similar to the in vivo study, bones bearing breast cancer tumors exhibited p16^Ink4a mRNA upregulation, higher CTX, and lower P1NP in the culture media (Fig. 7E–G). Treatment with DQ did not affect tumor burden or P1NP levels; however, it restored p16^Ink4a expression to control levels and decreased by ∼50% CTX levels in bones bearing breast cancer cells (Fig. 7F and G). These data suggest that breast cancer-induced senescence in the tumor niche is an initiating event that triggers bone loss during breast cancer bone metastasis by promoting bone resorption. Lastly, we evaluated the potential impact of DQ on the antitumor efficacy of chemotherapy on established breast cancer bone metastasis using the same ex vivo model. Treatment with doxorubicin decreased tumor burden by 50% in bones colonized by human MDA-MB-231-luciferase breast cancer cells, an effect that was not impacted by DQ coadministration (Supplementary Fig. S16A–S16C).

Further bioinformatic analyses of our scRNA-seq dataset revealed the presence of a higher number of p16^Ink4a+-Rankl⁺ cells in bones with established breast cancer metastasis, indicating the presence of senescent cells of the osteoblastic lineage with a proresorptive phenotype (Fig. 8A). Moreover, comparative transcriptomic analysis identified that Ots from bones with breast cancer had a gene signature associated with osteoclastogenesis, with changes in the expression of pro- and anti-osteoclastogenic factors Rankl, Mmp13, Lgals3, Serpine3, Cd9, Vegfa, and Cthrc1 (Fig. 8B). On the basis of these observations, we investigated the effect of breast cancer cells on the pro-resorptive potential of Ot-like cell lines. Treatment with breast cancer CM, but not from noncancerous breast epithelial cells, upregulated Rankl in Ot-like cells (Fig. 8C–E; Supplementary Fig. S17A–S17D). Breast cancer CM also increased Mmp13 and Il6 and decreased Cthrc1 mRNA expression in MLO and Ocy454 cells (Fig. 8F) Treatment with DQ prevented Rankl and Il6 upregulation and reduced by 80% the elevated Mmp13 in Ocy454 cells treated with EO771-CM (Fig. 8F). Lastly, we investigated if premature cellular senescence in Ots contributes to osteoclastogenesis by assessing in vitro if factors derived from breast cancer-induced senescent Ots affect osteoclast precursor differentiation in vitro. We found that senescent Ot-CM led to the formation of more osteoclasts than CM from the control Ocy454 cells (Fig. 8G). This set of experiments demonstrated that upon premature cellular senescence induced by breast cancer cells, Ots acquire a unique pro-osteoclastogenic SASP, which acts on osteoclast lineage cells in a paracrine manner to support osteoclastogenesis and bone destruction.

In this study, we mapped the transcriptome of Ots from normal bones and bones with established breast cancer metastasis and defined a signature of 850 genes that distinguishes healthy from diseased Ots. We found that the transcriptome of Ots from bones bearing breast cancer cells was significantly enriched in genes associated with cellular senescence. Further, we demonstrated that breast cancer cells induce premature cellular senescence and a distinctive pro-osteoclastogenic SASP in Ots, using various mouse models and histological sections from breast cancer patients with bone metastasis. Moreover, we showed that senescent Ots promote osteoclastogenesis and that pharmacological depletion of senescent cells mitigates bone resorption and preserves bone mass and microarchitecture in a mouse model of established breast cancer skeletal metastasis. Collectively, these findings reveal the profound cellular and molecular reprogramming that Ots undergo in established breast cancer metastasis and illustrate the contribution of Ots to the lytic bone disease induced by breast cancer cells.

scRNA-seq has emerged as a transformative tool in bone research, enabling a detailed exploration of gene expression profiles at the individual cell level. However, previous studies failed to capture a significant number of Ots. Our study is one of the first to leverage scRNA-seq to study the Ot transcriptome in vivo in the dynamic context of established breast cancer bone metastasis. Our approach using a NuTRAP reporter mouse, a novel tool in bone cell studies, crossed with the Dmp1-8kb-Cre strain allowed us to isolate and sequence a significantly higher number of Ots compared to other methods (TdTomato or microbeads; refs. 29, 30, 44). Through comparative analysis with existing datasets, this study began to delineate a set of common markers defining primary Ots under physiological conditions, which goes beyond the more traditional markers previously described (45). In addition, our findings revealed a distinctive gene signature of diseased Ots specific to breast cancer bone metastasis, which is a critical step in understanding the molecular intricacies involved in the reprogramming of Ots by breast cancer cells. Our method also captured OB and pre-OB populations and revealed an accumulation of OB precursors accompanied by a reduction in OBs in bones with breast cancer metastasis. These in vivo results corroborate previous in vitro findings showing that breast cancer cells suppress the differentiation of OB precursors (11, 46–49), which might accumulate in the bone marrow niche over time. Bioinformatic analysis of our dataset also unveiled a new cell population that descended from Dmp1⁺ cells characterized by high expression of genes associated with cancer-associated fibroblasts and increased in bones colonized by breast cancer cells. Future studies beyond the scope of the current work are warranted to investigate the biological implications of this cell population.

The current study showed that breast cancer cells in the bone promote a rapid increase in senescent Ots. Cellular senescence is a cell fate that involves irreversible proliferative arrest, altered chromatin organization, and resistance to apoptosis. A rigorous orthogonal approach was used to confirm the presence of Ot senescence. First, we generated a bioinformatic senescence score algorithm, allowing for a more comprehensive assessment of senescent cells at the transcriptome level. This approach bypasses the limitations of relying solely on the expression of p16^Ink4a transcript levels to identify senescent cells because of its low and variable expression in senescent cells (27). Second, we confirmed the utility of our bioinformatics tool in predicting senescence using in vitro, in vivo, and novel ex vivo models of established breast cancer in bone. Consistent with the elevated senescence score in Ots from breast cancer bone metastasis, we detected an accumulation of senescent Ots using a combination of gold standard methods to detect senescence: SA-β-Gal staining, in situ hybridization (p16^Ink4a), immunostaining (SASP), and FISH (SADS/TAF). Lastly, we used AI-assisted histological phenotyping to confirm the presence of senescent Ots close to metastatic breast cancer cells in patients. This combination of in silico, preclinical, and clinical data shows that our bioinformatics senescence score is a powerful tool for assessing cellular senescence in transcriptomic datasets from bone and demonstrates that Ot cellular senescence is prematurely induced by breast cancer cells in the tumor niche.

Accumulation of senescent Ots has been described in various models of bone loss induced by aging, radiotherapy, or diabetes (36, 42, 43). In these models, senolytics or Ot-specific genetic depletion of senescent cells preserve bone mass by suppressing bone resorption and maintaining or increasing bone formation (43, 50). Our study shows that senolytic therapy also protects against bone loss in established breast cancer metastasis. Our in vivo and in vitro studies support that, in the context of breast cancer, senescent Ots have increased osteoclastogenic potential (increased Rankl, Mmp13, Il6) and are key drivers of cancer-induced bone loss. This observation coincides with a previous study showing that senescent osteoblastic cells have increased osteoclastogenic potential and overexpress RANKL during skeletal aging (32). While our research demonstrates that senescence plays a significant role in driving these transcriptional changes, additional factors, including parathyroid hormone-releasing protein produced by breast cancer cells, are likely to play a contributory role in enhancing the osteoclastogenic potential of Ots irrespective of cellular senescence (51). The bone protection seen with DQ could be explained by a transient decrease in bone resorption, as detected in our ex vivo system, which, in our experience, often predicts in vivo responses (41, 52). However, the ex vivo results should be interpreted in the context of different drug delivery dynamics and lack of hormonal systemic influences compared to in vivo models. In contrast, senolytic therapy did not affect bone formation in vivo or ex vivo, suggesting that mechanisms other than cellular senescence are responsible for reduced OB function caused by breast cancer. Notably, global genetic clearance of senescent cells results in greater benefits than Ot-specific clearance during aging, suggesting the involvement of additional senescent cell populations in bone preservation (50). Thus, it is likely that the clearance of other senescent cell populations (OBs, immune cells, or pre-OBs) contributes to the reduced bone loss seen with senolytics. Further studies are warranted to determine the specific contribution of senescent Ots versus other senescent cells in the tumor niche.

Many chemotherapeutic interventions trigger senescence in cancer cells, resulting in stable cell cycle arrest and reduced tumor growth. However, this process eventually leads to SASP induction and the creation of a pro-inflammatory and immunosuppressive microenvironment that can support tumor progression (53–55). Senolytics have antitumor efficacy in cancer cell lines when combined with various senescence-inducing chemotherapies (53–55). Further, high doses of senolytic agents alone can decrease breast cancer cell proliferation and tumor growth, although the results appear to be concentration- and cell-dependent (56–58). Moreover, genetic or pharmacological induction of senescence in osteoblastic cells, stromal cells, or fibroblasts within the tumor microenvironment has been shown to promote breast cancer growth and bone disease (59, 60). This body of work suggests that inhibition of cellular senescence could affect tumor growth via direct or indirect mechanisms. While we cannot discard that the senolytics led to early decreases in tumor burden not captured in our study, DQ did not significantly affect overall tumor progression in our model. A similar observation was recently reported using the senolytic ABT-263 (61). The absence of an overall antitumor effect can be attributed to either the dose of senolytic therapy used in our studies, the early but temporary inhibition of bone resorption by DQ, or the inability of malignant tumors to undergo spontaneous senescence without external interventions, such as chemotherapy or radiation therapy. Therefore, additional research is needed to determine if senescent Ots affect tumor progression in bone. Additionally, while our ex vivo studies suggest that short-term coadministration of senolytics does not compromise the anticancer efficacy of doxorubicin in established bone metastasis, we acknowledge that future studies are warranted to explore the long-term implications of such a combination or a one-two-punch sequential treatment approach employing senescence-induced chemotherapy followed by senolytic therapy in both tumor and bone health.

In summary, by combining transcriptomic, bioinformatics, and pharmacological approaches, we identified that established breast cancer bone metastases transcriptionally reprogram and induce premature senescence in Ots. The results of our study extend beyond the existing work focusing on the interactions between breast cancer cells and OBs/osteoclasts by shedding light on the potential role of senescent Ots as mediators of bone resorption in breast cancer–induced bone disease. In addition, this work identifies cancer-induced senescence as an initiating event in bone metastases and underscores the therapeutic potential of targeting senescent cells, including Ots, in the tumor niche to mitigate bone loss in patients with cancer with established bone metastasis.

M. Diaz-delCastillo reports grants from NIH during the conduct of the study and grants from Takeda outside the submitted work. C.A. O’Brien reports grants from NIH during the conduct of the study. No disclosures were reported by the other authors.

J. Kaur: Data curation, formal analysis, investigation, methodology, writing–original draft, writing–review and editing. M. Adhikari: Data curation, formal analysis, investigation, methodology, writing–original draft, writing–review and editing. H.M. Sabol: Formal analysis, validation, investigation, methodology, writing–review and editing. A. Anloague: Formal analysis, validation, investigation, methodology, writing–review and editing. S. Khan: Formal analysis, validation, investigation, methodology, writing–review and editing. N. Kurihara: Data curation, formal analysis, methodology. M. Diaz-delCastillo: Resources, data curation, formal analysis, investigation, methodology, writing–review and editing. C.M. Andreasen: Resources, data curation, formal analysis, investigation, methodology, writing–review and editing. C.L. Barnes: Formal analysis, investigation, writing–review and editing. J.B. Stambough: Formal analysis, investigation, writing–review and editing. M. Palmieri: Investigation, methodology, writing–review and editing. O. Reyes-Castro: Formal analysis, investigation, methodology, writing–review and editing. J. Zarrer: Investigation, methodology, writing–review and editing. H. Taipaleenmäki: Investigation, methodology, writing–review and editing. E. Ambrogini: Investigation, writing–review and editing. M. Almeida: Formal analysis, supervision, writing–review and editing. C.A. O’Brien: Resources, formal analysis, supervision, writing–review and editing. I. Nookaew: Data curation, formal analysis, investigation, visualization, methodology, writing–review and editing. J. Delgado-Calle: Conceptualization, resources, formal analysis, supervision, funding acquisition, writing–original draft, project administration, writing–review and editing.

This work was supported by the National Institutes of Health (NIH) R37CA251763, R01CA209882, R01CA241677 to J. Delgado-Calle, P20GM125503 to C.A. O’Brien, F31CA284655 to H.M. Sabol, AG075227 to M. Diaz-delCastillo, and by German Research Foundation (DFG) TA1154/1-2 and TA1154/2-2 to H. Taipaleenmäki, the UAMS Winthrop P. Rockefeller Cancer Institute Seeds of Science Award, Voucher Program Award, and Arkansas Breast Cancer Research Program Award to J. Delgado-Calle. The authors would like to acknowledge the services provided by the TBAPS of the UAMS Winthrop P. Rockefeller Cancer Institute, the Flow Cytometry Core (supported in part by the Center for Microbial Pathogenesis and Host Inflammatory Responses NIGMS P20GM103625), Genomics, Histology, and Microscopy Cores at UAMS. The authors thank Drs. Khosla and Farr (Mayo Clinic, Rochester, MN) and Dr. Dole (UAMS) for providing technical support in the analysis of SASD/TAF and MMP13 using immunofluorescence, respectively, and Kaja Laursen (Aarhus University) and Malene H. Nielsen (University of Southern Denmark) for providing technical support in the analysis of human bone histological sections.