Single Cell Transcriptome Sequencing Reveals the Potential Mechanism of Heterogeneity in Immunoregulatory Function Between Mesenchymal Stromal Cells


 Background: Mesenchymal stromal cells (MSCs) could be applied for the treatment of immune-related diseases. However, some studies have found there is enormous heterogeneity in immunomodulatory function of MSCs isolated from different tissue. At present, the underlying mechanism of heterogeneity in immunoregulatory function is still unclear. Methods: In this study, the foreskin mesenchymal stromal cells (FSMSCs) and human umbilical cord mesenchymal stromal cells (HuMSCs) were isolated and cultured to the 3rd passage. Cell types were confirmed according to the standard of International Association for Cell Therapy. Then, FSMSCs and HuMSCs were co-cultured with human peripheral blood mononuclear cells (PBMC) stimulated by lipopolysaccharides (LPS) in viro respectively. And the supernatant was sampled for enzyme-linked immunosorbent assay to investigate the secretion of IL-1β, IL-6, IL-10, TNF-α and TGF-β. Finally, single cell transcriptome sequencing was performed in order to elucidate the mechanism for the difference of immunomodulatory function.Results: FSMSCs and HuMSCs are successfully identified as MSCs. When co-cultured with LPS pre-treated PBMC, FSMSCs and HuMSCs could effectively reduce the secretion of IL-1β and TNF-α. But FSMSCs were able to stimulate the PBMC to secrete more IL-10, TGF-β and IL-6. Furthermore, 4 MSCs subsets in integrated data were identified, including Proliferative MSCs , Pericyte, Immune MSCs and Progenitor Proliferative MSCs. Among them, the proportions of Immune MSCs in FSMSCs and HUMSCs were 56% and 10% respectively. Varieties of immune-related genes, gene sets and regulons were enriched in Immune MSCs. And Immune MSCs with powful transcriptional activity were found to be near to Pericyte at the degree of differentiation and closed to other cell subsets. Finally, the foreskin tissue might be an ideal source of isolating Immune MSCs when comparing the subset composition of MSCs derived from adipose tissue and bone marrow from public database. Conclusions: Immune MSCs may play a key role in the heterogeneity of immunoregulatory function. It is a new insight that Immune MSCs could be isolated and better applied for the treatment of immune-related diseases without being limited by the heterogeneity of immunomodulatory function derived from different tissues.

