Single-cell RNA landscape of the osteoimmunology microenvironment in periodontitis

Abstract

Single-cell RNA sequencing (scRNA-seq) enables specific profiling of cell populations at single-cell resolution. The osteoimmunology microenvironment in the occurrence and development of periodontitis remains poorly understood at the single-cell level. In this study, we used single-cell transcriptomics to comprehensively reveal the complexities of the molecular components and differences with counterparts residing in periodontal tissues. Methods: We performed scRNA-seq to identify 51248 single cells from healthy controls (n=4), patients with severe chronic periodontitis (n=5), and patients with severe chronic periodontitis after initial periodontal therapy within 1 month (n=3). Uniform manifold approximation and projection (UMAP) were further conducted to explore the cellular composition of periodontal tissues. Pseudotime cell trajectory and RNA velocity analysis, combined with gene enrichment analysis were used to reveal the molecular pathways underlying cell fate decisions. CellPhoneDB were performed to identify ligand-receptor pairs among the major cell types in the osteoimmunology microenvironment of periodontal tissues. Results: A cell atlas of the osteoimmunology microenvironment in periodontal tissues was characterized and included ten major cell types, such as fibroblasts, monocytic cells, endothelial cells, and T and B cells. The enrichment of TNFRSF21⁺ fibroblasts with high expression of CXCL1, CXCL2, CXCL5, CXCL6, CXCL13, and IL24 was detected in patients with periodontitis compared to healthy individuals. The fractions of CD55⁺ mesenchymal stem cells (MSCs), APOE⁺ pre-osteoblasts (pre-OBs), and IBSP⁺ osteoblasts decreased significantly in response to initial periodontal therapy. In addition, CXCL12⁺ MSC-like pericytes could convert their identity into a pre-OB state during inflammatory responses even after initial periodontal therapy confirmed by single-cell trajectory. Moreover, we portrayed the distinct subtypes of monocytic cells and abundant endothelial cells significantly involved in the immune response. The heterogeneity of T and B cells in periodontal tissues was characterized. Finally, we mapped osteoblast/osteoclast differentiation mediators to their source cell populations by identifying ligand-receptor pairs and highlighted the effects of Ephrin-Eph signaling on bone regeneration after initial periodontal therapy. Conclusions: Our analyses uncovered striking spatiotemporal dynamics in gene expression, population composition, and cell-cell interactions during periodontitis progression. These findings provide insights into the cellular and molecular underpinning of periodontal bone regeneration.

Keywords: Mesenchymal stem cells; Osteoimmunology; Periodontitis; Single-cell RNA-seq; alveolar bone.

Figure 1

Overview of the 51248 single cells from periodontal tissues of HCs, PDs, and PDTs. A. Study overview. B. Uniform Manifold Approximation and Projection (UMAP) of the 51248 cells, colored by cell-type annotation from left to right: the total corresponding donors and the different conditions (HC, PD, and PDT). MDSCs: Myeloid-derived suppressor cells. C. UMAPs as in (B) but colored by expression of key cell-type marker genes. D. The box plots showing the percentage of cells for each of ten clusters as in (C) from HC (blue, n=14552), PD (red, n=19865), and PDT (green, n=16831) samples with plot center, box, whiskers and point corresponding to the median, IQR, 1.5 × IQR, and >1.5 × IQR respectively. Statistical analysis was performed using unpaired two-tailed t-tests.

Figure 2

