Single-cell transcriptome analysis of stem cells from human exfoliated deciduous teeth investigating functional heterogeneity in immunomodulation

Identification of SHED

SHEDs form spindle-like morphology (Fig. 1A) and maintain proliferation activity and morphology from the third to sixth generations. Flow cytometric analysis demonstrates that cells express CD105, CD73, and CD90 (Fig. 1B–D), which are mesenchymal stem cell markers. Immuno-related markers were also analyzed, and those cells showed a negative signal of CD34 and CD45 (Fig. 1E–F), the traditional surface proteins known as hematopoietic stem cells and leukocytes. Those flow cytometric results indicate that SHEDs are cultured in a low differentiation state.

Fig. 1
figure 1

Identification of SHED (A) Brightfield image of SHED. (BF), Flow cytometric analysis showing the expression of CD105, CD73, CD90, CD34 and CD45.

SHEDs lie along a continuum of differentiation states

To investigate the cell-to-cell and immunomodulatory heterogeneity in SHED at the single-cell transcriptome level, cells isolated from one human deciduous exfoliated tooth and passaged to the third generation were collected and used for scRNA-seq. An atlas with 18,211 single‐cell transcriptomes was constructed after stringent quality control (Supplementary Fig. 1A–F). The UMAP was applied to visualize the 14 clusters identified by the expression of top variable genes (TVGs) (Supplementary Fig. 1G). The top 10 differentially expressed genes are shown in Supplementary Table 1. In order to understand the temporal relationships within the SHED population, we performed the pseudotime analysis of all SHED cells using the Monocle2 package. Eleven cell subpopulations of SHED were identified, indicating that SHED was a heterogeneous cell population. The curated instructive gene list for pseudotime analysis are provided in Supplementary Table 2. For characterizing the stemness of 11 subpopulations, the two genes, including Thy-1 membrane glycoprotein (THY1, CD90), seen as one of the classical markers for MSCs24, and cell surface glycoprotein MUC18 (MCAM, CD146), which was a classical marker for dental stem cells, were chosen as the markers to explore the differentiation trajectory of SHED. We found that the alterations of THY1 and MCAM expression existed in 11 SHED subpopulations. The highest expression of both two markers was exhibited at SHED in state 7 (S7), representing that S7 was the SHED subpopulation of strong stemness (Fig. 2A–F). Therefore, the possible developmental trajectory of the SHED was further constructed, and the 11 subpopulations were reordered along a pseudotime axis by designating the S7 as the root (Fig. 2G–H). The multiple branches in the trajectory represented the multidirectional differentiation potential of SHED. The SHED in state 1 (S1) was observed in trajectory termini. These results further demonstrated that SHED lies along a continuum of differentiation states. S7 was in the low differentiation state, whereas S1 was in the relatively high differentiation state. To present the results more intuitively, we will refer to S1 as the “late state” and S7 as the “early state” in the following sections.

Fig. 2
figure 2

The expression variations of MCAM and THY1 in SHED differentiation. (A) Scatter plots showing the distribution levels of MCAM in SHED subpopulations at different states. (B) Violin plots displaying the expression levels of MCAM. (C) Scatter plots showing the expression change of MCAM in SHED subpopulations ordered by the relative expression level of MCAM. (D) Scatter plots showing the distribution levels of THY1 in SHED subpopulations at different states. (E) Violin plots displaying the relative expression levels of THY1. (F) Scatter plots showing the expression change of THY1 in SHED subpopulations ordered by the relative expression level of THY1. (G) Reconstructed trajectory plot showing the SHED differentiation direction by designating S7 as the root. The degree of differentiation from low to high was represented by color from the deeper blue to the lighter blue. (H) The principal component graph of cell differentiation trajectory of SHED showing the S1 was trajectory termini. The clusters were colored by identifies.

The heterogeneity of immunomodulatory capacity in SHED

