Identification of FSMSCs and HuMSCs
According to the minimum identification standard about MSCs of the ISCT, the cells isolated from foreskin tissue and umbilical cord were confirmed as MSCs, respectively. First, FSMSCs and HuMSCs were able to attach plastic surfaces, and they exhibited a typical long, spindle-shaped fibroblast-like morphology under the microscope. Second, positive surface markers, such as CD73, CD90, and CD105, were expressed in more than 95% of the cells, while negative surface markers, such as CD45, CD34, CD11b, CD19, and HLA-DR, were expressed in fewer than 2% of the FSMSCs and HuMSCs (Fig. 1a, b). Red spherical vacuoles were observed by Oil Red O staining after induction of adipogenic differentiation in FSMSCs or HuMSCs (Fig. 1c). This finding indicated that lipid droplets formed in these cells. After the induction of osteogenic differentiation, scattered round crimson nodules were observed in FSMSCs or HuMSCs by 2% ARS staining (Fig. 1c). Massive amounts of calcium were deposited in the cytoplasm. Finally, blue acid mucopolysaccharides were visible in FSMSCs or HuMSCs by Alcian Blue staining after the induction of chondrocyte differentiation (Fig. 1c).
Immunomodulatory abilities of FSMSCs and HuMSCs
The secretion of IL-1β, IL-6, IL-10, TNF-α, and TGF-β1 was measured to compare the immunomodulatory abilities of FSMSCs and HuMSCs through ELISA. First, we analyzed the differences in autocrine inflammation-related cytokines among FSMSCs, HuMSCs, and PBMCs (Fig. 2a, Additional file 1: Tables S1–S4). The differences in the basal secretion of IL-1β, IL-10, and TGF-β1 in these three cell types were not statistically significant at any time point. The basal secretion of IL-6 was significantly higher in HuMSCs than in FSMSCs and PBMCs, while there was no significant difference between FSMSCs and PBMCs. The basal secretion of TNF-α was significantly higher in PBMCs than in FSMSCs and HuMSCs, while there was no significant difference between FSMSCs and HuMSCs. In addition, we analyzed whether there were differences in the basal secretion of inflammation-related cytokines among these three cell types at different time points (Fig. 2a, Additional file 1: Table S1–S4). The basal secretion level of IL-6 increased with time, while the basal secretion level of TGF-β1 showed a fluctuating pattern. The differences in the secretion of IL-1β, IL-10, and TNF-α were not statistically significant at different time points.
Next, we analyzed the differences in the secretion of inflammation-associated cytokines among FSMSCs, HuMSCs, and PBMCs under LPS stimulation. The differences in the secretion of IL-10 and TGF-β1 were not statistically significant. The secretion of IL-1β and TNF-α was highest in PBMCs, while the secretion level of IL-6 was highest in FSMSCs. In addition, we analyzed whether there were differences in the secretion of inflammation-associated cytokines among MSCs at different time points after stimulation with LPS (Fig. 2b, Additional file 2: Table S5–S8). The secretion of IL-6 increased with time, while the differences in the secretion of IL-10 were not statistically significant. The secretion of IL-1β and TNF-α in PBMCs increased with time. The secretion of TGF-β1 in FSMSCs increased with time.
Moreover, we analyzed the differences in the secretion of inflammation-related cytokines among FSMSCs and HuMSCs cocultured with PBMCs and PBMCs cultured alone (Fig. 2c, Additional file 3: Table S9–S12). Overall, the secretion levels of IL-6, TNF-α, IL-10, and TGF-β1 were higher in FSMSCs cocultured with PBMCs than in HuMSCs cocultured with PBMCs or in PBMCs cultured alone. There was no significant difference in the secretion of IL-1β among these three conditions. In addition, we analyzed whether there were differences in the secretion of inflammation-related cytokines at different time points among these three conditions (Fig. 2c, Additional file 3: Table S9–S12). The secretion of IL-6 increased with time. The secretion levels of IL-10 and TNF-α in FSMSCs cocultured with PBMCs increased with time.
Finally, we mimicked the inflammatory response by stimulating PBMCs with LPS in vitro. FSMSCs and HuMSCs were cocultured with PBMC prestimulation. Then, we compared the immunoregulatory abilities by detecting the secretion of inflammation-related cytokines (Fig. 2d, Additional file 4: Table S13–S16). The secretion of the proinflammatory cytokines IL-1β and TNF-α was lower in the FSMSCs and HuMSCs than in the control cells. The secretion of the bidirectional immunomodulatory cytokine IL-6 and the anti-inflammatory cytokines IL-10 and TGF-β1 was higher in FSMSCs than in HuMSCs. In addition, we investigated whether there were significant differences in the secretion levels of inflammation-associated cytokines in the FSMSCs and HuMSCs at different time points (Fig. 2d, Additional file 4: Table S13–S16). IL-6, IL-10, and TGF-β1 secretion increased with time in FSMSCs cocultured with prestimulated PBMCs. IL-10 secretion increased with time in HuMSCs cocultured with prestimulated PBMCs. In summary, FSMSCs and HuMSCs both inhibited the secretion of the proinflammatory cytokines IL-1β and TNF-α from PBMCs. FSMSCs were more significantly able to upregulate the secretion of the bidirectional immunomodulatory cytokine IL-6 and the anti-inflammatory cytokines IL-10 and TGF-β1 than HuMSCs.
Quality control, clustering and biological annotation
Quality control, clustering, and biological annotation were performed for FSMSCs and HuMSCs. First, the Seurat objects for FSMSCs and HuMSCs were created separately. Second, quality control was performed, and low-quality cells with fewer than 200 genes, more than 20,000 genes, more than 10% mitochondrial genes and more than 10% erythrocyte genes were filtered. Then, we obtained 7335 FSMSCs, in which the median UMI was 33202 and the median gene number was 5440. Moreover, we obtained 12542 HuMSCs, in which the median UMI was 15191 and the median gene number was 3831. This result indicated that the quality of the scRNA-seq data was high. Afterward, normalization was run, and the effect of the mitochondrial genes on scaling was removed. Next, we executed dimensionality reduction by PCA through the first 2000 genes with high heterogeneity and evaluated the first 50 PCs. We evaluated the first 50 PCs with the JackStraw and ElbowPlot functions of the R package Seurat (Additional file 11: Plot S1, Additional file 12: Plot S2). We also calculated the percentages of each principal component in the population variance (Additional file 10). We found that the first 20 PCs with standard deviations could explain most of the variability in the data (Additional file 11: Plot S1b, Additional file 12: Plot S2b). Thus, the first 20 PCs with the most significant p values (Additional file 11: Plot S1a, Additional file 12: Plot S2a) were used for UMAP. The first 20 PCs were also used to cluster cells with a resolution of 0.5, and we ultimately obtained seven clusters of FSMSCs and 10 clusters of HuMSCs (Fig. 3a).
Then, the clusters were annotated according to the surface markers of MSCs [1]. Among FSMSCs, those in C0, C1, C2, C3, and C5 highly expressed the positive marker genes CD73, CD90 and CD105 and lacked expression of the negative marker genes HLA-DRA, CD11b, CD19, CD34, and CD45 (Fig. 3a). This finding indicated that the cell in these clusters were possible MSCs. However, C4 cells highly expressed HLA-DRA, and C6 cells lacked expression of CD73 and CD105. This result indicated that C4 and C6 cells did not satisfy the identification criteria for MSCs. Similarly, among HuMSCs, the cells in C0, C1, C2, C3, C5, and C6 highly expressed the positive marker genes CD73, CD90 and CD105 and lacked expression of the negative marker genes HLA-DRA, CD11b, CD19, CD34, and CD45, suggesting that the cells in these clusters might be MSCs. Nevertheless, C4, C7, C8, and C9 cells lacked expression of CD73 and CD105, indicating that the cells in these clusters did not satisfy the identification criteria for MSCs (Fig. 3a). In summary, the cells in C0, C1, C2, C3, and C5 were MSCs among FSMSCs and accounted for 93.05% of the FSMSCs. The cells in C0, C1, C2, C3, C5, and C6 were MSCs among HuMSCs and accounted for 95.85% of the HuMSCs.
Integration analysis
We integrated the cell subsets annotated as MSCs among FSMSCs and HuMSCs to reveal the potential heterogeneity of MSCs (Fig. 3a). Then, the batch effect was removed by canonical correlation analysis and evaluated with a PCA plot. The FSMSCs and HuMSCs were able to overlap well after integration (Fig. 3d), indicating that the batch effect was corrected well. Afterward, the integrated data were normalized and scaled again. We also executed dimensionality reduction by PCA through the first 2000 genes with high heterogeneity and evaluated the first 50 PCs. We evaluated the first 50 PCs with the JackStraw and ElbowPlot functions of the R package Seurat (Additional file 13: Plot S3). We also calculated the percentages of each principal component in the population variance (Additional file 10). We found that the first 20 PCs with standard deviations could explain most of the variability in the data (Additional file 13: Plot S3b). Thus, the first 20 PCs with the most significant p values (Additional file 13: Plot S3a) were used for TSNE. The first 20 PCs were also used to cluster cells with a resolution of 0.5, through which seven major clusters were obtained (Fig. 3e). Furthermore, we calculated the proportions of these clusters between FSMSCs and HuMSCs, as shown on the bar plot (Fig. 3f). Moreover, the cell cycle phase of each cluster was inferred with the CellCycleScoring function of the R package Seurat (version 4.0.0) [25]. There were elevated proportions of G1 phase cells in C0, C2, and C6, suggesting a resting state. However, G2/M or S phase cells occupied large proportions of C2, C3, C4, and C5 cells (Fig. 3f). In addition, MKI67, a gene that encodes a nuclear protein associated with cell proliferation [27], was highly expressed in C2, C3, C4, and C5 cells (Fig. 3b). This result suggested that the proliferative activity of the cells in C2, C3, C4, and C5 was more substantial than that of the cells in the other subsets. Intriguingly, this result is consistent with the results of previous studies on MSCs cultured in vitro. Three morphologies of MSCs cultured in vitro have been reported: spindle-shaped, large and flattened, and small and prototypical. The first two kinds of cells are more extensive and proliferate more slowly. The latter cells, known as rapid self-renewing cells, [28, 29]. The results indicated that some cell subsets were more active in proliferation than others. Therefore, the cells in C2, C3, C4, and C5 were named proliferative MSCs (Fig. 3e).
Second, C6 cells highly expressed the positive marker genes CD146, PDGFRB and NG2 and lacked expression of the negative marker genes CD31, CD45 and CD34 (Fig. 3b). These genes have been demonstrated to be markers of pericytes, but pericytes can also express the markers of MSCs [30]. Therefore, C6 cells were regarded as possible pericytes (Fig. 3e).
In addition, CXCL12 and PTGIS were highly expressed in C0 cells, while MKI67 and CD146 lacked expression in C0 cells (Fig. 3b). Among these, CXCL12 plays an essential role in the immunoregulatory function of MSCs [48]. PTGIS, prostacyclin I2 synthase, can catalyze the conversion of prostaglandin H2 to prostaglandin I2, which plays an essential role in immunoregulatory function [31]. We also observed that C1 and C6 cells expressed CXCL12 and PTGIS, but the expression levels in these cells were much lower than those in C0 cells. This result suggested that C0, C1, and C6 cells all possessed immunomodulatory abilities. However, the C0 cells exhibited the strongest immunomodulatory function among these cells. C6 cells were identified as pericytes, which are known to exhibit immunomodulatory function [49]. Furthermore, C1 cells expressed low levels of CXCL12, and PTGIS and lacked MKI67 and CD146 expression (Fig. 3b). In addition, a few C1 cells expressed MKI67 (Fig. 3b). Therefore, we constructed a clustering tree based on the similarities of MSC subsets and visualized it with a dendrogram plot. Then, we found that the expression patterns of C1 and C5 cells were closer than those of other cells (Fig. 3c). Moreover, the C5 cells belonged to the proliferative MSCs. The C1 cells were possible precursor cells of C5 cells, while the C1 cells might have represented a transitional stage from C0 cells to proliferative MSCs. Therefore, C1 cells were named progenitor proliferative MSCs, and C0 cells were named immune MSCs (Fig. 3e).
Finally, the proportion of subsets of FSMSCs and HuMSCs were calculated (Fig. 3f). Among them, proliferative MSCs made up 16% of FSMSCs and 68% of HuMSCs, pericytes made up 6% of FSMSCs and 4% of HuMSCs, immune MSCs made up 56% of FSMSCs and 10% of HuMSCs, and progenitor proliferative MSCs made up 22% of FSMSCs and 18% of HuMSCs.
Trajectory inference
The CytoTRACE package was used to infer the differentiation degrees of four MSC subsets (Fig. 4b). A CytoTRACE score nearer to 1.0 indicated a lower degree of differentiation, while a CytoTRACE score nearer to 0.0 indicated a higher degree of differentiation. The degrees of differentiation of proliferative MSCs and progenitor proliferative MSCs were higher than those of immune MSCs and pericytes (Fig. 4b). Some studies have reported that MSCs might originate from pericytes [50]. Therefore, we defined pericytes as the starting point of the differentiation trajectory of MSCs. Next, we constructed two differentiation trajectories based on transcriptional similarity with the R package slingshot (Fig. 4a). In the first trajectory, the cells started from pericytes and differentiated into proliferative MSCs. In the second trajectory, the cells started from pericytes, passed through the immune MSC and progenitor proliferative MSC stages, and ultimately differentiated into proliferative MSCs.
RNA velocity analysis
The transcriptional activity of MSC subsets was inferred by RNA velocity analysis, in which a longer arrow indicates stronger transcriptional activity. We found that most of the arrows for immune MSCs were longer than those for proliferative MSCs (Fig. 4c). This finding indicated that the overall transcriptional activity of immune MSCs was stronger than that of proliferative MSCs.
Gene set enrichment analysis
All cell subsets were subjected to gene set enrichment analysis to assess the potential biological functions of the subsets. First, the enrichment score and adjusted p value of the hallmark gene set of different MSC subsets were calculated (Additional files 6, 7) and then visualized with a heatmap plot (Fig. 4d). The asterisks on the heatmap plot indicate that the p value of the hallmark gene set was less than 0.05.
First, we found 15 statistically significant gene sets with adjusted p values ≤ 0.05 that were filtered from immune MSCs (Fig. 4d). Among them, the immune-related gene sets “IL6 JAK STAT3 SIGNALING”, “INTERFERON GAMMA RESPONSE” and “INTERFERON ALPHA RESPONSE” were upregulated in immune MSCs. The mitosis-related gene sets “MITOTIC SPINDLE”, “G2M CHECKPOINT” and “E2F TARGETS” were downregulated in immune MSCs. Second, 13 and 19 statistically significant gene sets were filtered by a corrected p value ≤ 0.05 from progenitor proliferative MSCs and pericytes, respectively. The immune-related gene set “TNFA SIGNALING VIA NFKB” was upregulated in both progenitor proliferative MSCs and pericytes. Furthermore, in progenitor proliferative MSCs, gene sets related to adipose differentiation, such as “ADIPOGENESIS” and “FATTY ACID METABOLISM”, were preferentially enriched. However, in pericytes, gene sets related to mesenchymal histogenesis, such as “MYOGENESIS”, “ANGIOGENESIS” and “EPITHELIAL MESENCHYMAL TRANSITION”, were preferentially enriched. Finally, 20 statistically significant gene sets with corrected p values ≤ 0.05 were filtered for proliferative MSCs. These cells lacked enrichment of the hypoxia-related gene set “HYPOXIA” and the immune-related gene sets “IL6 JAK STAT3 SIGNALING”, “INTERFERON GAMMA RESPONSE” and “INTERFERON ALPHA RESPONSE”. Nevertheless, the mitosis-related gene sets “MITOTIC SPINDLE”, “G2M CHECKPOINT” and “E2F TARGETS” were significantly upregulated in proliferative MSCs.
Immune-related differential gene expression analysis
First, we compared the gene expression between FSMSCs and HuMSCs and filtered 205 differentially expressed genes with adjusted p values ≤ 0.05 and absolute log2FC values ≥ 2 (Additional file 8). Then, we intersected the immune-related genes with the 205 differentially expressed genes and finally obtained a total of 37 immune-related differentially expressed genes. We divided these 37 intersecting genes into seven immune categories: “Antigen processing and presentation”, “Antimicrobials”, “Cytokine receptors”, “Cytokines”, “Cytoskeleton”, “Negative regulation of inflammatory response” and “Positive regulation of inflammatory response”. Then, we visualized the expression of these 37 immune-related differentially expressed genes between FSMSCs and HuMSCs with a heatmap plot (Fig. 4e). Second, we compared the gene expression among different MSC subsets and filtered 84 differentially expressed genes in immune MSCs with adjusted p values ≤ 0.05 and absolute log2FC values ≥ 0.8 (Additional file 9). Then, we also intersected the immune-related genes with the 84 differentially expressed genes and finally obtained a total of 14 immune-related differentially expressed genes. We visualized the expression of these 14 intersecting genes among different MSC subsets (Fig. 3g).
Gene regulatory networks analysis
A gene regulatory network was constructed, and 22 cell type specific regulons (CTSRs) were identified in immune MSCs. The immune-related regulons IRF2 and IRF4 were highly expressed in immune MSCs (Fig. 5a). Interestingly, IRF2 and IRF4 might regulate the expression of target gene CXCL12 based on the integrated data (Fig. 5b). Finally, we showed the expression of these 22 CTSRs among various MSC subsets with a heatmap plot (Fig. 5c).
Cell–cell communication analysis
We performed cell–cell communication analysis to explore the interactions among different MSC subsets. We identified a total of 2309 ligand–receptor pairs based on the integrated data (Fig. 6a). Moreover, these ligand–receptor pairs were mapped to 47 signaling pathways (Fig. 6d). We also observed the interaction of CXCL signaling pathways in all MSC subsets and determined the role of each cell subset at the signaling pathway level by calculating the network centrality index (Fig. 6b, c). Furthermore, the mediator score of immune MSCs was high, but the influence score was low, suggesting that immune MSCs were not a network node in the CXCL signaling pathway. However, both the sender and receiver scores of immune MSCs were high, indicating that immune MSCs played an autocrine role in the CXCL signaling pathway. In addition, we calculated the ligand–receptor pairs that contributed to the CXCL signaling pathway. CXCL12-ACKR3 made the most significant contribution to the CXCL signaling pathway. This result suggests that CXCL12-ACKR3 might play an autocrine role in immune MSCs. Finally, we quantified the functional similarity of all 47 signaling pathways and identified four signaling pathway groups (Fig. 6d).
Public database analysis
The proportions of different MSC subsets were assessed among ADMSCs, BMSCs, FSMSCs and HuMSCs. Immune MSCs accounted for 24.36%, 43.46%, 55.65%, and 9.99% of cells among ADMSCs, BMSCs, FSMSCs, and HuMSCs respectively. Progenitor proliferative MSCs accounted for 43.25%, 31.49%, 22.42%, and 17.83% of ADMSCs, BMSCs, FSMSCs, and HuMSCs respectively. Proliferative MSCs accounted for 31.67%, 8.17%, 18.35%, and 67.88% of cells in ADMSCs, BMSCs, FSMSCs, and HuMSCs respectively. Pericyte accounted for 0.72%, 16.89%, 5.58%, and 4.3% of ADMSCs, BMSCs, FSMSCs, and HuMSCs respectively. (Fig. 6f).