Distinct subclusters of the osteoblastic lineage in fibroblast cells cluster. A. Left panels: Uniform Manifold Approximation and Projection (UMAP) of 2194 fibroblasts (as in Figure 1A), annotated and colored by the sample type of origin (HCs, PDs, and PDTs) and clustering. Right panels: UMAPs color-coded for expression (gray to red) of key cell-type markers to define the clusters. Red contours (from left to right): Fibroblasts (cluster 1/2/3/5/8/9), Pericytes (cluster 4/7), Myofibroblasts (cluster 6), and Proliferative cell (cluster 10). B. Left panel: UMAP of 1358 fibroblasts (cluster 1/2/3/5/8/9 in A), annotated and colored by clustering. Center panels: violin plots showing distinct expressions of the selected marker genes (Row) in each subcluster. Right panel: identified subpopulations of Fibroblasts (subcluster 1-7) with the percentages shown. PLF: Periodontal Ligament Fibroblasts; MSC: Mesenchymal Stem Cell; pre-OB: pre-Osteoblast; OB: Osteoblast; FB: Fibroblast. C. The box plots showing the percentage of cells for each of seven subclusters as in (B) from HC (blue, n=655), PD (red, n=1095), and PDT samples (green, n=158) with plot center, box, whiskers, and points corresponding to the median, IQR, 1.5 × IQR and >1.5× IQR, respectively. Statistical analysis was performed using unpaired two-tailed t-tests. D. Representative images of in vitro assays to determine multiple differentiation potential of CD55⁺NT5E⁺LepR⁺ cells. Periodontal cells of a clinically healthy donor were collected and stained with various antibodies and subjected to flow cytometry sorting. CD55⁺NT5E⁺LepR⁺ cells were collected and subjected to the following assays. For osteogenic induction, colony-forming units (CFU-F, left)/alkaline phosphatase (ALP)-positive colony-forming units (CFU-ALP, middle)/nodules (right) formation were examined. For adipogenic induction, Oil red staining was performed to examine adipocytes. For chondrogenic induction, Alcian blue staining was performed to examine chondrocytes. E. Representative images of immunofluorescence staining (left panels) and the quantification (right panel) of periodontal tissues from HC (left) and PD samples (right) for double fluorescent analysis of CD55 (green) and LepR (red) expression. Top panel: Merged; bottom left: CD55; bottom right: LepR. Scale bar = 100μm. F. Upper panel: RNA velocity analysis of MSCs (subcluster 2), pre-OBs (subcluster 3), and OBs (subcluster 4) with velocity field projected onto the UMAP plot of fibroblast subclusters from Figure 2B. Arrows show the local average velocity evaluated on a regular grid and indicate the extrapolated future states of cells. Bottom panel: Monocle pseudotime analysis revealing the progression of MSCs (subcluster 2), pre-OBs (subcluster 3), and OBs (subcluster 4). Trajectory reconstruction of all single cells revealing three branches: Pre-branch, Fate 1, and Fate 2. G. Heatmap showing the scaled expression of differently expressed genes in three branches as in (F), cataloged into three major gene clusters (labels on left) that vary as a function of pseudotime, highlighting specific representative genes in each gene cluster along the right margin. From the center to the left of the heatmap, the kinetic curve from the pre-branch along the trajectory to fate 2 branch. From the center to the right, the curve from pre-branch to successful branch. H. GO analysis of differently expressed genes associated with three gene clusters as in (G) identified unique response pathways for each branch.

Figure 3

The multilineage differentiation capacity of pericytes. A. Left panel: UMAP of 521 cells as pericytes (cluster 4/7) and myofibroblasts (cluster 6) in Figure 2A, annotated and colored by clustering. Center panels: violin plots showing distinct expressions of the selected marker genes (Row) in each subcluster. Right panel: identified subpopulations of pericytes (subcluster 1-6) with the percentages shown. B. The box plots showing the percentage of cells for each of six subclusters as in (A) from HC (blue, n=42), PD (red, n=373), and PDT (green, n=106) samples with plot center, box, whiskers, and points corresponding to the median, IQR, 1.5 × IQR and >1.5× IQR, respectively. Statistical analysis was performed using unpaired two-tailed t-tests. C. Upper panel: Monocle pseudotime analysis revealing the progression of six subclusters as in (A). Trajectory reconstruction of all single cells revealing three branches: Pre-branch, Fate 1, and Fate 2. Bottom panel: RNA-Velocity analysis of the pseudotime trajectory of pericyte and myofibroblast. Direction of the arrows points to the cells are heading toward; length of the arrows reflects how fast the cells are heading toward a particular fate. D. The expression dynamics of selected marker genes go from CXCL12⁺ MSC-like pericyte and bifurcating into two branches with respect to pseudotime coordinates. All single cells in the six subclusters are colored based on (A) and ordered based on pseudotime. E. Heatmap showing the scaled expression of differently expressed genes in three branches as in (C), cataloged into three major gene clusters (labels on left) that vary as a function of pseudotime, highlighting specific representative genes in each gene cluster along the right margin. From the center to the left of the heatmap, the kinetic curve from the root along the trajectory to Osteoblastogenesis. From the center to the right, the curve from the root to Myofibrogenesis. F. GO analysis of differently expressed genes associated with three gene clusters as in (E) identified unique response pathways for each branch.

Figure 4

