Single-cell RNA landscape of osteoimmune microenvironment in osteoporotic vertebral compression fracture and Kümmell's disease

Abstract

Background: Single-cell RNA sequencing (scRNA-seq) enables specific analysis of cell populations at single-cell resolution; however, there is still a lack of single-cell-level studies to characterize the dynamic and complex interactions between osteoporotic vertebral compression fractures (OVCFs) and Kümmell's disease (KD) in the osteoimmune microenvironment. In this study, we used scRNA-seq analysis to investigate the osteoimmune microenvironment and cellular composition in OVCFs and KD. Methods: ScRNA-seq was used to perform analysis of fractured vertebral bone tissues from one OVCF and one KD patients, and a total of 8,741 single cells were captured for single-cell transcriptomic analysis. The cellularity of human vertebral bone tissue was further analyzed using uniform manifold approximation and projection. Pseudo-time analysis and gene enrichment analysis revealed the biological function of cell fate and its counterparts. CellphoneDB was used to identify the interactions between bone cells and immune cells in the osteoimmune microenvironment of human vertebral bone tissue and their potential functions. Results: A cellular profile of the osteoimmune microenvironment of human vertebral bone tissue was established, including mesenchymal stem cells (MSCs), pericytes, myofibroblasts, fibroblasts, chondrocytes, endothelial cells (ECs), granulocytes, monocytes, T cells, B cells, plasma cells, mast cells, and early erythrocytes. MSCs play an immunoregulatory function and mediate osteogenic differentiation and cell proliferation. The differentiation trajectory of osteoclasts in human vertebral bone tissue was also revealed. In addition, ECs actively participate in inflammatory infiltration and coupling with bone cells. T and B cells actively participate in regulating bone homeostasis. Finally, by identifying the interaction of ligand-receptor pairs, we found that immune cells and osteoclasts have bidirectional regulatory characteristics, have the effects of regulating bone resorption by osteoclasts and promoting bone formation, and are essential for bone homeostasis. It is also highlighted that CD8-TEM cells and osteoclasts might crosstalk via CD160-TNFRSF14 ligand-receptor interaction. Conclusion: Our analysis reveals a differential landscape of molecular pathways, population composition, and cell-cell interactions during OVCF development into KD. OVCFs exhibit a higher osteogenic differentiation capacity, owing to abundant immune cells. Conversely, KD results in greater bone resorption than bone formation due to depletion of MSCs and a relatively suppressed immune system, and this immune imbalance eventually leads to vertebral avascular necrosis. The site of action between immune cells and osteoclasts is expected to be a new therapeutic target, and these results may accelerate mechanistic and functional studies of osteoimmune cell types and specific gene action in vertebral avascular necrosis and pathological bone loss diseases, paving the way for drug discovery.

Keywords: Kümmell’s disease; osteoimmunology; osteoporosis; osteoporotic vertebral compression fracture; single-cell RNA sequencing.

FIGURE 1

Overview of 8,741 single cells from human vertebral bone tissue from KD and OVCF samples. (A) Study overview. (B) Based on uniform manifold approximation and projection (UMAP), dimensionality reduction visualization is performed and coloring is marked according to the cell type. From left to right are 8,741 cells of KD and OVCF samples, 4,489 cells of the KD sample, and 4,252 cells of the OVCF sample. (C) UMAPs are shown in Figure 1B, showing differentially highly expressed classical marker genes for each cell type and coloring according to their expression. (D) Proportional pie chart for each cell type in the respective sample.

FIGURE 2

Distinct subclusters derived from MSC differentiation. (A) All preliminary groups of MSCs, pericytes, myofibroblasts, fibroblasts, and chondrocytes were combined to reconstruct the UMAP profile and labeled for staining according to the cell type. From left to right are the KD and OVCF samples, the KD sample, and the OVCF sample. (B) UMAPs are shown in Figure 2A and show differentially highly expressed classical marker genes in the distinct cell cluster of fibroblasts and colored according to their expression. (C) Proportional pie chart of subclusters of cells derived from the differentiation of MSCs. (D) The upper left panel shows the analysis of single-cell trajectories of all MSCs, OLCs, osteoprogenitor cells, osteoblasts, chondrocytes, pericytes, and myofibroblasts in the KD and OVCF samples; the lower left panel shows the differentiation trajectories of the above cells in KD versus OVCF; and the right panel shows the single-cell trajectories of each cell type. (E) The left panel is a heatmap of branch point 1 shown in Figure 2D, with the GO terms enriched by the pre-branch on the right side of the heatmap and the genes differentially highly expressed based on the GO terms enriched by the pre-branch on the left side of the heatmap. (F) The left panel is a heatmap of branch point 2 shown in Figure 2D, the right panel is the GO terms enriched in cell fate 1 and cell fate 2, and the right side of the heatmap is the genes differentially highly expressed based on the GO terms enriched in cell fate 1 and cell fate 2.