MSCs, which exist in varieties of tissues, are a kind of stromal cells with multi-directional differentiation potential. MSCs express surface markers CD73, CD90 and CD105, but do not express surface markers HLA-DR, CD11b, CD19, CD34 and CD45 [1]. MSCs can differentiate into different types of cells including osteoblasts, adipocytes, chondrocytes, nerve cells, germ cells and so forth. In addition to the huge potential for tissue differentiation and regeneration, MSCs with the powerful immunoregulatory function are also able to applied for the treatment of multiple immune-related diseases, such as graft-versus-host disease (GVHD) [2], multiple sclerosis [3], sepsis [4], systemic lupus erythematosus [5], crohn disease[6], osteoarthritis [7] and acute lung injury with COVID-19 infection [8]. At present, MSCs are already demonstrated to exert immunomodulatory function by direct cell-cell contact and the secretion of immune regulatory molecules [9]. MSCs were rst isolated from bone marrow. But the clinical application of MSCs derived from bone marrow is greatly limited due to ethical disputes, invasive operation and low yield cultured in vitro. Therefore, it is urgent to nd a new source of MSCs. In recent years, MSCs have been successfully obtained from other tissues such as adipose tissue [10], foreskin tissue [11], umbilical cord [12] and synovial uid [13]. Among them, foreskin tissue and umbilical cord are considered as biological wastes without ethical issues. Meanwhile, foreskin tissue and umbilical cord can be donated from circumcision or natural birth in which the operations are less invasive and do not cause other dangers to the donors. In addition, FSMSCs and HuMSCs have not only high clone proliferation potential in vitro but also huge immunomodulatory ability [11,14]. Therefore, foreskin tissue and umbilical cord are considered to replace bone marrow and become the main source of MSCs. However, the difference between FSMSCs and HuMSCs remains elusive.
Previous studies have shown MSCs with powerful immunomodulatory function are able to be applied in the therapy of lots of immune-associated diseases [2][3][4][5][6][7][8]. However, the clinical application of MSCs still encounters different problems, in which the heterogeneity of MSCs is a major concern [9]. The heterogeneity of immunomodulatory function plays a key role in the heterogeneity of MSCs. It has been shown the difference of immunoregulatory capacity is signi cant between MSCs derived from different tissues [14]. For example, Ivanova-Todorova et al found that the immunomodulatory capacity of adipose MSCs was stronger than bone marrow MSCs [15], whereas Xishan et al found the opposite[16], while Yoo et al showed that the adipose MSCs and bone marrow MSCs had equal immunomodulatory capacity [17]. Even MSCs isolated from the same tissue may exhibit mixed results in immunotherapy. For example, bone marrow MSCs can both inhibit arthritis by upregulating IL- 10[18] and promote arthritis by upregulating IL-6 [19]; some studies have shown that bone marrow MSCs can signi cantly alleviate the clinical symptoms of GVHD [20] while some studies have reported bone marrow MSCs are not effective in improving the outcome of GVHD [21,22]. These controversies indicate that there is great heterogeneity in the immunomodulatory function of MSCs. And the potential mechanism of the difference in immune regulation function of MSCs from different sources is still unclear.
The main intervention for the heterogeneity of immunomodulatory function of MSCs is to treat MSCs with IFN-γ and/or TNF-α in virto before application, this intervention can differentiate MSCs into antiin ammatory MSCs [23]. This is based on the theory of immunomodulatory plasticity that MSCs can be induced to differentiate into pro-in ammatory MSCs or anti-in ammatory MSCs under different immune microenvironment [24]. Some researches have shown that the heterogeneity is signi cantly reduced and immunomodulatory capacity is signi cantly enhanced when MSCs are pretreated with IFN-γ and/or TNFα [25]. However, the pretreatment will also promote the expression of immunogenicity-related genes in MSCs that are unfavorable for clinical application [25]. Therefore, we need to explore new ways to solve the problem of immunomodulatory heterogeneity. Recent evidence suggests that NES + MSCs can signi cantly reduce macrophage in ltration as well as induce macrophage conversion to an anti- in ammatory M2 phenotype[26]; CD271 + MSCs can signi cantly inhibit T lymphocyte proliferation [27] and have stronger chondrogenic potential[28]; CD106 + MSCs can more effectively regulate helper T cells and secrete more immune-related cytokines [29]; CD146 + MSCs can inhibit the activation of Th17 cells [30] while they also promote the conversion of M1-type macrophages to M2-type and improve the in ammation or brosis of the knee [31]. These ndings indicate that MSCs are mixed cell population, and gene markers can help us to distinguish cell subsets with different biological functions in MSCs. A single gene marker may not fully de ne the biological functions of MSCs subsets. Thus, we need to integrate or combine multiple gene markers to distinguish cell subsets with different biological functions in MSCs. Single-cell transcriptome sequencing, an emerging technology, can detect the difference of gene expression between cells and help us to explore MSCs subsets with different biological functions [32]. In this study, we compared the difference between FSMSCs and HuMSCs in immunomodulatory functions in vitro. MSCs subsets with different biological functions in FSMSCs and HuMSCs were identi ed by single cell transcriptome sequencing. Finally, the potential mechanism of the difference in immunomodulatory function between FSMSCs and HuMSCs was investigated.

Methods
Isolation, culture and expansion of mesenchymal stromal cell The foreskin tissue and umbilical cord were collected in the sterile tube containing sterile phosphate buffer solution (PBS) respectively. Under sterile conditions, the foreskin tissue and umbilical cord were cleaned with iodophor and PBS in turn, and then transferred to a petri dish. Next, we separated dermis tissue from foreskin tissue and cut up the dermis tissue, and then transfer in a 10-cm petri dish. In umbilical cord, we separated the amniotic membrane, artery and vein. And then, we collected the white connective tissue between the amniotic membrane and blood vessel in a 10-cm petri dish too. Afterwards, the minced dermis tissue and white connective tissue were cultured in incubator (37℃, 5% CO 2 ) with 5 mL of DMEM/F12 (Gibco, USA) containing 10% fetal bovine serum (Gibco, USA) respectively until primary cells migrated from tissue fragments and attached to the plate. When cells became 80% con uent, they were trypsinized and passaged. Finally, cells from the 3rd passage were used for further study and cell morphology was observed under an inverted microscope (Leica, Germany).