Trajectory of the osteoclast maturation in periodontal tissues. A. UMAP of 2874 monocytic clusters, annotated and colored by the sample type of origin (HCs, PDs, and PDTs) and clustering. Red contours: Macrophages (cluster 3, 4, 5; left); Monocytes (cluster 1, 2; right). B. UMAPs color-coded for expression (gray to red) of key cell-type markers to define the nine cell groups: Monocytes (cluster 1/2), Macrophages (cluster 3/4/5), OCs (cluster 6), Proliferative cells (cluster 7), CD1C⁺ DCs (cluster 8), LAMP3⁺ DCs (cluster 9), CD5⁺ DCs (cluster 10), CLEC9A⁺ DCs (cluster 11) and Plasmacytoid DCs (cluster 12). OCs: osteoclasts; DCs: Dendritic cells. C. The box plots showing the percentage of cells for each of nine groups as in (B) from HC (blue, n=627), PD (red, n=1577), and PDT samples (green, n=670) with plot center, box, whiskers, and points corresponding to the median, IQR, 1.5 × IQR and >1.5× IQR, respectively. Statistical analysis was performed using unpaired two-tailed t-tests. D. Monocle pseudotime analysis revealing the progression of three major groups: monocyte (cluster 1, 2), macrophage cell (cluster 3, 4, 5), and osteoclast (cluster 6). Trajectory reconstruction of all single cells revealing a continuous lineage path without a branch colored according to clusters as in (A) (top). The distribution of single cell from each cluster mapped in a continuous lineage path, respectively (bottom). E. Heatmap showing the scaled expression of differentially expressed genes in the cells of each major group as in (D), cataloged into three major gene clusters (labels on right) that vary as a function of pseudotime, highlighting representative differentially expressed genes in each gene cluster along the right margin. From the left to the right of the heatmap, a continuous lineage path from monocytes to OCs. OCs: osteoclasts; DEGs: differentially expressed genes. F. Pathway enrichment analysis associated with three gene clusters as in (E). G. Heatmap showing the scaled expression of differentially expressed transcription factor genes along with the pseudotime curve as in (D), cataloged into three major gene clusters (labels on right) that vary as a function of pseudotime, highlighting the top 50 transcription factors along the right margin. From the left to the right of the heatmap, a continuous lineage path from monocytes to OCs. OCs: osteoclasts; TFs: transcription factors. H. The expression dynamics of representative genes that encoded secretion protein. These genes are differentially expressed across pseudotime corresponding to three gene clusters as in (E). All single cells in the six subclusters are colored based on (A) and ordered based on pseudotime.

Figure 5

Abundant endothelial cells involved in immune response. A. UMAP of 4461 endothelial clusters, annotated and colored by the sample type of origin (HCs, PDs, and PDTs) and clustering. Red contour: Venous ECs (cluster 1/2/3/5/6/7/8). EC: Endothelial Cell. B. UMAPs color-coded for expression (gray to red) of key cell-type markers to define the four cell groups: Venous ECs (cluster 1/2/3/5/6/7/8), Arterial ECs (cluster 4), Lymphatic ECs (cluster 10), and Proliferative ECs (cluster 9). EC: Endothelial Cell. C. The box plots showing the percentage of cells for each of four cell groups as in (B) from HC (blue, n=153), PD (red, n=2999), and PDT samples (green, n=1309) with plot center, box, whiskers, and points corresponding to the median, IQR, 1.5 × IQR and >1.5× IQR, respectively. Statistical analysis was performed using unpaired two-tailed t-tests. D. Heatmap of the mean value of area under the recovery curve (Aucell) value of expression regulation by the transcription factors, as estimated using SCENIC, for the indicated three ECs clusters (Row). Shown the top 25 transcription factors genes (TF) (labels on right) having the highest difference in expression regulation estimates among HCs, PDs, and PDTs (Row). E. UMAP of endothelial cells, dots coded the for the Aucell value of the estimated regulon activity of selected TFs and the corresponding targeted genes in Arterial ECs (left), Lymphatic ECs (center), Venous ECs (right). EC: Endothelial Cell. F. Volcano plots displaying the differentially expressed genes in Venous ECs among PD VS. HC, PDT VS. PD, and PDT VS. HC (from left to right). Each dot represented one gene. Representative differentially expressed genes (blue) are indicated. Green dots, differentially down-regulated genes with log_eFC < -0.25 and FDR < 0.05; red dots, differentially up-regulated genes with log_eFC > 0.25 and FDR < 0.05; gray dots, non-differentially expressed genes.

Figure 6