In order to explore the shared and distinct biological function among SHED subpopulations at different differentiation states, the late state and early state, which presented the most significant differences in cellular differentiation, were selected as subjects for subsequent analyses. After identifying upregulated differentially expressed genes (DEGs) between late state and early state, the enrichment analysis using the Reactome database was performed. The top 25 enriched pathways for upregulated DEGs in late state and early state are shown in Fig. 3. To clarify the functions of each population, the late state and early state were further analyzed separately. A significant difference between late state and early state populations was realized in the results of pathways enrichment. Upregulated DEGs in late state and early state had no common pathway enrichment. The enrichment for pathways in late state was related to growth factor signaling, including signaling by Platelet-derived growth factor (PDGF), signaling by bone morphogenetic protein (BMP), signaling by transforming growth factor beta (TGF-β), signaling by nerve growth factor (NGF), integrin cell surface interaction, and signaling by epidermal growth factor receptor (EGFR). However, the enrichment for pathways in early state, which were mainly associated with metabolism, cell proliferation, and maintaining stemness, included metabolism of proteins, gene expression, cell cycle mitotic, cell cycle checkpoints, and signaling by Wnt (Fig. 3A, highlighted in the red font).

Fig. 3
figure 3

Functional characterization between SHED in late state and early state. (A) Reactome enrichment analysis for upregulated differentially expressed genes (DEGs) of late state (state 1, S1) or/and early state (state 7, S7) in different cells number. Dot plot showing the terms with significant. The gene ratio (enriched genes/total number of genes) was represented by the size of dot. The adjusted p-value for enrichment analysis was indicated by the color of dot. Red label highlights the important terms related to SHED functions. Gene set enrichment analysis (GSEA) exhibiting that some upregulated DEGs of S1 (B), and S7 (C), were significantly enriched cytokine-cytokine receptor interaction function and leukocyte transendothelial migration function, respectively.

To further investigate the immunological functional significance of corresponding genes from late state and early state, the GSEA was carried out using MSigDB C5: BP GO biological process collection (Supplementary Table 3, Fig. 3B–C). To be noticed, we detected that the upregulated DEGs in late state are significantly enriched in cytokine-cytokine receptor interaction function (NES = 1.492 q-value = 0.023). Nevertheless, the upregulated DEGs in early state are significantly enriched in leukocyte transendothelial migration function (NES = − 1.599, q-value = 0.035). These results implied that SHED is a highly heterogeneous cell population with different immunoregulatory functions during differentiation. According to the results mentioned above, we speculated that late state may provide a microenvironment suitable for immune cell differentiation, whereas early state may have a solid ability to attract immune cell.

SHED in low differentiation state have a stronger chemotactic ability towards immune cells

To validate the above observation, the immune cell recruitment scores of each MSC population were calculated by single-sample GSEA (ssgsea) function in the GSVA package using an immune cell recruitment set23 (Supplementary Table 4). The immune cell recruitment scores were calculated to quantify the ability of late state and early state to induce the chemotaxis of immune cells (Fig. 4). A higher recruitment score indicated more vital chemotactic ability. We found that the early state had significantly higher recruitment scores in various immune cell types than late state, including B cell, CD4+T cell, CD8+T cell, Macrophage, monocyte, neutrophil, NK cell, Th17, Th22 and Treg (Supplementary Table 5), which support our previous results. The chemotaxis‐inducing effect of early state on various immune cells was more substantial than that of late state.

Fig. 4
figure 4

Comparison of immune cell recruitment scores between SHED in late state and early state. The immune cell recruitment scores with significant difference between late state (state 1, S1) and early state (state 7, S7) are shown as box plots with medians (lines inside boxes). Plots represent outliers. Wilcox test was used for p-value. The x-axis shows the different states, and the y-axis shows the normalized ssGSEA score.

Single-cell transcriptome and chemotactic ability comparison of MSCs across multiple tissue origins

To gain insight into the heterogeneity and immune cells chemotactic ability of MSCs derived from different tissues, we next integrated and compared our single-cell SHED transcriptome data with two previous scRNA-seq studies of human bone marrow mesenchymal stem cells (hBMSCs) and hUCMSCs20,21. The clustering results showed that the five samples were well integrated, indicating that batch effects between samples were removed (Fig. 5A). As expected, MSC clusters presented a tissue-type-dependent distribution in UMAP, while the projections overlapped in MSCs from the same tissue origin (Fig. 5B,C). It revealed that significant transcriptomic heterogeneity existed in MSCs from different tissues. Cell composition of the three datasets was visualized, and 16 clusters were identified (Fig. 5D). We further assessed the differently expressed genes in SHED, hBMSCs, and hUCMSCs (Fig. 5E–J). Among these, the expression of stemness genes THY1, MCAM, MET, and growth factors genes VEGFA, VEGFD, and TGF-β1 showed interesting patterns. SHED exhibited a lower level of THY1 and a higher level of MCAM than hBMSCs and hUCMSCs. In contrast, MET is barely expressed in SHED.

