Identification of FSMSCs and HuMSCs
According to the minimum identification standard about MSCs of the ISCT, we confirmed the type of cells which were isolated from foreskin tissue and umbilical cord respectively. First, FSMSCs and HuMSCs were able to attach plastic surfaces, and they were observed to be typical fibroblast-like and long spindle shaped under the microscope. Secondly, the result of flow cytometry showed positive surface markers, such as CD73, CD90 and CD105, were both expressed more than 95%, while negative surface markers, such as CD45, CD34, CD14, CD19 and HLA-DR, were expressed less than 2% in FSMSCs and HuMSCs (Figs. 1a-1b). Then, the red spherical vacuoles could be observed by Oil Red-O staining after the induction of adipogenic differentiation in FSMSCs or HuMSCs (Fig. 1c). It meant the lipid droplets were formed in these cells. In addition, after the induction of osteogenic differentiation, we observed scattered, crimson and round nodules in FSMSCs or HuMSCs by 2% ARS staining (Fig. 1c). It revealed that massive calcium was deposited in the cytoplasm. Finally, the blue acid mucopolysaccharides were visible in FSMSCs or HuMSCs by Alcian Blue staining after the induction of chondrocytic differentiation (Fig. 1c), suggesting that the FSMSCs and HuMSCs both had chondrogenic ability. Thus, these assays altogether indicated that the cells isolated from forskin tissue or umbilical cord were mesenchymal stromal cells according to the identification standard of ISCT[1].
Immunomodulatory function of FSMSCs and HuMSCs
In order to compare the immunomodulatory function between FSMSCs and HuMSCs, the secretion levels of IL-1β, IL-6, IL-10, TNF-α and TGF-β were measured through ELISA.
First, we compared the differences in autocrine inflammation-related cytokines among FSMSCs, HuMSCs and peripheral blood mononuclear cells (PBMC) (Fig. 2a, Supplementary Table 1). The differences in basal secretion levels of IL-1β, IL-10 and TGF-β in these three cells were not statistically significant at any time point; The basal secretion level of IL-6 was significantly higher in HuMSCs than that in FSMSCs and PBMC while it showed no significant differences between FSMSCs and PBMC; The basal secretion level of TNF-α was significantly higher in PBMC than that in FSMSCs and HuMSCs, while it showed no significant differences between FSMSCs and HuMSCs. In addition, we also compared whether there were differences in basal secretion levels of inflammation-related cytokines among these three cell types at different time points (Fig. 2a, Supplementary Table 1). The basal secretion level of IL-6 increased with time while the basal secretion level of TGF-β showed fluctuating pattern; The differences in the secretion of IL-1β, IL-10 and TNF-α were not statistically significant at different time points. In summary, FSMSCs and HuMSCs, like PBMC, were able to secrete inflammation-associated cytokines.
Next, we compared the differences in the secretion of inflammation-associated cytokines between FSMSCs, HuMSCs and PBMS under lipopolysaccharides (LPS) stimulation. The differences in secretion levels of IL-10 and TGF-β were not statistically significant; The secretion levels of IL-1β and TNF-α were higher in PBMCs while the secretion level of IL-6 was higher in FSMSCs. In addition, we compared whether there were differences in the secretion level of inflammation-associated cytokines among MSCs at different time points after stimulating with LPS (Fig. 2b, Supplementary Table 2). Among them, the secretion level of IL-6 all increased with time while the differences in secretion level of IL-10 was not statistically significant; The secretion levels of IL-1β and TNF-α in PBMC increased with time; The secretion level of TGF-β in FSMSCs increased with time. In summary, the secretion levels of pro-inflammatory cytokines IL-1β and TNF-α were significantly increased in PBMCs after stimulation by LPS. However, FSMSCs and HuMSCs stimulated by LPS showed no significant changes in the secretion levels of pro-inflammatory cytokines IL-1β and TNF-α.
Moreover, we compared the differences in the secretion of inflammation-related cytokines between FSMSCs and HuMSCs co-cultured with PBMC respectively and PBMC cultured alone (see Fig. 2c, Supplementary Table3). The secretion levels of IL-6, TNF-α, IL-10, and TGF-β were overall higher in FSMSCs co-cultured with PBMC than in HuMSCs co-cultured with PBMC or in PBMC cultured alone; There was no significant difference in the secretion level of IL-1β among these three conditions. In addition, we compared whether there were differences in the secretion of inflammation-related cytokines at different time points among these three conditions (Fig. 2c, Supplementary Table 3). Among them, the secretion level of IL-6 increased with time; The secretion levels of IL-10 and TNF-α in FSMSCs co-cultured with PBMC increased with time. In summary, FSMSCs were able to stimulate PBMC to secrete TNF-α, IL-6, IL-10 and TGF-β than the HuMSCs.
Finally, we mimicked the inflammatory response by stimulating PBMC with LPS in vitro. FSMSCs and HuMSCs were co-cultured with PBMC pre-stimulation respectively. And then, we compared the immunoregulatory ability through detecting the secretion of inflammation-related cytokines (Fig. 2d, Supplementary Table 4). The secretion levels of the pro-inflammatory cytokines IL-1β and TNF-α were lower in the FSMSCs or HuMSCs than that in the control; The secretion levels of the bidirectional immunomodulatory cytokines IL-6, anti-inflammatory factor IL-10 and TGF-β were higher in FSMSCs than that 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, Supplementary Table 4). Among them, the secretion levels of IL-6, IL-10 and TGF-β increased with time in FSMSCs co-cultured with PBMC pre-stimulation; The secretion level of IL-10 increased with time in HuMSCs co-cultured with PBMC pre-stimulation. In summary, FSMSCs and HuMSCs both inhibited the secretion of pro-inflammatory cytokines IL-1β and TNF-α from PBMC. It implied that FSMSCs and HuMSCs were both able to regulate immune function. Meanwhile, FSMSCs were more significantly able to upregulate the secretion of the bidirectional regulator inflammation IL-6 and the anti-inflammatory cytokines IL-10 and TGF-β than HuMSCs. It indicated that FSMSCs showed stronger immunomodulatory ability than HuMSCs, which was more suitable for the treatment of inflammation-related diseases.
Quality control, cluster and biological annotation
The quality control, cluster and biological annotation were used for FSMSCs and HuMSCs through bioinformatics analysis. First, we created the Seurat object for FSMSCs and HuMSCs separately. Secondly, we performed quality control and filtered the low-quality cells with less than 200 genes, more than 20000 genes, more than 10% percentage of mitochondrion genes and more than 10% percentage of erythrocyte genes. Then, we gained 7335 cells in FSMSCs in which the median of UMI was 33202 and the median of gene was 5440. We also gained 12542 cells in HuMSCs in which the median of UMI was 15191 and the median of gene was 3831. The result showed that the quality of the transcriptome sequencing data was high. Afterwards, we ran normalization and removed the effect of mitochondrial gene in the process of scale. Next, we executed dimensionality reduction by PCA through first 2000 genes with high heterogeneity and evaluated the first 50 PCs. The first 20 PCs with the most significant p values were used for TSNE and UMAP. We also used the first 20 PCs to cluster cells with a resolution of 0.5, and we finally obtained 7 clusters in FSMSCs and 10 clusters in HuMSCs respectively (Fig. 3a).
Then, we annotated the clusters according to the surface markers of MSCs[1]. In FSMSCs, C0, C1, C2, C3 and C5 highly expressed the positive marker genes CD73, CD90 and CD105, and lacked the expression of the negative marker genes HLA-DRA, CD11b, CD19, CD34 and CD45 (Fig. 3a). It meant that C0, C1, C2, C3 and C5 were possible MSCs. But C4 highly expressed HLA-DRA and C6 was lack of the expression of CD73 and CD105. It indicated that C4 and C6 did not satisfy the identification criteria for MSCs. Similarly, in HuMSCs, C0, C1, C2, C3, C5 and C6 highly expressed the positive marker genes CD73, CD90 and CD105, and lacked the expression of the negative marker genes HLA-DRA, CD11b, CD19, CD34 and CD45, suggesting that C0, C1, C2, C3, C5 and C6 might be MSCs. But C4, C7, C8 and C9 were lack of the expression of CD73 and CD105, therefore C4, C7, C8 and C9 did not satisfy the identification criteria for MSCs (Fig. 3a). In summary, these results suggested that only C0, C1, C2, C3 and C5 were MSCs in FSMSCs and accounted for 93.05% of FSMSCs. The C0, C1, C2, C3, C5 and C6 were MSCs in HuMSCs and accounted for 95.85% of HuMSCs.
Integration anlysis
In order to reveals the potential heterogeneity of MSCs, we integrated the subsets which were annotated as MSCs in FSMSCs and HuMSCs. Then, we evaluated the batch effect of integration data by PCA plot in which the cells of FSMSCs and HuMSCs were able to overlap well after integration (Fig. 3d), indicating that the batch effect was corrected well. Afterwards, the integrated data were re-normalized, downscaled and clustered (resolution=0.5) for which we gained 7 major clusters in integration data (Fig. 3e). At the same time, we calculated the proportion of these clusters between FSMSCs and HuMSCs and showed them on bar plot (Fig. 3f).
Moreover, we inferred the cell cycle of each cluster by the CellCycleScoring function of R package Seurat (version 4.0.0)[34]. There was a higher proportion of G1 phase cells in C0, C2 and C6, suggesting a resting state. But G2M or S phase cells occupied a large proportion in C2, C3, C4 and C5 (Fig. 3f). In addtion, MKI67, a gene encodes a nuclear protein associated with cell proliferation[52], was highly expressed in C2, C3, C4 and C5 (Fig. 3b). It suggested that the proliferative activity of C2, C3, C4 and C5 was stronger than other subsets. Intriguingly, the result of cellcycle inference were consistent with the result of MSCs cultured in vitro. It had been reported that there are three morphologies of MSCs in vitro, which are elongated spindle-shaped, large and flattened, and small and prototypical. The first two kinds of cells are larger and proliferated more slowly. The latter cells are smaller and proliferated more rapidly, which are known as rapid self-renewing cells[53, 54]. It meant there were some cells with more active proliferation in MSCs. Therefore, we named C2, C3, C4 and C5 as Proliferative MSCs.
Secondly, C6 highly expressed the positive marker genes CD146, PDGFRB and NG2, and lacked the expression of the negative marker genes CD31, CD45 and CD34 (Fig. 3b). These genes were proved to be the markers of pericyte which could expressed the markers of MSCs too in prior study[55]. Therefore, C6 were possible pericyte. In addition, CXCL12 and PTGIS were highly expressed in C0 while MKI67 and CD146 were lack of the expression in C0 (Fig. 3b). Among them, The CXCL12 was known to play an important role in immunoregulatory function of MSCs[56]. PTGIS, prostacyclin I2 synthase, is able to catalyze the conversion of prostaglandin H2 to prostaglandin I2, which has been demonstrated to play an important role in immunoregulatory function according to previous research[57]. Meanwhile, we also observed that C1 and C6 could express CXCL12 and PTGIS, but their expression level was lower than C0. Its suggested that C0, C1 and C6 both possessed immunomodulatory ability. But the immunomodulatory function of C0 was most powerful among them. Therefore, we named C0 as Immune MSCs.
Finally, C1 lowly expressed CXCL12 and PTGIS, and lacked expression of MKI67 and CD146 (Fig. 3b). In addition, a few cells in C1 could express MKI67 (Fig. 3b). And we constructed a clustering tree based on the similarity of MSCs subsets as well as visualized it on dendrogram plot. Then, the expression patterns of C1 and C5 were closer (Fig. 3c). In addition, C5 belonged to Proliferative MSCs. It revealed that C1 were possible the precursor cells of C5 while C1 might be a transitional stage from C0 to Proliferative MSCs. Therefore, we named C1 as Progenitor Proliferative MSCs.
In summary, there were 4 subsets in the integrated data (Fig. 3e), such as Proliferative MSCs (MKI67+, CD146low+, NG2+, PDGFRB-), Pericyte (CD146high+, PDGFRB+, MKI67-, CD31-, CD45-, CD34-), Immune MSCs (CXCL12high+, PTGIShigh+, PDGFRB+, CD146-, MKI67-) and Progenitor Proliferative MSCs (CXCL12low+, PTGISlow+, PDGFRB+, CD146-, MKI67-). In addition, we calculated the proportion of different MSCs subsets between FSMSCs and HuMSCs (Fig. 3f). Among them, Proliferative MSCs made up 16% of FSMSCs and 68% of HuMSCs; Pericyte made up 6% of FSMSCs and 4% of HuMSCs; Immune MSCs made up 56% of FSMSCs and 10% of HuMSCs; Progenitor Proliferative MSCs made up 22% of FSMSCs and 18% of HuMSCs. In previous analysis, Pericyte, Progenitor Proliferative MSCs and Immune MSCs could express CXCL12 and PTGIS. It indicated that Pericyte, Progenitor Proliferative MSCs and Immune MSCs all possessed immunoregulatory capacity. However, CXCL12 and PTGIS in Immune MSCs had the highest expression, suggesting Immune MSCs might have stronger immunoregulatory ability. Then, we also observed that the difference of the proportion of Immune MSCs between FSMSCs and HuMSCs was significant, and most of Immune MSCs were derived from FSMSCs. These results demonstrated that Immune MSCs might played a key role in immunomodulatory function between FSMSCs and HuMSCs.
Trajectory inference
First, we inferred the differentiation degree of different MSCs subsets. The cytotrace score which was near to 1.0 indicated a lower degree of differentiation while the cytotrace score was near to 0.0 indicated a higher degree of differentiation. And the degree of differentiation of Proliferative MSCs and Progenitor Proliferative MSCs were higher than Immune MSCs and Pericyte (Fig. 4b). Some studies had reported that MSCs might originate from Pericyte[58]. Therefore, we defined Pericyte as the starting point of the differentiation trajectory of MSCs. Next, we constructed and finally gained two differentiation trajectories (Fig. 4a). The first trajectory started from Pericyte and differentiated into Proliferative MSCs; the second trajectory started from Pericyte, passed through Immune MSCs and Progenitor Proliferative MSCs, and finally differentiated into Proliferative MSCs. Proliferative MSCs were located at the end of differentiation in both 2 trajectories. And the immunomodulatory capacity of Proliferative MSCs was lower than that of Immune MSCs and Pericyte. These analyses indicated that the immunomodulatory function of MSCs might be influenced by the degree of differentiation.
RNA velocity analysis
The expression of MKI67 in Immune MSCs was lower than that in Proliferative MSCs from above analysis, which meant the proliferative activity of Immune MSCs was much lower. But we find most of arrows in Immune MSCs were longer than that in Proliferative MSCs (Fig. 4c). And the longer the arrow, the stronger the transcriptional activity. It meant that the overall transcriptional activity of Immune MSCs was stronger than Proliferative MSCs. Although Immune MSCs were less differentiated and showed weaker proliferative activity, it still express a large number of genes to maintain its biological function.
Gene set enrichment analysis
In order to assess the potential biological functions, all cell subsets were subjected to gene set enrichment analysis. First, we calculated the enrichment score and adjusted P value of Hallmark gene set in different MSCs subsets (Supplementary Tables 5 and 6) and demonstrated them on heatmap plot (Fig. 4d). Then, we filtered 15 statistically significant gene sets in Immune MSCs with the adjusted P value ≤ 0.05. Among them, the immune-related gene sets "IL6 JAK STAT3 SIGNALING", "INTERFERON GAMMA RESPONSE" and "INTERFERON ALPHA RESPONSE" were up-regulated in Immune MSCs, indicating a strong immune regulatory ability. In contrast, the mitosis-related gene sets "MITOTIC SPINDLE", "G2M CHECKPOINT", "E2F TARGETS" were down-regulated in Immune MSCs, suggesting weak proliferative activity. Secondly, we filtered 13 and 19 statistically significant gene sets by corrected P value ≤ 0.05 in Progenitor Proliferative MSCs and Pericyte respectively. Among them, the immune-related gene set "TNFA SIGNALING VIA NFKB" was upregulated in both Progenitor Proliferative MSCs and Pericyte, suggesting that these MSCs also had immunoregulatory capacity. In addition, Progenitor Proliferative MSCs tended to enrich gene set related to adipose differentiation, such as "ADIPOGENESIS" and "FATTY ACID METABOLISM", thereby Progenitor Proliferative MSCs might possess a stronger lipogenic potential. However, in Pericyte, it preferred to enrich gene set related to mesenchymal histogenesis, such as "MYOGENESIS", "ANGIOGENESIS" and "EPITHELIAL MESENCHYMAL TRANSITION". It indicated that Pericyte had a strong potential for mesenchymal tissue transformation and regeneration. Next, we filter 20 statistically significant gene sets in Proliferative MSCs with corrected P value ≤ 0.05. It lacked the immune-related gene sets "IL6 JAK STAT3 SIGNALING", "INTERFERON GAMMA RESPONSE" and "INTERFERON ALPHA RESPONSE". However, the mitosis-related gene set "MITOTIC SPINDLE", "G2M CHECKPOINT", and "E2F TARGETS" MSCs were significantly up-regulated in Proliferative. It meant that Proliferative MSCs might not possess immune regulatory ability but their activity of proliferation was high. We also found that the hypoxia-related gene set "HYPOXIA" was significantly down-regulated in Proliferative MSCs while it was significantly up-regulated in Immune MSCs and Pericyte. In previous studies, hypoxia was proved to promote the immunoregulatory capacity of MSCs[59]. Thus, it further indicated Immune MSCs and Pericyte may possess strong immunoregulatory capacity.
In summary, "IL6 JAK STAT3 SIGNALING", "INTERFERON GAMMA RESPONSE" and " INTERFERON ALPHA RESPONSE" gene sets were significantly enriched in Immune MSCs (Fig. 4d). "IL6 JAK STAT3 SIGNALING", "INTERFERON GAMMA RESPONSE" and "INTERFERON ALPHA RESPONSE" gene sets were associated with immunomodulatory functions. We concluded that Immune MSCs played an important role in immunomodulatory functions of MSCs.
Differential gene expression analysis
We compared the gene expression between FSMSCs and HuMSCs and filtered 205 differential genes with the adjusted p value ≤ 0.05 and the absolute value of log2 fold change (Log2FC) ≥ 2 (Supplementary Table 7). Then, we intersected the immune-related genes with 205 differential genes, and obtained a total of 37 intersected genes. We annotated the intersected genes for immune functions based on the ImmPort database, Genecard database and relevant literature reports. These intersecting genes could be divided into 7 major categories, which were "Antigen processing and presentation", "Antimicrobials ", "Cytokine receptors", "Cytokines", "Cytoskeleton", "Negative regulation of inflammatory response" and "Positive regulation of inflammatory response". Among them, RARRES2, B2M, LGALS3 and ADM which were highly expressed in FSMSCs were in "Antimicrobials" category and reported to possess immunomodulatory and antimicrobial ability [60-63]. In addition , RARRES2 has been found to have anti-inflammatory ability[64]. Then, we showed the expression of these 37 intersected genes between FSMSCs and HuMSCs (Fig. 4e). Next, we calculated the marker genes in different MSCs subsets and filtered 84 differential genes in Immune MSCs with the adjusted P value ≤ 0.05 and the absolute value of Log2FC ≥ 0.8 (Supplementary Table 8). Then, we intersected the immune-related genes with 84 differential genes, and obtained a total of 14 intersected genes. Among them, BIRC5 and PBK were down-regulated while CXCL12, TNFRSF11B, ADM, LGALS3, NUPR1, PTGIS, IGFBP4, FGF7, B2M, CRLF1, ACKR4 and SERPINF1 were highly expressed in Immune MSCs. B2M, LGALS3 and ADM were proved to possess immunomodulatory and antimicrobial ability in the previous studies[61-63]. It indicated that Immune MSCs might have strong immunomodulatory functions. Finally, we visualized the expression of these 14 intersected genes between different MSCs subsets (Fig. 3g).
Gene regulatory networks analysis
Next, we constructed gene regulatory network and identified twenty-two cell type specific regulons (CTSRs). These CTSRs might play an important role in the function of Immune MSCs. Among them, regulons IRF2 and IRF4, were highly expressed in Immune MSCs (Fig. 5a). IRF2 and IRF4 belong to the interferon regulatory factor family which is associated with interferon regulation in response to viral infection and the regulation of interferon-inducible genes. Usually, members of this family are lymphocyte-specific and can negatively regulate Toll-like receptor signaling. The negative regulation is essential for the activation of innate and adaptive immunity. It had been reported that IRF4 was an inhibitor of TLR-induced inflammatory pathways[65]. Interestingly, IRF2 and IRF4 might regulate the expression of target gene CXCL12 based on single-cell RNA-seq data (Fig. 5b). It further indicates that IRF2 and IRF4 might play a key role in the immunomodulatory function of Immune MSCs. Finally, we showed the expression of all CTSRs among various MSCs subsets on heatmap plot (Fig. 5c).
Cell-Cell communication analysis
Afterwards, we performed Cell-Cell communication analysis to explore the interaction among different MSCs subsets. We identified a total of 2309 ligand-receptor pairs based on single cell transcriptome data (Fig. 6a). And these ligand-receptor pairs were mapped to 47 signaling pathways (Fig. 6e). Among them, there was a close interaction between Immune MSCs and other cell subsets (Fig. 6a). We observed the interaction of CXCL signaling pathways in all MSCs subsets and determined the role of each cell subset at the signaling pathway level by calculating the network centrality index (Figs. 6b-6c). And 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 signal pathway. However, both the Sender score and Receiver score of Immune MSCs were high, indicating that Immune MSCs played an autocrine role in CXCL signaling pathway. Meanwhile, we showed the ligand-receptor pairs which contributed to the CXCL signal pathway on the histogram plot. Among them, CXCL12-ACKR3 made the greatest contribution to the CXCL signal pathway (Fig. 6d). It suggested that CXCL12-ACKR3 might play an autocrine role in in Immune MSCs. ACKR3 is the scavenger receptor of CXCL12, which can regulate the concentration gradient of CXCL12[66]. The expression of ACKR3 in Immune MSCs can prevent excessive CXCL12 from causing damage to the cells themselves. In addition, some studies have found that blocking ACKR3 with antagonist AMD3100 can improve the immunosuppressive microenvironment of tumor[67]. It suggests that ACKR3 not only plays an important role in regulating the concentration gradient of CXCL12, but also maintains the immunosuppressive function. Therefore, CXCL12-ACKR3 might plays an important role in the Immune MSCs. Finally, we quantified the functional similarity of all 47 signal pathways and identified 4 signal pathway groups (Fig. 6e). Among them, CXCL, IGF, MIF, IL10, PARs, HGF, NT, GAS, EDN, ACTIVIN and GRN belonged to the same signal pathway, which implied that these signaling pathways might help CXCL pathway maintain the biological functions of Immune MSCs.
Public database analysis
We assessed the proportion of different MSCs subsets between ADMSCs, BMSCs, FSMSCs and HuMSCs which came from public database and private data. We found that Immune MSCs account for 24.36%, 43.46%, 55.65%, and 9.99% of cells in ADMSCs, BMSCs, FSMSCs, and HuMSCs respectively (Fig. 6f). It suggested that foreskin tissue might be an ideal tissue source for Immune MSCs.