The heterogeneity of T & B cells in periodontal microenvironment. A. UMAP of 12805 T cell clusters, annotated and colored by the sample type of origin (HCs, PDs, and PDTs) and clustering. Red contours in right panel: CD4⁺ T cells (cluster 1/2/4/8/12) (left), CD8⁺ T cells (cluster 3/5/6/7) (center), NK T cells (cluster 9/10) (right). B. UMAPs color-coded for expression (gray to red) of key cell-type markers to define the four cell groups: CD4⁺ T cells (cluster 1/2/4/8/12), CD8⁺ T cells (cluster 3/5/6/7), NK T cells (cluster 9/10), Proliferative T cells (cluster 11). C. Left panel: UMAP of 5066 CD4⁺ T cells, annotated and colored by clustering with the label of each subcluster. Right panel: identified six subclusters of CD4⁺ T cells with the percentages shown. D. Left panel: UMAP of 5951 CD8⁺ T cells, annotated and colored by clustering with the label of each subcluster. Right panel: identified seven subclusters of CD8⁺ T cells with the percentages shown. E. Heatmap of the QuSAGE activity for the enrichment pathways (labels on the right) that associated with cell differentially expressed genes among seven subclusters of CD8⁺ T cells from HC, PD, and PDT samples (labels on the top). Red indicates an increased average expression of genes in the modules. F. Bubble plot showing expressions (dots) of the selected immune checkpoint inhibitor genes (columns) in each subcluster of CD8⁺ T cells (rows, as in B). Dot colored by the average expression level, and dot size proportional to the percentage expression. G. UMAP of 10576 B & Plasma cell clusters, annotated and colored by the sample type of origin (HCs, PDs, and PDTs) and clustering. Red contours: B Cell (cluster 1/6/9) (left), Plasma Cell (cluster 2/3/4/5/7/10) (right). H. UMAPs color-coded for expression (gray to red) of key cell-type markers to define the three cell groups: B Cell (cluster 1/6/9), Plasma Cell (cluster 2/3/4/5/7/10), and Proliferative cell (cluster 8). I. Top left panel: UMAP of 3248 B cells, annotated and colored by clustering. Red contours: Naive B cells (left), Memory B cells (center), ABCs (right). Top right panel: identified subclusters of B cells with the percentages shown. Bottom panel: Violin plots showing distinct expressions of the selected marker genes (row) in each subcluster (labels on left) of B cells. ABCs: Activated B Cells. J. Pathway enrichment analysis associated with genes upregulated (top) between PD vs. HC and downregulated (bottom) genes between PDT vs. PD.

Figure 7

Cell-cell communication among all cell types in the periodontal osteoimmunology microenvironment. A. The circo plot showing the potential cell interactions among ten major cell types as in Figure 1B in periodontal tissues predicted by CellphoneDB. The node size represents the number of interactions. The width of the edge represents the number of significant ligand-receptor pairs between the two cell types. B. Venn diagram representing the interaction between the significant ligand-receptor pairs identified by CellPhone-DB analysis on scRNA-seq data from HC, PD, and PDT samples. C. The dot plot generated by CellPhoneDB showing potential ligand-receptor pairs associated with osteoblastogenesis between Pre-OB and all detected cellular types in PD group. Dots colored by mean expression of ligand-receptor pair between two clusters and dots size proportional to the value of -log10 (P Value). D. Stacked bar graph showing the number of all significantly expressed pairs between ECs and Pre-OB on different conditions (HC: blue bar; PD: red bar; PDT: green bar) and the Ephrin-EPH receptor pairs within these conditions (orange bar). E. The dot plot generated by CellPhoneDB showing potential ligand-receptor pairs associated with Ephrin-Eph receptor signaling pathway between Pre-OB and all detected cellular types in PDT group. Dots colored by mean expression of ligand-receptor pair between two clusters and dots size proportional to the value of -log10 (P Value). F. Representative images of paraffin sections from periodontal tissues of HC (left panel) and PD samples (right panel) were subjected to IF with anti-ALP (pre-OBs, green), anti-Eph A7 (red), anti-CD31 (ECs, blue), and anti-Ephrin A1 (white) Abs. Adjacent paraffin sections were stained with H & E staining to show the location. Top left: merged; top right: ALP; middle right: Eph A7; bottom left: CD31; bottom middle: Ephrin A1; bottom right: H & E. Scale bar = 50 µm. The quantifications of Eph A7⁺ALP⁺ Pre-OBs and Ephrin A1⁺CD31⁺ ECs per field were shown. Pre-OB: Pre-osteoblast, EC: Endothelial cell.

Figure 8

Predicted regulatory network centered on the osteoimmunology microenvironment of periodontal tissue after initial periodontal therapy. Sc-RNA seq analysis of the osteoimmunology microenvironment of periodontal tissue after initial periodontal therapy. In this Figure, the Ephrin-Eph signaling and the interactions regarding OB/OC differentiation are shown at the level of intercellular interactions. Pre-OB: Pre-osteoblast; OB: Osteoblast; EC: Endothelial cell; OC: Osteoclast; OCP: Osteoclast precursor cell.