FIGURE 3

Differentiation track of osteoclasts in vertebral bone tissue. (A) Preliminary clusters of monocyte lineages were merged to reconstruct UMAP profiles, with coloration labeled by the cell type. From left to right are KD and OVCF samples, the KD sample, and the OVCF sample. (B) Violin plots of marker genes for monocyte lineage cell subsets. (C) Proportional pie chart of monocyte lineage cell subsets in respective samples. (D) The left panel shows the analysis of the single-cell trajectories of all monocytes, macrophages, and osteoclasts in KD and OVCF samples; the right panel shows the single-cell trajectories of the above three cells in divided samples. (E) The left panel is a heatmap of branch point 1 shown in Figure 3D, the right panel is based on the GO terms of the three modules, and the right side of the heatmap is the genes with differentially high expression of the GO terms corresponding to each of the three modules.

FIGURE 4

Endothelial cells actively participate in inflammatory infiltration and coupling with osteopoiesis. (A) All endothelial cells initially grouped were pooled to reconstruct the UMAP profile, with coloration labeled by the cell type. From left to right are the KD and OVCF samples, the KD sample, and the OVCF sample. (B) Bubble plots of marker genes for the endothelial cell subcluster. (C) Bar chart of the endothelial cell subcluster in respective sample proportions. (D) Biological process of GO enrichment was performed for genes differentially upregulated in endothelial cells in KD samples compared with OVCF samples. The left panel shows arterial endothelial cells; the middle panel shows capillary endothelial cells; and the right panel shows venous endothelial cells. (E) SCENIC was used to analyze the heatmap of the number and intensity of genes regulated by each transcription factor in each of the three endothelial cell subsets in the two samples, with red representing stronger regulation and blue representing weaker regulation.

FIGURE 5

T and B cells actively participate in regulating bone homeostasis. (A) All T cells initially grouped were combined to reconstruct the UMAP profile, with coloration labeled by the cell type. From left to right are the KD and OVCF samples, the KD sample, and the OVCF sample. (B) Bubble plots of marker genes for the T-cell subcluster, with larger dots representing stronger expression and dots colored to distinguish KD and OVCF samples. (C) Bar chart of the T-cell subcluster in respective sample proportions. (D) Enriched signaling pathways were analyzed using QuSAGE analysis for two T-cell subclusters in two samples, with horizontal axis labels representing cell types in different samples, vertical axis labels representing signaling pathway gene sets, square colors representing statistical probability levels, and darker colors indicating greater significance. The orange bubble size indicates relative activation, and blue indicates relative inactivity. (E) All B cells initially grouped were combined to reconstruct the UMAP profile, with coloration labeled by the cell type. From left to right are KD and OVCF samples, the KD sample, and the OVCF sample. (F) Bubble plots of marker genes for B-cell subclusters, with larger dots representing stronger expression and dots colored to distinguish KD and OVCF samples. (G) Bar chart of B-cell subclusters in respective sample proportions. (H) Corresponding volcano plots were drawn for the top five differentially highly expressed genes up- and downregulated in the KD samples compared with the OVCF samples, with red representing differentially upregulated genes and blue representing differentially downregulated genes.

FIGURE 6

Intercellular communication reveals the osteoimmune microenvironment in human vertebral bone tissue. (A) CellphoneDB analysis was performed on all cell types in KD samples, with the transverse longitudinal axis as the ligand cell group and the longitudinal axis as the receptor cell group, and the colored square represents the number of protein interaction relationships between cell groups, with red representing the more significant communication and blue representing the weaker interaction. (B) Network relationship of osteoclasts with other cells when they act as a ligand. (C) Osteoclasts act as ligand–receptor pairs with other cells when they are ligands, with the transverse longitudinal axis representing cell groups with interactions, red representing ligands, blue representing receptors, the longitudinal axis representing protein interaction pairs, red representing ligands, blue representing receptors, dots representing significant sizes, larger dots representing more significance, and dot color representing the magnitude of interaction intensity, the greater the intensity, the more red, and the smaller the intensity, the more blue. (D) Network relationship of osteoclasts with other cells when they act as receptors. (E) Osteoclasts act as ligand–receptor pairs with other cells when they are receptors, with the transverse longitudinal axis representing cell groups with interaction relationships, red representing ligands, blue representing receptors, the longitudinal axis representing protein interaction pairs, red representing ligands, blue representing receptors, dots representing significant sizes, larger dots representing significant sizes, dot colors representing interaction intensity sizes, greater intensities being more biased toward red, and smaller intensities being more biased toward blue.