The surface markers of mesenchymal stromal cell
According to the identi cation standard about MSCs of the International Association for Cell Therapy (ISCT) [1], we detected the expression of CD73, CD90, CD105, CD45, CD34, CD14, CD19 and HLA-DR in the third passage of HuMSCs and FSMSCs by ow cytometry. First, we discarded medium and washed cells twice with PBS, and then used trypsin to digest the cells. Next, we stopped the digestion and adjusted the cell concentration to 1x10 6 cells/mL. We prepared 8 tubes and added 1 mL of suspended cells in each tube. Afterwards, we added 5 uL different antibodies, such as mouse anti-human FITC-CD90 (BD, Germany), mouse anti-human PerCP-CyTM5.5-CD105 (BD, Germany), mouse anti-human APC-CD73 (BD, Germany), mouse anti-human PE-HLA-DR (BD, Germany), mouse anti-human PE-CD11b (BD, Germany), mouse anti-human PE-CD19 (BD, Germany), mouse anti-human PE-CD45 (BD, Germany) and mouse antihuman PE-CD34 (BD, Germany), in different tubes and incubated them at 37℃ for 30 minutes in dark.
After incubation, we added 1 mL of PBS to wash out the excess antibodies and then centrifuged cells at 250 g for 10 minutes. Finally, the cells were resuspended in 0.5 mL of PBS and analyzed in Canto II cytometer (BD Bioscience).

Comparsion of immunomodulatory function
HuMSCs and FSMSCs were cultured to 3rd passage in vitro. The cells in the logarithmic growth phase were digested with trypsin and then seeded in 12-well plates (1x10 6 cells per well). The cells were cultured at 37℃ until the con uence reaches 50%. Then, we discarded the supernatant and replaced with 1 mL fresh serum-free medium. PBMC were isolated by Ficoll-Hypaque gradient centrifugation from whole blood the health donor providing with informed consent. At the same time, we resuspend the pbmc to 1x10 6 cells/mL in serum-free medium and then PBMC were mixed with MSCs at a 10:1 ratio during the co-culture. The experiment was divided into 4 groups: a) PBMC/HuMSCs/FSMSCs cultured alone, b) PBMC/HuMSCs/FSMSCs stimulated by LPS, c) HuMSCs/FSMSCs co-cultured with PBMC without stimulation, d) HuMSCs/FSMSCs co-cultured with PBMC stimulated by LPS. And the nal concentration of LPS (Sigma) was 10 ng/mL. Cells were incubated for 24 hours (37℃, 5% CO 2 ). Then, the supernatant was centrifuged (1000 g, 5 min) and collected at 2h, 4h, 12h and 24h respectively. Finally, we detected the secretion levels of IL-1β, IL-6, IL-10, TNF-α and TGF-β in the supernatant according to the manual of enzyme-linked immunosorbent assay (ELISA) kit (NeoBiosience, China Second, the concentration of cell suspension were adjusted to 1×10 6 cells/ml for 10X Genome Single Cell Platform (Single Cell 3' library and Gel Bead Kit V3). Then, the samples were processed as follows: GEMs creation, thermal cycling, post cycling cleanup, cDNA amplication, library preparation, library quanti cation and sequencing on Illumina Hiseq3000. Finally, the result of sequencing was mapped to GRCh38 human genome, and reads were counted by cellranger software (version 3.1.0) [33].

Data processing
We performed cell quality control and cluster analysis for single cell transcriptome data by R package Seurat (version 4.0.0) [34]. First, we used the Read10X function of R package Seurat (version 4.0.0) [34] to load the UMI count matrix into R and create the Seurat objec. Secondly, we ltered 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, in order to correct the sequencing depth of every cell, we ran normalization for single cell transcriptome data and removed the effect of mitochondrial gene during the process of scaling. Afterwards, we executed principal component analysis (PCA) through rst 2000 genes with high heterogeneity and evaluated the rst 50 PCA dimensions by JackStraw and ElbowPlot function of R package Seurat (version 4.0.0) [34]. And the rst 20 PCA dimensions with the most signi cant p values were used for t-distributed stochastic neighbor embedding (TSNE) and uniform manifold approximation and projection (UMAP). We also used the rst 20 principal components (PCs) to cluster cells, and the resolution of the clustering was 0.5. Finally, we annotated the clusters based on the expression results of the surface markers of MSCs [1].

Integration analysis
First, we integrated the subsets annotated as MSCs in FSMSCs and HuMSCs by the   [40], for which we were able to assess the potential biological function for every subset.

Differential gene expression analysis
We compared the gene expression between FSMSCs and HuMSCs by the Findallmarker function of R package Seurat (version 4.0.0) [34] and lter the differential genes with the corrected p value ≤ 0.05 and the absolute value of log 2 fold change (Log 2 FC) ≥ 2. Meanwhile, we also calculated the marker genes in different MSCs subsets by the Findallmarker function of R package Seurat (version 4.0.0) [34] with the corrected P value ≤ 0.05 and Log 2 FC ≥ 0.8. In order to analyze the potential immune function of these highly expressed genes, we extracted 2073 immune-related genes from ImmPort database [41] and MSigDB database (Gene Ontology database) [42]. Then, we intersected the immune-related genes with the differential genes and annotated these intersected genes in immune function according to the related literature, ImmPort database and GeneCards database [43]. Finally, we visualized the expression of these intersected genes by R packages ComplexHeatmap (version 2.6.2) [44] and circlize (version 0.4. 12) [45].
Gene regulatory networks analysis pySCENIC (version 0.10.3) is a python tool which enables us to infer single cell gene regulatory networks from single cell RNA-seq data [46]. The work ow was as follows: a). We inferred the co-expression modules containing different transcription factors and their target genes from single-cell RNA-seq matrix; b). We retained the signi cantly enriched motifs in the co-expression module by motif enrichment analysis, and removed the target genes which were less correlated with the signi cant motif in the coexpression module. Then, the transcription factor and remaining target genes of co-expression module were called regulon; c). We assessed the regulon activity score (RAS) of different regulons in each cell, as well as calculated the RAS threshold for each regulon. This regulon was considered to be activated in cell when RAS was greater than threshold and vice versa. We performed the "0/1" conversion for RAS matrix based on threshold and created binary matrix in order to show the difference of regulon in various cell types best and eliminate technical bias; d). We calculated the regulon speci c score (RSS) in different MSCs subsets by using R package philentropy (version 0.4.0) [47]. The regulons with RSS > 0.6 were considered as cell-type-speci c regulons (CTSRs).