Fig. 5
figure 5

Characterization of multiple‐tissue derived MSCs populations. (A) UMAP visualization of hBMSCs, SHED and hUCMSCs, respectively. (B) UMAP visualization of the Harmony integrated BMSCs, SHED and hUCMSCs, colored by sample source. (C) Distributions of the five samples visualized in UMAP. Populations were colored by tissue type. (D) Visualization of 17 color-coded clustering of hBMSCs, SHED and hUCMSCs (n = 24,887 cells) using the UMAP. (EJ) The distribution of selected stem cell surface makers genes (MCAM, THY1, and MET) and growth factors gene (VEGF-A, VEGF-D, TGF-β) expression in hBMSCs, SHED and hUCMSCs. The cells were respectively marked with gray and purple color based on the low and high expression of selected genes.

Moreover, SHED had increased expressions in all growth factor genes. The differential expression of these important genes may imply the different functional potentials of MSC populations. Following function enrichment analysis showed that the regulated DEGs of SHED were mainly enriched in hemostasis, signaling by PDGF, and integrin cell surface interaction when compared with hBMSCs or hUCMSCs. The three mentioned functions were not enriched in hBMSCs or hUCMSCs (Fig. 6A,B). Next, we calculated the immune cell recruitment scores of each MSC population. All scores are compared (Fig. 6C) and visualized in a radar chart (Fig. 6D). We found significant differences in immune cell recruitment scores, including CD4+T cell, CD8+T cell, macrophage, neutrophils, NK cell, Th17, Th22, and Treg. These observations suggested the chemotaxis ability heterogeneity in MSCs from different tissues to recruit immune cells. As shown in the radar graph, SHED and hBMSCs played wider and stronger chemotaxis to various immune cells than hUCMSCs. Despite this, SHED had more attractive roles on Th17, Th22, and Treg, whereas hBMSCs had on macrophage. hUCMSCs only presented the enhanced recruitment ability of macrophages compared with SHED and hBMSCs.

Fig. 6
figure 6

Functional characterization among hBMSCs, SHED and hUCMSCs. (A) Reactome enrichment analysis for upregulated DEGs of hBMSCs and SHED in different cells number. (B) Reactome enrichment analysis for upregulated DEGs of hUCMSCs and SHED in different cells number. Dot plot shows the terms with significant. The gene ratio (enriched genes/total number of genes) was represented by the size of dot. The adjusted p-value for enrichment analysis was indicated by the color of dot. Red label highlights the important terms related to SHED functions. (C) The immune cell recruitment scores with significant difference among hBMSCs, SHED and hUCMSCs are shown as box plots with 25th and 75th quartiles (limits inside boxes). The y-axis shows the x-axis shows the normalized ssGSEA score, and the different immune cell types. Bars with different colors correspond to different groups. Kruskal–Wallis rank test was used for p-value. *** p < 0.001. (D) The percentages calculated by the median of each significant immune cell recruitment score in the corresponding maximum median of score among hBMSCs (red line), SHED (blue line) and hUCMSCs (green line) visualized as a radar plot, and the values range from 1 to 100%.

Immunoregulatory characterization of SHED

We put the SHEDs in a low differentiation state co-cultured with activated or inactivated PBMC to validate the strong chemotactic ability and another effect on immune cells. The result shows a tendency to enhance the proliferation of Treg (CD4+CD25+CD127) (Fig. 7A) under an appropriate concentration of IL-2. Although there are different original ratios and promote efficiency among the diverse PBMCs, a clear increase of Treg was observed in inactivated PBMC between control and SHEDs treated. (Fig. 7B) we also co-cultured activated PBMC and SHEDs to mimic an inflammatory environment; the SHEDs can down-regulate inflammatory factors such as TNF-α (Fig. 7C).

Fig. 7
figure 7

Immunoregulatory characterization of SHED (A) Flow cytometric profiles of CD25 and CD127 in CD4+ cells in the PBMC of the control and SHED-treated groups. (B) Numerical values denote the percentage of CD4+CD25+CD127 in each group. (C) Numerical values denote the concentration of TNF-α in each group.