Cell-Cell communication analysis
CellChat (version 1.0.0) is an R package that can infer and visualize intercellular communication network from single-cell RNA-seq data [48]. The work ow was as follows: a). We created a CellChat object and set to human ligand-receptor interactions database; b). We identi ed overexpressed ligands or receptors and projected gene expression data onto the protein-protein interaction network; c). We assigned a probability value to each ligand-receptor interaction and then performed permutations to infer biologically meaningful interactions; d). We summarized the probabilities of all interactions associated with each signaling pathway and calculated the probability of communication for each signaling pathway; e). We identi ed the roles of each cell subsets at signaling pathways level by calculating network centrality; f).
We nally identi ed several signaling groups through quantifying the functional similarity of all signi cant signaling pathways.

Comparsion of MSCs isolated from different tissue sources
In order to evaluate the ideal tissue source for isolating Immune MSCs, we downloaded single-cell transcriptome data from the NCBI public database for adipose mesenchymal stromal cells (ADMSCs) (SRP148833) [10] and bone marrow mesenchymal stromal cells (BMSCs) (GSE115149, GSE162692) [49,50]. Among them, SRP148833 database was processed by 10x Genomics cellranger software (version 3.1.0) [33] for which we were able to acquire the gene expression matrix of ADMSCs. GSE115149 and GSE162692 database only kept the gene expression matrix of untreated BMSCs. After obtaining the gene expression matrix, they were processed according to the same quality control process as above. Finally, we evaluated the proportions of Immune MSCs in ADMSCs and BMSCs based on the integrated data above by R package SingleR (version 1.4.0) [51].

Identi cation of FSMSCs and HuMSCs
According to the minimum identi cation standard about MSCs of the ISCT, we con rmed 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 broblast-like and long spindle shaped under the microscope. Secondly, the result of ow 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 identi cation 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 in ammation-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 signi cant at any time point; The basal secretion level of IL-6 was signi cantly higher in HuMSCs than that in FSMSCs and PBMC while it showed no signi cant differences between FSMSCs and PBMC; The basal secretion level of TNF-α was signi cantly higher in PBMC than that in FSMSCs and HuMSCs, while it showed no signi cant differences between FSMSCs and HuMSCs. In addition, we also compared whether there were differences in basal secretion levels of in ammation-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 uctuating pattern; The differences in the secretion of IL-1β, IL-10 and TNF-α were not statistically signi cant at different time points. In summary, FSMSCs and HuMSCs, like PBMC, were able to secrete in ammation-associated cytokines.
Next, we compared the differences in the secretion of in ammation-associated cytokines between FSMSCs, HuMSCs and PBMS under lipopolysaccharides (LPS) stimulation. The differences in secretion levels of IL-10 and TGF-β were not statistically signi cant; 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 in ammation-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 signi cant; 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 proin ammatory cytokines IL-1β and TNF-α were signi cantly increased in PBMCs after stimulation by LPS. However, FSMSCs and HuMSCs stimulated by LPS showed no signi cant changes in the secretion levels of pro-in ammatory cytokines IL-1β and TNF-α.
Moreover, we compared the differences in the secretion of in ammation-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 signi cant difference in the secretion level of IL-1β among these three conditions. In addition, we compared whether there were differences in the secretion of in ammation-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 in ammatory 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 in ammation-related cytokines (Fig. 2d, Supplementary Table 4). The secretion levels of the pro-in ammatory 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-in ammatory factor IL-10 and TGF-β were higher in FSMSCs than that in HuMSCs. In addition, we investigated whether there were signi cant differences in the secretion levels of in ammation-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-in ammatory cytokines IL-1β and TNF-α from PBMC. It implied that FSMSCs and HuMSCs were both able to regulate immune function. Meanwhile, FSMSCs were more signi cantly able to upregulate the secretion of the bidirectional regulator in ammation IL-6 and the anti-in ammatory 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 in ammation-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 ltered 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 rst 2000 genes with high heterogeneity and evaluated the rst 50 PCs. The rst 20 PCs with the most signi cant p values were used for TSNE and UMAP. We also used the rst 20 PCs to cluster cells with a resolution of 0.5, and we nally obtained 7 clusters in FSMSCs and 10 clusters in HuMSCs respectively (Fig. 3a).

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 attened, and small and prototypical. The rst 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)

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 de ned Pericyte as the starting point of the differentiation trajectory of MSCs. Next, we constructed and nally gained two differentiation trajectories (Fig. 4a). The rst trajectory started from Pericyte and differentiated into Proliferative MSCs; the second trajectory started from Pericyte, passed through Immune MSCs and Progenitor Proliferative MSCs, and nally 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 in uenced 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 nd 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 ltered 15 statistically signi cant gene sets in Immune MSCs with the adjusted P value ≤ 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 signi cantly down-regulated in Proliferative MSCs while it was signi cantly 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 signi cantly 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 ltered 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 in ammatory response" and "Positive regulation of in ammatory 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-in ammatory 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 ltered 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 identi ed twenty-two cell type speci c 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-speci c 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 in ammatory 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 identi ed 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 In uence 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 quanti ed the functional similarity of all 47 signal pathways and identi ed 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.

Discussion
MSCs are a kind of stromal cells in varieties of tissues, which they are proved to possess the immunomodulatory function and can be used to treat lots of immune diseases. However, MSCs derived from different tissue are signi cantly heterogeneous in immunomodulatory function. At present, the major way to reduce the heterogeneity is pre-treating MSCs with IFN-γ and/or TNF-α. But it also increases the immunogenicity of MSCs, accelerate the clearance of MSCs in vivo, as well as enhances the potential risk of immune rejection. This is detrimental to the clinical application of MSCs. In this study, we compared the immunomodulatory function between FSMSCs and HuMSCs. We found FSMSCs, which had better immunomodulatory capacity, are more suitable for the treatment of in ammation than In addition, we assessed the proportion of Immune MSCs between ADMSCs, BMSCs, FSMSCs and HuMSCs, for which the foreskin tissue might be an ideal source of isolating Immune MSCs. We believe MSCs are immunoplastic and susceptible because of external environmental in uences. And the foreskin tissue is usually exposed to bacterial and viral attack due to its surrounding dirty environment. Moreever, foreskin tissue is prone to undergo ischemia or hypoxia because of insu cient blood supply. In previous studies, hypoxia is proved to promote the immunoregulatory capacity of MSCs [59]. Thus, MSCs derived from foreskin tissue are more likely to differentiate toward anti-in ammatory MSCs which may be comprised mainly of immune MSCs in our study. In summary, there is a great heterogeneity in the MSCs derived from different tissues. An the Immune MSCs may be the potential mechanism of heterogeneity in immunoregulatory function through single cell transcriptome sequencing. Therefore, our study provides a new insight to explain and address the immunomodulatory heterogeneity of MSCs.

Conclusions
In summary, in the study, FSMSCs and HuMSCs are successfully isolated and identi ed. We compare the immunoregulatory function between FSMSCs and HuMSCs in viro. Then, the single cell transcriptome sequencing is performed and we nd Immune MSCs may play a key role in the heterogeneity of immunoregulatory function through bioinformatics analysis. It is a new insight that Immune MSCs could be isolated and better applied for the treatment of immune-related diseases without being limited by the heterogeneity of immunomodulatory function derived from different tissues.     The results of trajectory inference, RNA velocity analysis, gene set enrichment analysis and differential gene expression analysis.

Supplementary Files
This is a list of supplementary les associated with this preprint. Click to download.