Skip to main content
Download PDF

Deconstructing transcriptional variations and their effects on immunomodulatory function among human mesenchymal stromal cells

Stem Cell Research & Therapy volume 12, Article number: 53 (2021) Cite this article

Abstract

Background

Mesenchymal stromal cell (MSC)-based therapies are being actively investigated in various inflammatory disorders. However, functional variability among MSCs cultured in vitro will lead to distinct therapeutic efficacies. Until now, the mechanisms behind immunomodulatory functional variability in MSCs are still unclear.

Methods

We systemically investigated transcriptomic variations among MSC samples derived from multiple tissues to reveal their effects on immunomodulatory functions of MSCs. We then analyzed transcriptomic changes of MSCs licensed with INFγ to identify potential molecular mechanisms that result in distinct MSC samples with different immunomodulatory potency.

Results

MSCs were clustered into distinct groups showing different functional enrichment according to transcriptomic patterns. Differential expression analysis indicated that different groups of MSCs deploy common regulation networks in response to inflammatory stimulation, while expression variation of genes in the networks could lead to different immunosuppressive capability. These different responsive genes also showed high expression variability among unlicensed MSC samples. Finally, a gene panel was derived from these different responsive genes and was able to regroup unlicensed MSCs with different immunosuppressive potencies.

Conclusion

This study revealed genes with expression variation that contribute to immunomodulatory functional variability of MSCs and provided us a strategy to identify candidate markers for functional variability assessment of MSCs.

Background

Mesenchymal stromal cells (MSCs), also well-known as mesenchymal stem cells, are one of the most comprehensively studied multipotent cells with highly promising applications in immune modulation and regenerative medicine. To date, there are more than 1000 clinical trials registered in ClinicalTrials.gov (http://www.clinicaltrials.gov) attempting to utilize MSCs for cellular therapy [1, 2]. However, the clinical outcomes vary significantly among different cellular therapy trials. For example, in clinical usage of MSCs for treating therapy-refractory graft-versus-host disease (GvHD), only a small proportion of patients had favorable outcomes while many others did not show any apparent efficacy [3]. Functional variation and heterogeneity of MSCs, plus the lack of efficient assays for assessing MSC potency, could be some of the leading reasons resulting in inconsistent clinical outcomes [4].

MSCs have been isolated from various tissues, such as bone marrow [5], adipose tissue [6], umbilical cord [7, 8], and placenta [9, 10]. These cells comply with the minimal criteria defined by International Society for Cellular Therapy (ISCT) in 2006 [11] based on their morphological, phenotypic, and functional characteristics. However, recently, increasing number of MSC-based studies demonstrated that MSCs derived from different donors, tissues, and even sub-clones from the same cell line differed in their functional properties, such as immunomodulatory function, which will impact their applications [12,13,14,15]. Besides cell origins, the heterogeneity of MSCs could also be introduced by different isolation methods, culture media and methods, passage numbers, and/or freezing processes and lead to changes in proliferation and differentiation capacities, as well as in immunosuppression potency [16,17,18].

The unique immunomodulatory plasticity of MSCs makes them an invaluable cell type. MSCs exert their therapeutic effects through forming a balanced inflammatory and regenerative microenvironment, and their immunomodulatory capabilities are not constitutive but rather are licensed by inflammatory cytokines, such as INFγ and TNFα [1]. Licensed MSCs could release various cytokines (such as TGFβ, IL-10, CCL2, IL-6, IL-7, CXCL9, and CXCL10) [19, 20], growth factors (such as HGF and LIF) [21, 22], immunosuppressive molecules (such as NO, PGE2, TSG6, HO1, and galectins) [23, 24], and/or MSC-derived exosomes [25,26,27] to modulate differentiation, maturation, and inflammatory state of immune cells, such as dendritic cells (DCs), macrophages, and monocytes, and promote the formation of regulatory T (Treg) cells and prevent the activation of effector T cells [24]. In addition, MSCs responding to inflammatory environment could also upregulate immunosuppressive molecule PD-L1 which inhibits T cell activation [28] and FASLG which induces T cell apoptosis [29] through cell-to-cell interaction. Recent studies greatly improved our understanding in the immunoregulatory mechanisms of MSCs. However, why different MSC samples differ in immunomodulatory potency before licensing remains unclear and needs to be further investigated [1].

Here, we comprehensively investigated transcriptomic variations among MSC samples derived from multiple tissues. According to the expression patterns, we categorized these samples into 7 groups exhibiting distinct functional enrichment. To identify potential molecular mechanisms that result in distinct MSCs with different immune modulatory potency, we analyzed transcriptome changes of MSCs licensed with INFγ. Differential expression analysis indicated that different groups of MSCs deployed common regulation networks in response to inflammatory stimulation while expression variation of genes in the networks could trigger their different immunosuppressive capability. We also found that these different responsive genes showed high expression variability among unlicensed MSC samples. Finally, a gene panel was derived from these different responsive genes and was able to regroup unlicensed MSCs with different immunosuppressive potency.

Methods

Cell isolation and culture

MSCs were obtained from consenting subjects according to the Institutional Review Board of BGI (BGI-IRB for short) approved guidelines. Human placentas (n = 4) and umbilical cords (n = 7) were collected from naturally delivered full-term newborns. Dental pulps were obtained from donated wisdom teeth (n = 4, age 15–30).

Placenta-derived MSCs (PL-MSC) were isolated from chorionic plate (CP) after removing the amniotic membrane from the fetal surface of the placenta. After washing with PBS, PL were mechanically dissociated into tissue explants of approximately 2 mm2, which were then seeded into T75 flasks and cultured in MSC medium at 37 °C with 5% CO2 in a humidified atmosphere. MSC medium composed of MEMα (Hyclone), 1× L-glutamine (Thermo Fisher Scientific), 10% FBS (Hyclone), and 1× Antibiotic-Antimycotic (Thermo Fisher Scientific).

Umbilical cord-derived MSCs (UC-MSCs) were isolated from Wharton’s jelly (WJ) within the umbilical cord after dissection and removal of the arteries, vein, and amniotic epithelium. Tissue explants were applied to isolate and culture UC-MSCs using the same method as PL after tissue dissection.

Human dental pulp-derived MSCs (DP-MSCs) were isolated from DP uncovered by means of bone forceps to fracture the dental crown in several parts, as described previously [30]. Then, the dental pulp was enzymatically treated with 1 mg/ml type I collagenase (Sigma) and 3 mg/ml type II dispase (Sigma) for 1 h and cultured in MSC medium at 37 °C with 5% CO2 in a humidified atmosphere.

When cell density reached about 80% confluence, cells were dissociated with TrypLE™ Select (ThermoFisher Scientific) incubated at 37 °C for 2 min. The collected cells were passaged about every 3–5 days at a seeding density of 5000 cells/cm2. All assays were performed using MSCs between passages 2 and 5.

IFNγ treatment and cell collection

MSCs were seeded into 6-well plates at a density of 5000 cells/cm2 and then cultured in MSC medium at 37 °C with 5% CO2 in a humidified atmosphere. When cell density reached about 70% confluence, MSCs were stimulated with IFNγ (5 ng/mL, R&D) for 24 h; meantime, parallel untreated wells were used as paired control. After treated for 24 h, the cells were lysed by adding of 1 ml TRIzol (Invitrogen) into each well after removing the medium and washed 3 times with PBS. For total mRNA extraction, each sample was pooled from 2 wells of 6-well plates cultured at the same time.

RNA-Seq library construction and sequencing

Total mRNA was extracted using TRIzol (Invitrogen) reagent, as described previously [31]. Briefly, cells lysed by TRIzol were centrifuged and chloroform was added to the supernatant and mixed well. After spin, supernatant was mixed with chloroform/isopropanol (24:1) and centrifuged again. The same volume of isopropanol was added to the supernatant and stored at − 20 °C for 1 h, and then samples were centrifuged to precipitate RNA. RNA was washed by 75% alcohol twice and dissolved in nuclease-free water. The purity, integrity, and density of RNA were detected by Nanodrop, Agarose gel electrophoresis, and Agilent 2100 BioanaLyzer respectively, and then cDNA was synthesized and PCR was used to construct RNA-Seq library. All protocols for BGISEQ-500 library construction, preparation, sequencing, and quality control were provided by BGI.

To enhance the repeatability of our experiments, 13 cell lines (UC (n = 7), PL (n = 3), AD (n = 2), and PD (n = 1)) banked in the National Institutes for Food and Drug Control were independently cultured and treated with IFNγ with the similar methods as mentioned above. Then, the cells were lysed by TRIzol and shipped to our labs for RNA-seq.

RNA-Seq data processing and quality control

To get public available RNA-seq data of MSCs, we searched in Gene Expression Omnibus using keywords “(MSC OR Mesenchymal stem cell OR Mesenchymal stromal cell) AND “Homo sapiens”[porgn:__txid9606]” with “Expression profiling by high throughput sequencing” type for data information collection. Then, we manually removed samples cultured with certain treatment, from donors with certain disease, with gene modification, or differentiated from ESCs, iPSCs, or other cell types (Fig. S1). We used Illumina platforms for sequencing and retained samples with reads not less than 10,000,000. Totally, we obtained 120 samples, of which raw files were downloaded from NCBI SRA database [32].

Quality control for each sample was performed by FastQC; adaptors and poor quality bases at read ends were trimmed by cutadapt [33] before mapping. Reads were mapped to the human genome (GRCh38) using HISAT2 with default parameters [34]. Raw counts of sequencing reads for the feature of genes were extracted by featureCounts [35]. MSCs RNA-seq data sequenced in our lab was also processed by the same mapping and feature counts extraction methods processes as mentioned above.

After read mapping and raw count extraction, we further compared the percentage of aligned exactly 1 time and median pairwise correlation coefficient for each sample, and we considered samples with percentage of aligned exactly 1 time that are lower than 60% and samples with median pairwise correlation coefficient r less than 0.9 as outliers (Fig. S2a) and left out for the further analysis. Finally, 69 downloaded samples together with our 20 untreated samples were used to investigate transcriptomic variation among MSC samples (Table S1).

Filtering and normalization

Expressed genes in MSCs were defined as genes with counts per million (CPM) value more than 1 in at least 10% of the total samples; others were considered as none (not detected in all samples) or lowly expressed genes (CPM < 1 in at least 90% of the total samples), which were filtered out before normalization. The trimmed mean of M values (TMM) normalization method was used to estimate scale factors between samples and normalize RNA composition by calcNormFactors function in R package edgeR [36, 37].

Variable-expressed gene identification

To quantify variability of gene expression across MSCs, distance-to-median (DM) statistic was used as a corrected version of coefficients of variation (CV), which is independent of the mean expression and gene length, as previously described [38]. Briefly, counts per million (CPM) were computed for the mean-corrected residual of the squared CV of each gene to its corresponding rolling median calculation. To correct for the effect of gene length on the mean corrected residual, DM was defined as the difference between the mean corrected residual of the squared CV of each gene and its expected residual from gene length. We computed the rolling median in 50 windows and set the number of overlapping genes between adjacent windows to 25 [38].

Data dimension reduction, visualization, and clustering

R function prcomp with default parameters was used to perform principal component analysis (PCA) for expression of selected gene sets. R function Rtsne from Rtsne package was applied for visualizing high-dimensional data into a two-dimensional map in Fig. 3a, b with initial_dims = 10, and before running, we set seed to 1. A graph-based clustering approach [39] was used to cluster the samples into different groups. The first 10 PCs in the data were applied to construct an SNN matrix using the FindNeighbors function in Seurat v3 with k.param set to 10. We then identified clusters using the FindClusters command with the resolution parameter set to 2.

Differential expression analysis

To identify differentially expressed genes (DEGs), we used R package edgeR to organize, filter, and normalize the data, and quasi-likelihood F tests were applied to identify DEGs according to the guide [37, 40]. Genes that differed in expression by two folds and with a false discovery rate (FDR) < 0.1 were assigned as DEGs.

GO and pathway enrichment analysis

To find the GO and KEGG terms enriched in defined gene sets, we used the DAVID web-tool [41]. For figures, we only reported the top-ranked terms illustrated in the legends.

Gene set enrichment analysis

Gene set enrichment analysis (GSEA) [42] was performed with 1000 permutations by GSEA_4.0.2 desktop application. Gene lists were ranked by the significant score defined as −log10(FDR) multiplied by log-transformed fold-change between two conditions or DM value for Fig. 2a, g, S3c, and S3d. Gene sets from the Molecular Signature Database (MSigDb) were used for GO and KEGG analysis. Gene sets that contained between 15 and 300 genes were included to provide more biologically meaningful results and reduce false positives.

Pathway enrichment analysis and visualization

Pathway enrichment analysis was achieved according to protocol described in [43]. We downloaded the pathway gene set database Human_GOBP_AllPathways_ no_GO_iea_October_01_2019_symbol.gmt from the Bader lab dated October 01, 2019, for all pathway enrichment analyses. Gene lists were ranked by the significant score defined as −log10(FDR) multiplied by log-transformed fold-change between two conditions. After gene set enrichment analysis (GSEA), pathway sets that contained between 15 and 300 genes were included to provide more biologically meaningful results and reduce false positives.

For map visualization, pathway enrichment analysis results were interpreted in Cytoscape using its EnrichmentMap, AutoAnnotate, WordCloud, and clusterMaker2 applications [43, 44]. The pathway enrichment map was created with parameters FDR Q value < 0.001 and combined coefficient > 0.375 with combined constant = 0.5 [44].

MSC and PBMC co-culture for immunosuppressive potency assessment

The suppressive effect of MSCs on leukocyte expansion was confirmed as described previously [45]. Briefly, MSCs were seeded into 24-well plates at a density of 1 × 105 cells per well, and CFSE (Sigma)-labeled human PBMCs were added to each MSCs well at a 1:5 (cell number) co-culture ratio of MSCs to PBMCs. Then, 10 μg/mL phytohemagglutinin (PHA) (Sigma) was used to activate PBMC cells. PBMCs at a same density without MSCs and PHA were used as negative control. For positive control, we plated the same number of PBMCs/well with PHA to active leukocytes. After 5 days of co-culture, PBMCs were collected and measured using a FACSCalibur platform (BD Biosciences). The suppression of T cell proliferation by MSCs was calculated as [100% − (T cell proliferation after co-culture with MSCs divided by positive control × 100)%]. Negative control was applied to define a threshold of the CFSE signal of non-proliferating T cells.

Results

Data selection and quality control

To comprehensively investigate transcriptomic variations among MSC cell lines cultured in vitro, RNA-seq data of total 102 samples were integrated for gene expression analysis, of which 69 samples were selected from public database and 33 samples were newly sequenced in this study (Fig. S1 and S2; Table S1). Overall, the MSC samples were derived from 6 tissues (Fig. 1a) in 17 studies (Fig. S2b), including adipose tissue (AD), bone marrow (BM), dental pulp (DP), endometrial (ED), placenta (PL), and umbilical cord (UC). Of these samples, number of reads were mostly between 10,000,000 and 60,000,000, and more than 60% aligned exactly 1 time to the transcriptome (Fig. S2c), indicating high quality of these collected RNA-seq data.

Fig. 1
figure 1

Overview of data collected for transcriptome variation analysis. a Barplot showing number of samples derived from different tissue origins used for transcriptome variation analysis. b, c Boxplots showing distribution of expression (up) and ranking (down) for positive markers (b) and negative markers (c) across samples, respectively. Some MSC samples, including AM CM, CP, and CV, derived from different anatomical parts of placenta. AD, adipose tissue; AM, amniotic membrane; BM, bone marrow; CM, chorionic membrane; CP, chorionic plate; CV, chorionic villi; ED, endometrial; DP, dental pulp; PL, placenta; UC, umbilical cord

Full size image

Minimal criteria of defining MSC claimed that MSCs must express three positive markers, i.e., CD105 (ENG), CD73 (NT5E), and CD90 (THY1), and lack expression of several negative markers, including CD45 (PTPRC), CD34, CD14 or CD11b (ITGAM), CD79a (CD79A) or CD19, and HLA-DR (HLA-DRA and HLA-DRB1 etc.) [11]. Indeed, gene expression level ranked by TPM (transcripts per kilobase million) showed that those positive markers were highly expressed (Fig. 1b) while negative markers were weakly or not expressed in our samples, except HLA-DR molecules showed highly variable expression across samples (Fig. 1c). Considering that MSCs express HLA-DR surface molecules not only in response to stimulations, such as IFNγ, but also under some normal expansion culture conditions [46, 47], therefore, we did not remove these samples with higher expression of HLA-DR for the following gene expression analysis.

Expression variations and characteristics among MSC samples

To analyze transcriptomic variations across MSC samples, the DM value was used to rank variations of genes expressed in MSCs. A larger DM value indicated a greater variation (Fig. S3a and S3b; Table S2). Interestingly, we noticed that growth signaling genes such as BMP2/4/6, WNT2/4/11, NOG, and TGFB2, as well as transcript factors such as RUNX3, GATA2, SOX4, SOX11, HES1, EGR2, FOS, and FOSB, which are involved in differentiation process of MSCs [48,49,50,51,52,53,54], were among the top genes with large DM values (Fig. 2a, b). Meanwhile, genes that regulated immunomodulation [1, 19, 55], such as CD274 (PD-L1), TNFAIP6 (TSG6), CCL2, CCL5, IL6, CSF3, HGF, and ICAM1, also exhibited highly variable expression among those samples (Fig. 2c). In contrast, genes involved in several other biological processes, such as metabolic process, gene expression, RNA process, and RNA binding, had very small gene expression changes (Fig. 2d). Therefore, we hypothesized that gene expression variations among MSCs may be used to select MSCs with different immunomodulation potency and differentiation propensity.

Fig. 2
figure 2

Transcriptomic variations across MSC samples. a–d Expression distribution of certain signaling-pathway related genes (a), transcription factors (b), and immune-related genes (c) with higher DM value compared to metabolic process-related genes (d). e, f GESA positive results showing enrichment in GO biological process (e) and KEGG pathway (f) gene sets database based on ranked genes list in descending order by the DM value. Only the top 20 terms with highest NES were presented (p < 0.001)

Full size image

To further investigate whether expression variations among the MSC samples were related to specific MSC biological functional properties or not, we next performed gene set enrichment analysis (GSEA) with above pre-ranked gene list. Notably, genes with highly variable expression were significantly enriched in immune modulation and developmental process, such as humoral immune response, response to chemokine, nephron development, cardiac chamber morphogenesis, and digestive system development (Fig. 2e, f). Expression variability was also overrepresented in coding proteins located in extracellular matrix and cell surface, involving in cytokine, receptor regulator, and binding activity (Fig. S3c). On the other hand, genes with stable expression were housekeeping genes involved in basic cellular function, such as RNA processing (Fig. S3d). These results together demonstrated that transcriptome-wide sample-to-sample variations among MSCs were associated with their functional properties.

Grouping MSC samples based on highly variable genes

To identify candidate groups of the MSC samples based on gene expression pattern, genes with DM value more than 1 as highly variable genes (HVGs) were selected for data dimension reduction and clustering. The results showed that the samples we collected could be clustered into 7 groups (G0–G6) (Fig. 3a). Among these groups, G1 included samples mostly derived from BM while G3 included all AD-MSCs plus some BM-MSCs (Fig. 3b; Table S3). Other groups included MSC samples derived from different origins as well (Fig. 3b; Table S3). Meanwhile, the HVGs we selected grouped into five clusters with distinct functional enrichment, such as system development, tube development, metabolic process, and response to cytokine (Fig. 3b), indicating that different groups of MSC samples may have different differentiation propensity and immunomodulatory potency. In addition, different expression analysis among these groups demonstrated that functional enrichment of DEGs among these groups are associated with MSC function-related properties as well, such as angiogenesis, nervous system development, cell migration, and inflammatory response. Taken together, these results demonstrated that MSCs from different tissue origins could be classified into the same group with similar functional enrichment based on expression of HVGs, although tissue origins have been reported to impact greatly on functions of MSCs [56, 57].

Fig. 3
figure 3

Groups and functional enrichment in MSCs. a, b tSNE visualizing the results of samples clustering (a) and tissue origins (b). c Heatmap showing expressive characteristics of HVGs among the groups. d Results of GO biological process enrichment analysis for each genes cluster. Top 5 terms with lowest p-adjusted value were presented

Full size image

To illustrate expression difference in these groups, we presented results obtained from comparison among G0 (n = 22, 15 downloaded, 7 newly sequenced), G2 (n = 15, 1 downloaded, 15 newly sequenced), G3 (n = 13, 11 downloaded, 2 newly sequenced), and G4 groups (n = 12, 3 downloaded, 9 newly sequenced) (Table S4), to which our 33 MSC samples were assigned (Table S3). GO biological process enrichment analysis demonstrated that G0 significantly upregulated genes involved in response to stimulus and inflammatory response, including cytokines such as CXCL2, CXCL3, CXCL5, and CXCL20 (Fig. S4A and B). The upregulated genes in G2, G3, and G4 were overrepresented in developmental process with distinct developmental cell types. For example, genes related to nervous system development, circulatory system development, and nephron development were respectively upregulated in the G2, G3, and G4 groups (Fig. S4C-S4H). Altogether, we found significant gene expression variations existing among MSCs that could potentially influence their functional properties. Therefore, we hypothesized that quantitative RNA analysis of selected genes from HVGs could serve as a candidate matrix assay for characterizing MSC potency [58].

Characteristics of expression changes in MSCs upon IFNγ licensing

Although above functional enrichment analysis demonstrated that some inflammatory response-related genes were over-represented in G0 (Fig. S4), it was not clear how these differences would affect MSC immunomodulatory behavior. Due to the critical role of IFNγ in licensing MSC-mediated T cell suppression [28], IFNγ, which can be used as an alternative for human PBMCs as responder cells in a MSC potency assay [19], was used to treat 27 MSCs within the G0 (n = 4), G2 (n = 15), and G4 (n = 8) groups to study response variation on transcriptomic level. For clarity, we assigned them as ssG0, ssG2, and ssG4 due to their small size compared to the number of samples in each group.

Comparing to their paired untreated samples, we identified 902 genes that were differentially expressed (655: 72.62% upregulated vs 247: 31.26% downregulated) (Fig. 4a, Table S5). In line with previous studies, IDO1, the dominant determinant of MSC-mediated inhibition of T cell proliferation, and chemokines, including CCL5, CXCL9, CXCL10, and CXCL11, were upregulated in MSCs upon IFNγ licensing, which could potentially form a chemokine-IDO axis to exert immunoregulatory effects on various immune cells (Fig. 4b) [19, 59,60,61]. In addition, cytokines, including CCL2, CCL7, and IL6, apoptosis inducer TNFSF10 (TRAIL), immune checkpoint proteins, including CD274 (PD-L1) and PDCD1LG2 (PD-L2), cell adhesion molecules, including ICAM1 and VCAM1, and class II major histocompatibility complex (MHC), including HLA-DRA, HLA-DRB1, and HLA-DRB5, were also overexpressed in IFNγ activated MSCs (Fig. 4b and Table S5). Similarly, these genes have also been shown to play critical roles in MSC-mediated immunosuppression in previous studies [58]. Meanwhile, IFNγ licensing triggered specific signaling pathways in MSCs, such as upregulation of JAK2, JAK3, STAT1, STAT2, SOSC1, SOSC3, TLR2, and TLR3, to orchestrate their immune response. Pathway enrichment map further illustrated that IFNγ-licensed MSCs upregulated several gene clusters linked to response to interferon, immune response, and antigen and protein degradation (Fig. 4c). Interestingly, sterol biosynthetic process was downregulated in IFNγ-treated MSCs (Fig. 4c). Sterols are the major component of the cellular membranes and are essential for mammalian cell growth [62]. Decreased sterol synthesis could partially explain why IFNγ leads to a cytostatic response in MSCs [60]. Overall, a panel of immunomodulation-related genes was upregulated upon IFNγ licensing, and here we define these upregulated genes as common response genes (CRGs).

Fig. 4
figure 4

Gene expression changes in MSCs treated with IFNγ. a A mean difference plot showing DEGs identified in MSCs treated with IFNγ. b Representative genes upregulated in IFNγ-licensed MSCs. c Enrichment map showing pathways enriched in INFγ-treated MSCs (red) and untreated MSCs (blue). d Boxplot showing immunosuppressive scores for each group calculated based on VEGF, IFNa, CXCL10, GCSF, CXCL9, IL-7, and CCL2 expression level. e Results showing MSCs from ssG0, ssG2, and ssG4 co-cultured with PBMC for immunosuppressive potency assessment

Full size image

To predict T cell suppression potency, the sum of log normalized expression of VEGF, IFNa, CXCL10, GCSF, CXCL9, IL-7, and CCL2 genes, which were correlated with T cell suppression capacity of MSCs according to previous studies [19, 63], was calculated and served as MSC immunosuppressive score for the 27 MSCs within ssG0, ssG2, and ssG4 treated with IFNγ. Our results showed that MSC immunosuppressive score in ssG0 was significantly lower than those of ssG2 and ssG4 (Fig. 4d). MSC and PBMC co-culture experiments in vitro were performed on 16 MSCs (4 from ssG0, 8 from ssG2, and 4 from ssG4), and the results demonstrated that T cell proliferation inhibitory rate of G0 was significantly lower than those of ssG2 and ssG4. Taken together, these results demonstrated that different groups of MSCs clustered by their expression patterns of HVGs across unlicensed MSC samples could have distinct immunosuppressive capability, which may reflect on their different responses to inflammatory environment.

Transcriptional variations of MSCs in response to inflammatory environment imitated by IFNγ treatment

To identify genes that respond differentially to IFNγ licensing among the ssG0, ssG2, and ssG4 groups, we performed differentially expressed analysis and obtained a total of 472 genes defined as different response genes (DRGs). Most of the DRGs downregulated in ssG0 while upregulated in ssG2 and ssG4 are enriched in immune response pathways (Fig. S5A-S5D), including several well-known immune-modulating genes, such as CXCL9, CXCL10, CXCL11, CCL2, CCL7, CCL8, CD74, CXCL16, CD7, CD14, CD83, and LGALS9 (Fig. 5a; Table S6). These genes are tightly involved in immunomodulatory processes, such as regulation of immune cell migration, T cell development and differentiation, T cell chemotaxis, immune activation, and cell survival [64,65,66,67].

Fig. 5
figure 5

Expression variability of MSCs in response to IFNγ. a Representative genes showing different responsive to IFNγ licensing. b Venn diagram showing number of genes overlapped between CRGs and DRGs. Genes set1, Genes set2, and Genes set3 represented unique genes in CRGs, shared genes, and unique genes in DRGs respectively. c Venn diagram showing number of genes overlapped between genes downregulated in MSCs treated with IFNγ and DRGs. d Results of GO BP enrichment analysis for the Genes set1, Genes set2, and Genes set3

Full size image

To further identify the genes related to common immune-regulatory pathways and to differences in immunomodulatory capacities, we performed comparison between DRGs and CRGs. The genes unique to CRGs, shared between CRGs and DRGs, and unique to DGRs were designated as Genes set1, Genes set2, and Genes set3, respectively (Fig. 5b). Among the 472 DRGs, 205 genes were shared with CRGs and fall within Gene set2 (Fig. 5b and c; Table S7). Interestingly, functional enrichment analysis revealed that Genes set1 and Genes set2 were significantly involved in immunomodulatory functions, while Genes set3 was involved in regulation of developmental process (Fig. 5d), implying that IFNγ-treated MSCs may partly influence on their developmental behaviors. Taken together, these results suggested that MSCs exert immunosuppressive effects through common mechanisms, while variations in expression of some immunomodulatory genes upon inflammatory priming could result in their distinct immunosuppressive potency.

Refining gene panel for grouping unlicensed MSCs with predictable immunosuppressive potency

According to our results, differential expression of certain genes under IFNγ-licensed state could potentially explain for differences in their immunomodulatory activity. However, there is still a lack of reports on genes whose expression in MSCs under unlicensed condition will be related to MSCs’ immunomodulatory potency.

Considering that genes in Genes set2 were not only related to common immune-regulatory pathways but also to differences in immunomodulatory capacities (Fig. 5b, c), differential expression of these genes under unlicensed condition may contribute to distinct immunomodulatory behavior of MSCs in response to inflammatory environment. Interestingly, when we compared DM distribution of genes in Genes set1, Genes set2, and Genes set3 in the total unlicensed MSC samples, genes in Genes set2 demonstrated significantly higher variation than Genes Set1 (p = 2.06e−11) and Genes Set3 (p = 1.80e−02) (Fig. 6a). Several immune response-related genes, such as CCL2, CCL7, CD74, TNFSF10, LGALS9, IFIT1, VCAM1, and ICAM1, fall within Genes set2 and were among the top highly variable genes (Fig. 2c and Table S2). These results indicated that expression variation of genes in Genes set2 may exert greater influence during immune activation of unlicensed MSCs.

Fig. 6
figure 6

Candidate markers for MSC immunomodulatory potency assay. a Violin plot showing DM values distribution for the Genes set1, Genes set2, and Genes set3. p value were presented for one set compared to another performed by Wilcoxon rank sum test. b Scatter plot visualizing PCA results using top 100 genes with the highest DM values in the Genes set2 expressed in the untreated samples

Full size image

Then, we applied the top 100 genes with highest DM values in Genes set2 as a gene panel and utilized their expression for data dimension reduction. The results of principal component analysis (PCA) showed that all MSC samples in the G0 group laid in the third quadrant while majority of the samples from the G1 to G5 groups laid in the first and second quadrant based on principal component 1 and principal component 2 (Fig. 6b). Since the G0 group exhibited the lowest immunosuppressive capacity compared to the G2 and G4 groups in both in silico and in vitro analysis (Fig. 5e), samples fall within the same quadrant as G0 may have lower immunosuppressive capability as well. Meanwhile, with INFγ-treated samples included, PCA analysis showed that the G6 group, of which all samples laid in the fourth quadrant, were closer to INFγ-treated samples (Fig. S6), indicating pre-licensed state of these samples. Taken together, our results demonstrated that the gene panel we selected here would be valuable for characterizing MSC immunomodulatory potency based on quantitative RNA-seq analysis.

Discussion

Since the first report on the characteristics of MSCs derived from human BM [68], studies revealed that MSCs can be isolated from both pericytes and adventitial progenitor cells from nearly all tissues [69,70,71]. Although MSCs have been widely accepted as one of the most promising cell products to treat various degenerative and inflammatory disorders, such as graft versus host disease (GvHD), Crohn’s disease, multiple sclerosis, and diabetes, there are still clinical challenges, such as why the outcomes of advanced clinical trials were not as encouraging as pre-clinical animal data in a wide array of disease models [17]. In addition to limited understanding of mechanisms of action which MSCs deployed to regulate their anti-inflammatory and tissue repair functionalities, functional variability and heterogeneity could hinder development of effective assays for MSCs as potency release criterion for the advanced clinical trials [1, 17, 58]. These variability and heterogeneity manifest among donors, among tissue sources, as well as within cell populations [72,73,74,75]. Besides, distinct cell separation and preservation methods, culture media, and number of passages can affect cell functionality. For example, human umbilical cord blood mononuclear cells tested before and after cryopreservation showed different abilities to treat stroke [76], aged MSCs underwent morphological, phenotypic, and differentiation potential changes [77], and long-term culture increased genetic instability in MSCs [78]. These studies indicated that MSCs with distinct cell preparation, fitness, culture methods, and expansion levels could differ in their tissue-protective and immunomodulatory properties [17, 19]. However, the molecular contributors to the functional variability and heterogeneity remain unclear. Here, we analyzed RNA-seq from 102 MSC samples derived from 6 different tissues, including AD, BM, ED, PD, PL, and UC, to investigate their expression variations and relationship with their functional variability and identify candidate markers differently expressed influencing on their immunosuppressive potency.

Several studies have been done to compare gene expression similarity and variability among MSC samples [79,80,81,82]. These studies demonstrated that different gene expression profiles could reflect the ontogenetic sources of MSCs and indicate their distinct differentiation potential or other functional properties. In line with these studies, our results also showed that MSCs were mostly grouped together with the same tissue origin according to the expression pattern of HVGs (Fig. 3a, b; Table S3). However, these researches largely focused on expression difference among MSCs with distinct tissue origins while ignoring that MSCs from the same tissue might also exhibit functional variability. The functional differences could come from a variety of cues, including chemical, physical, and biological factors, expansion level, and characters of donors, which may result in changes of MSC functional characteristics [80, 83]. For example, compared to younger counterparts, aged MSCs from the same tissue with identical culture method display delayed clonogenic capacity and pro-inflammatory SASP-like phenotype, and their immunomodulatory properties were significantly reduced [83]. To address this, we performed data dimension reduction and clustering based on a nonparametric clustering technique (see the “Methods” section) to group these collected samples in the present study. Our results demonstrated that MSCs can be clustered into groups with diverse functional properties characterized by enrichment analysis (Fig. 3 and S4). Besides, MSCs from different tissues can be classified into the same group while MSCs from the same tissue as well can be clustered into different groups according to the expression patterns of HVGs (Fig. 3a, b; Table S3), indicating the importance of potency assays for MSCs before clinical trials or application.

Despite different tissue sourcing, our results are in line with that MSCs likely share fundamental mechanisms of action mediating their anti-inflammatory processes [58]. Our data are in agreement with reports that immunosuppression related molecules, such as IDO1, CCL5, CXCL9, CXCL10, CXCL11, CD274, TNSF10, CCL2, and FLT3LG (FLT3L), were upregulated upon INFγ licensing (Fig. 4b; Table S6) [60, 61], of which some were lowly or not expressed in unlicensed MSCs. Activated MSCs produce chemokines CCL5, CXCL9, CXCL10, and CXCL11, which could recruit T cells to the proximity of MSCs and suppress the proliferation and activity of T cells in their vicinity by expressing tryptophan catabolism rate-limiting enzyme IDO1 through metabolite kynurenic acid and/or by expression the immune checkpoint protein CD274 through cell-to-cell interaction [1, 84, 85]. In addition, a recent study demonstrated that MSCs might utilize IFNγ-FLT3L-FLT3 axis to suppress inflammation in lupus through upregulating tolerogenic DCs [86]. Different groups of MSCs should deploy shared regulation networks to exert immunosuppressive function upon IFNγ licensing, including JAT-STAT, NF-kappaB, IL-12/IL23, response to interferon, immune response, antigen and protein degradation, extrinsic apoptosis pathway, and complement system signaling pathways, which could form a regulatory network to orchestrate MSC immunomodulatory function (Fig. 4c). However, gene expression variations could result in different responses among MSC groups treated with INFγ, and these immunomodulation related genes were mainly in the shared regulation networks, including abovementioned TNSF10, CXCL9, CXCL10, CXCL11, and CCL2. Furthermore, human MSCs licensed by INFγ have been tested in NOD-SCID mice showing enhanced immunosuppressive properties to significantly reduce the symptoms of GvHD [61], indicating potential clinical application of INFγ primed MSCs. Nevertheless, expression variability among MSCs, which could lead to different expression levels of immunomodulatory genes among the licensed MSCs, implied their different immunomodulatory potency after priming. Therefore, conditions for MSC priming, such as optimum priming time, the concentration of INFγ, and whether different MSCs could be adjusted to similar immunomodulatory potency by priming design, need to be refined. And experimental immunosuppression and immunomodulation strategies could also be applied to enhance the predictive value of preclinical studies with MSCs [87].

In essential, genetic and epigenetic variations contribute to functional variability among MSCs. Identification of functional markers of potency in unlicensed MSCs could facilitate our understanding of MSCs’ mechanisms of action and development of release potency assays for them as potency release criterion [58]. IFNγ stimulation of MSCs recapitulates the molecular genetic changes that are observed in MSCs co-cultured with activated PBMCs [19]; however, it is necessary to notice that therapeutic effects of MSCs are multifaceted synergistic responses, which could form a balanced inflammatory and regenerative micro-environment in the presence of vigorous inflammation [1], and assays just focused on one functional aspect of MSCs, such as immunosuppressive potency, may ignore other functional capabilities of MSCs, which may also in some extend link to clinical results. Our analysis demonstrated that transcriptome-wide sample-to-sample variations among MSCs are associated with various functional properties (Fig. 2). Besides, their functional similarity and disparity can be classified based on expression of HVGs (Fig. 3), and thus, we speculated that these HVGs are valuable to serve as candidate matrix assays for potency analysis of MSCs without licensing. Comparison of response patterns to IFNγ among MSCs further showed that genes shared between CRGs and DRGs are significantly more variable than the other two sets (Fig. 6). Based on these genes’ expression, we established a primary model, which faithfully assessed immunosuppressive potency of unlicensed MSCs (Fig. 6). Beyond this, we inferred that RNA-seq technology combined with our model method can be extended to other functional variations of MSCs, such as interaction with innate immune cells and differentiation propensity.

Conclusions

In summary, our study demonstrated that MSC samples can be classified into groups exhibiting distinct functional properties, such as immune modulatory potency, according to the expression pattern of HVGs. We also highlighted that MSCs deployed common regulation networks to exert immunosuppressive function while expression variability of genes in the networks could result in distinct immunosuppressive potency in MSCs. Finally, we found these different responsive genes showed high expression variability among unlicensed MSC samples as well, from which candidate markers were refined for development of matrix assays to quantify the immunosuppression potency of human unlicensed MSCs. In the future, with increased number of MSC samples, our analysis approach can be extended beyond immune modulatory potency to characterize other functional variations and related genes.

Availability of data and materials

The data that support the findings of this study have been deposited in the CNSA (https://db.cngb.org/cnsa/) of CNGBdb with accession number CNP0000966. All data analyzed during this study have been included in this published article and its supplementary information files.

Abbreviations

AD:

Adipose tissue

BM:

Bone marrow

DCs:

Dendritic cells

DEGs:

Differentially expressed genes

DM:

Distance-to-median

DP:

Dental pulp

DRGs:

Different response genes

ED:

Endometrial

GSEA:

Gene set enrichment analysis

GvHD:

Graft-versus-host disease

INFγ:

Interferon gamma

ISCT:

International Society for Cellular Therapy

MSCs:

Mesenchymal stromal cells

NO:

Nitrogen oxide

NOG:

Noggin

NT5E:

5′-Nucleotidase

PBMCs:

Peripheral blood mononuclear cells

TPM:

Transcripts per kilobase million

tSNE:

t-distributed stochastic neighbor embedding

UC:

Umbilical cord

References

  1. Shi Y, Wang Y, Li Q, Liu K, Hou J, Shao C, Wang Y. Immunoregulatory mechanisms of mesenchymal stem and stromal cells in inflammatory diseases. Nat Rev Nephrol. 2018;14(8):493–507.

    Article  CAS  PubMed  Google Scholar 

  2. Moll G, Ankrum JA, Kamhieh-Milz J, Bieback K, Ringden O, Volk HD, Geissler S, Reinke P. Intravascular mesenchymal stromal/stem cell therapy product diversification: time for new clinical guidelines. Trends Mol Med. 2019;25(2):149–63.

    Article  PubMed  Google Scholar 

  3. Elgaz S, Kuci Z, Kuci S, Bonig H, Bader P. Clinical use of mesenchymal stromal cells in the treatment of acute graft-versus-host disease. Transfusion Med Hemotherapy. 2019;46(1):27–34.

    Article  Google Scholar 

  4. Wilson A, Hodgson-Garms M, Frith JE, Genever P. Multiplicity of mesenchymal stromal cells: finding the right route to therapy. Front Immunol. 2019;10:1112.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  5. Prockop DJ. Marrow stromal cells as stem cells for nonhematopoietic tissues. Science. 1997;276(5309):71–4.

    Article  CAS  PubMed  Google Scholar 

  6. Zuk PA, Zhu M, Ashjian P, De Ugarte DA, Huang JI, Mizuno H, Alfonso ZC, Fraser JK, Benhaim P, Hedrick MH. Human adipose tissue is a source of multipotent stem cells. Mol Biol Cell. 2002;13(12):4279–95.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  7. Romanov YA, Svintsitskaya VA, Smirnov VN. Searching for alternative sources of postnatal human mesenchymal stem cells: candidate MSC-like cells from umbilical cord. Stem Cells. 2003;21(1):105–10.

    Article  PubMed  Google Scholar 

  8. Wang HS, Hung SC, Peng ST, Huang CC, Wei HM, Guo YJ, Fu YS, Lai MC, Chen CC. Mesenchymal stem cells in the Wharton’s jelly of the human umbilical cord. Stem Cells. 2004;22(7):1330–7.

    Article  PubMed  Google Scholar 

  9. Fukuchi Y, Nakajima H, Sugiyama D, Hirose I, Kitamura T, Tsuji K. Human placenta-derived cells have mesenchymal stem/progenitor cell potential. Stem Cells. 2004;22(5):649–58.

    Article  CAS  PubMed  Google Scholar 

  10. In ‘t Anker PS, Scherjon SA, Kleijburg-van der Keur C, de Groot-Swings GM, Claas FH, Fibbe WE, Kanhai HH: Isolation of mesenchymal stem cells of fetal or maternal origin from human placenta. Stem Cells 2004, 22(7):1338–1345.

  11. Dominici M, Le Blanc K, Mueller I, Slaper-Cortenbach I, Marini F, Krause D, Deans R, Keating A, Prockop D, Horwitz E. Minimal criteria for defining multipotent mesenchymal stromal cells. The International Society for Cellular Therapy position statement. Cytotherapy. 2006;8(4):315–7.

    Article  CAS  PubMed  Google Scholar 

  12. Jin HJ, Bae YK, Kim M, Kwon SJ, Jeon HB, Choi SJ, Kim SW, Yang YS, Oh W, Chang JW. Comparative analysis of human mesenchymal stem cells from bone marrow, adipose tissue, and umbilical cord blood as sources of cell therapy. Int J Mol Sci. 2013;14(9):17986–8001.

    Article  PubMed  PubMed Central  Google Scholar 

  13. Ma J, Wu J, Han L, Jiang X, Yan L, Hao J, Wang H. Comparative analysis of mesenchymal stem cells derived from amniotic membrane, umbilical cord, and chorionic plate under serum-free condition. Stem Cell Res Ther. 2019;10(1):19.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  14. Wang L, Tran I, Seshareddy K, Weiss ML, Detamore MS. A comparison of human bone marrow-derived mesenchymal stem cells and human umbilical cord-derived mesenchymal stromal cells for cartilage tissue engineering. Tissue Eng A. 2009;15(8):2259–66.

    Article  CAS  Google Scholar 

  15. Yoo KH, Jang IK, Lee MW, Kim HE, Yang MS, Eom Y, Lee JE, Kim YJ, Yang SK, Jung HL, et al. Comparison of immunomodulatory properties of mesenchymal stem cells derived from adult human tissues. Cell Immunol. 2009;259(2):150–6.

    Article  CAS  PubMed  Google Scholar 

  16. Corotchi MC, Popa MA, Remes A, Sima LE, Gussi I, Lupu Plesu M. Isolation method and xeno-free culture conditions influence multipotent differentiation capacity of human Wharton’s jelly-derived mesenchymal stem cells. Stem Cell Res Ther. 2013;4(4):81.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  17. Galipeau J, Sensebe L. Mesenchymal stromal cells: clinical challenges and therapeutic opportunities. Cell Stem Cell. 2018;22(6):824–33.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  18. Pollock K, Sumstad D, Kadidlo D, McKenna DH, Hubel A. Clinical mesenchymal stromal cell products undergo functional changes in response to freezing. Cytotherapy. 2015;17(1):38–45.

    Article  PubMed  Google Scholar 

  19. Chinnadurai R, Rajan D, Qayed M, Arafat D, Garcia M, Liu Y, Kugathasan S, Anderson LJ, Gibson G, Galipeau J. Potency analysis of mesenchymal stromal cells using a combinatorial assay matrix approach. Cell Rep. 2018;22(9):2504–17.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  20. Maffioli E, Nonnis S, Angioni R, Santagata F, Cali B, Zanotti L, Negri A, Viola A, Tedeschi G. Proteomic analysis of the secretome of human bone marrow-derived mesenchymal stem cells primed by pro-inflammatory cytokines. J Proteome. 2017;166:115–26.

    Article  CAS  Google Scholar 

  21. Benkhoucha M, Molnarfi N, Dunand-Sauthier I, Merkler D, Schneiter G, Bruscoli S, Riccardi C, Tabata Y, Funakoshi H, Nakamura T, et al. Hepatocyte growth factor limits autoimmune neuroinflammation via glucocorticoid-induced leucine zipper expression in dendritic cells. J Immunol. 2014;193(6):2743–52.

    Article  CAS  PubMed  Google Scholar 

  22. Cao W, Yang Y, Wang Z, Liu A, Fang L, Wu F, Hong J, Shi Y, Leung S, Dong C, et al. Leukemia inhibitory factor inhibits T helper 17 cell differentiation and confers treatment effects of neural progenitor cell therapy in autoimmune disease. Immunity. 2011;35(2):273–84.

    Article  CAS  PubMed  Google Scholar 

  23. Dyer DP, Salanga CL, Johns SC, Valdambrini E, Fuster MM, Milner CM, Day AJ, Handel TM. The anti-inflammatory protein TSG-6 regulates chemokine function by inhibiting chemokine/glycosaminoglycan interactions. J Biol Chem. 2016;291(24):12627–40.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  24. Ma S, Xie N, Li W, Yuan B, Shi Y, Wang Y. Immunobiology of mesenchymal stem cells. Cell Death Differ. 2014;21(2):216–25.

    Article  CAS  PubMed  Google Scholar 

  25. Li T, Yan Y, Wang B, Qian H, Zhang X, Shen L, Wang M, Zhou Y, Zhu W, Li W, et al. Exosomes derived from human umbilical cord mesenchymal stem cells alleviate liver fibrosis. Stem Cells Dev. 2013;22(6):845–54.

    Article  CAS  PubMed  Google Scholar 

  26. Nargesi AA, Lerman LO, Eirin A. Mesenchymal stem cell-derived extracellular vesicles for renal repair. Current gene therapy. 2017;17(1):29–42.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  27. Morrison TJ, Jackson MV, Cunningham EK, Kissenpfennig A, McAuley DF, O’Kane CM, Krasnodembskaya AD. Mesenchymal stromal cells modulate macrophages in clinically relevant lung injury models by extracellular vesicle mitochondrial transfer. Am J Respir Crit Care Med. 2017;196(10):1275–86.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  28. Sheng H, Wang Y, Jin Y, Zhang Q, Zhang Y, Wang L, Shen B, Yin S, Liu W, Cui L, et al. A critical role of IFNgamma in priming MSC-mediated suppression of T cell proliferation through up-regulation of B7-H1. Cell Res. 2008;18(8):846–57.

    Article  CAS  PubMed  Google Scholar 

  29. Akiyama K, Chen C, Wang D, Xu X, Qu C, Yamaza T, Cai T, Chen W, Sun L, Shi S. Mesenchymal-stem-cell-induced immunoregulation involves FAS-ligand-/FAS-mediated T cell apoptosis. Cell Stem Cell. 2012;10(5):544–55.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  30. Pierdomenico L, Bonsi L, Calvitti M, Rondelli D, Arpinati M, Chirumbolo G, Becchetti E, Marchionni C, Alviano F, Fossati V, et al. Multipotent mesenchymal stem cells with immunosuppressive activity can be easily isolated from dental pulp. Transplantation. 2005;80(6):836–42.

    Article  PubMed  Google Scholar 

  31. Sun C, Zhang J, Zheng D, Wang J, Yang H, Zhang X. Transcriptome variations among human embryonic stem cell lines are associated with their differentiation propensity. PLoS One. 2018;13(2):e0192625.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  32. Leinonen R, Sugawara H, Shumway M. International Nucleotide Sequence Database C: the sequence read archive. Nucleic Acids Res. 2011;39(Database issue):D19–21.

    Article  CAS  PubMed  Google Scholar 

  33. Martin M: Cutadapt removes adapter sequences from high-throughput sequencing reads. 2011 2011, 17(1):3.

  34. Pertea M, Kim D, Pertea GM, Leek JT, Salzberg SL. Transcript-level expression analysis of RNA-seq experiments with HISAT, StringTie and Ballgown. Nat Protoc. 2016;11(9):1650–67.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  35. Liao Y, Smyth GK, Shi W: featureCounts: an efficient general purpose program for assigning sequence reads to genomic features. Bioinformatics 2014, 30(7):923–930.

  36. Robinson MD, Oshlack A. A scaling normalization method for differential expression analysis of RNA-seq data. Genome Biol. 2010;11(3):R25.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  37. Robinson MD, McCarthy DJ, Smyth GK. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 2010;26(1):139–40.

    Article  CAS  PubMed  Google Scholar 

  38. Kolodziejczyk AA, Kim JK, Tsang JC, Ilicic T, Henriksson J, Natarajan KN, Tuck AC, Gao X, Buhler M, Liu P, et al. Single cell RNA-sequencing of pluripotent states unlocks modular transcriptional variation. Cell Stem Cell. 2015;17(4):471–85.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  39. Macosko EZ, Basu A, Satija R, Nemesh J, Shekhar K, Goldman M, Tirosh I, Bialas AR, Kamitaki N, Martersteck EM, et al. Highly parallel genome-wide expression profiling of individual cells using nanoliter droplets. Cell. 2015;161(5):1202–14.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  40. Lun AT, Chen Y, Smyth GK. It’s DE-licious: a recipe for differential expression analyses of RNA-seq experiments using quasi-likelihood methods in edgeR. Methods Mol Biol. 2016;1418:391–416.

    Article  PubMed  Google Scholar 

  41. Huang da W, Sherman BT, Lempicki RA: Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat Protoc 2009, 4(1):44–57.

  42. Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, Paulovich A, Pomeroy SL, Golub TR, Lander ES, et al. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci U S A. 2005;102(43):15545–50.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  43. Reimand J, Isserlin R, Voisin V, Kucera M, Tannus-Lopes C, Rostamianfar A, Wadi L, Meyer M, Wong J, Xu C, et al. Pathway enrichment analysis and visualization of omics data using g: Profiler, GSEA, Cytoscape and EnrichmentMap. Nat Protoc. 2019;14(2):482–517.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  44. Merico D, Isserlin R, Stueker O, Emili A, Bader GD. Enrichment map: a network-based method for gene-set enrichment visualization and interpretation. PLoS One. 2010;5(11):e13984.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  45. Hsu PJ, Liu KJ, Chao YY, Sytwu HK, Yen BL. Assessment of the immunomodulatory properties of human mesenchymal stem cells (MSCs). J Visualized Experiments. 2015;106:e53265.

    Google Scholar 

  46. Grau-Vorster M, Rodriguez L, Torrents-Zapata S, Vivas D, Codinach M, Blanco M, Oliver-Vila I, Garcia-Lopez J, Vives J. Levels of IL-17F and IL-33 correlate with HLA-DR activation in clinical-grade human bone marrow-derived multipotent mesenchymal stromal cell expansion cultures. Cytotherapy. 2019;21(1):32–40.

    Article  CAS  PubMed  Google Scholar 

  47. Grau-Vorster M, Laitinen A, Nystedt J, Vives J. HLA-DR expression in clinical-grade bone marrow-derived multipotent mesenchymal stromal cells: a two-site study. Stem Cell Res Ther. 2019;10(1):164.

    Article  PubMed  PubMed Central  Google Scholar 

  48. van de Peppel J, Strini T, Tilburg J, Westerhoff H, van Wijnen AJ, van Leeuwen JP. Identification of three early phases of cell-fate determination during osteogenic and adipogenic differentiation by transcription factor dynamics. Stem Cell Reports. 2017;8(4):947–60.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  49. James AW. Review of signaling pathways governing MSC osteogenic and adipogenic differentiation. Scientifica. 2013;2013:684736.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  50. Soung do Y, Dong Y, Wang Y, Zuscik MJ, Schwarz EM, O’Keefe RJ, Drissi H: Runx3/AML2/Cbfa3 regulates early and late chondrocyte differentiation. J Bone Mineral Research 2007, 22(8):1260–1270.

  51. Kamata M, Okitsu Y, Fujiwara T, Kanehira M, Nakajima S, Takahashi T, Inoue A, Fukuhara N, Onishi Y, Ishizawa K, et al. GATA2 regulates differentiation of bone marrow-derived mesenchymal stem cells. Haematologica. 2014;99(11):1686–96.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  52. Zhang J, Liang Q, Lei Y, Yao M, Li L, Gao X, Feng J, Zhang Y, Gao H, Liu DX, et al. SOX4 induces epithelial-mesenchymal transition and contributes to breast cancer progression. Cancer Res. 2012;72(17):4597–608.

    Article  CAS  PubMed  Google Scholar 

  53. Xu L, Huang S, Hou Y, Liu Y, Ni M, Meng F, Wang K, Rui Y, Jiang X, Li G. Sox11-modified mesenchymal stem cells (MSCs) accelerate bone fracture healing: Sox11 regulates differentiation and migration of MSCs. FASEB J. 2015;29(4):1143–52.

    Article  CAS  PubMed  Google Scholar 

  54. Lei T, Bi Y, Gao MJ, Gao SM, Zhou LL, Zheng HL, Chen XD. HES1 inhibits adipogenesis of porcine mesenchymal stem cells via transcriptional repression of FAD24. Domest Anim Endocrinol. 2013;45(1):28–32.

    Article  CAS  PubMed  Google Scholar 

  55. Kyurkchiev D, Bochev I, Ivanova-Todorova E, Mourdjeva M, Oreshkova T, Belemezova K, Kyurkchiev S. Secretion of immunoregulatory cytokines by mesenchymal stem cells. World J Stem Cells. 2014;6(5):552–70.

    Article  PubMed  PubMed Central  Google Scholar 

  56. Klingemann H, Matzilevich D, Marchand J. Mesenchymal stem cells - sources and clinical applications. Transfusion Med Hemotherapy. 2008;35(4):272–7.

    Article  Google Scholar 

  57. Hass R, Kasper C, Bohm S, Jacobs R. Different populations and sources of human mesenchymal stem cells (MSC): a comparison of adult and neonatal tissue-derived MSC. Cell Communication Signaling. 2011;9:12.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  58. Galipeau J, Krampera M, Barrett J, Dazzi F, Deans RJ, DeBruijn J, Dominici M, Fibbe WE, Gee AP, Gimble JM, et al. International Society for Cellular Therapy perspective on immune functional assays for mesenchymal stromal cells as potency release criterion for advanced phase clinical trials. Cytotherapy. 2016;18(2):151–9.

    Article  CAS  PubMed  Google Scholar 

  59. Meisel R, Zibert A, Laryea M, Gobel U, Daubener W, Dilloo D. Human bone marrow stromal cells inhibit allogeneic T-cell responses by indoleamine 2,3-dioxygenase-mediated tryptophan degradation. Blood. 2004;103(12):4619–21.

    Article  CAS  PubMed  Google Scholar 

  60. Chinnadurai R, Copland IB, Ng S, Garcia M, Prasad M, Arafat D, Gibson G, Kugathasan S, Galipeau J. Mesenchymal stromal cells derived from Crohn’s patients deploy indoleamine 2,3-dioxygenase-mediated immune suppression, independent of autophagy. Molecular Therapy. 2015;23(7):1248–61.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  61. Kim DS, Jang IK, Lee MW, Ko YJ, Lee DH, Lee JW, Sung KW, Koo HH, Yoo KH. Enhanced immunosuppressive properties of human mesenchymal stem cells primed by interferon-gamma. EBioMedicine. 2018;28:261–73.

    Article  PubMed  PubMed Central  Google Scholar 

  62. Xu F, Rychnovsky SD, Belani JD, Hobbs HH, Cohen JC, Rawson RB. Dual roles for cholesterol in mammalian cells. Proc Natl Acad Sci U S A. 2005;102(41):14551–6.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  63. Guan Q, Li Y, Shpiruk T, Bhagwat S, Wall DA. Inducible indoleamine 2,3-dioxygenase 1 and programmed death ligand 1 expression as the potency marker for mesenchymal stromal cells. Cytotherapy. 2018;20(5):639–49.

    Article  CAS  PubMed  Google Scholar 

  64. Zhu C, Anderson AC, Schubart A, Xiong H, Imitola J, Khoury SJ, Zheng XX, Strom TB, Kuchroo VK. The Tim-3 ligand galectin-9 negatively regulates T helper type 1 immunity. Nat Immunol. 2005;6(12):1245–52.

    Article  CAS  PubMed  Google Scholar 

  65. Gil-Yarom N, Radomir L, Sever L, Kramer MP, Lewinsky H, Bornstein C, Blecher-Gonen R, Barnett-Itzhaki Z, Mirkin V, Friedlander G, et al. CD74 is a novel transcription regulator. Proc Natl Acad Sci U S A. 2017;114(3):562–7.

    Article  CAS  PubMed  Google Scholar 

  66. Fujimoto Y, Tu L, Miller AS, Bock C, Fujimoto M, Doyle C, Steeber DA, Tedder TF. CD83 expression influences CD4+ T cell development in the thymus. Cell. 2002;108(6):755–67.

    Article  CAS  PubMed  Google Scholar 

  67. Tokunaga R, Zhang W, Naseem M, Puccini A, Berger MD, Soni S, McSkane M, Baba H, Lenz HJ. CXCL9, CXCL10, CXCL11/CXCR3 axis for immune activation - a target for novel cancer therapy. Cancer Treat Rev. 2018;63:40–7.

    Article  CAS  PubMed  Google Scholar 

  68. Pittenger MF, Mackay AM, Beck SC, Jaiswal RK, Douglas R, Mosca JD, Moorman MA, Simonetti DW, Craig S, Marshak DR. Multilineage potential of adult human mesenchymal stem cells. Science. 1999;284(5411):143–7.

    Article  CAS  PubMed  Google Scholar 

  69. Crisan M, Corselli M, Chen WC, Peault B. Perivascular cells for regenerative medicine. J Cell Mol Med. 2012;16(12):2851–60.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  70. James AW, Zara JN, Corselli M, Askarinam A, Zhou AM, Hourfar A, Nguyen A, Megerdichian S, Asatrian G, Pang S, et al. An abundant perivascular source of stem cells for bone tissue engineering. Stem Cells Transl Med. 2012;1(9):673–84.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  71. Ankrum JA, Ong JF, Karp JM. Mesenchymal stem cells: immune evasive, not immune privileged. Nat Biotechnol. 2014;32(3):252–60.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  72. Sun C, Wang L, Wang H, Huang T, Yao W, Li J, Zhang X. Single-cell RNA-seq highlights heterogeneity in human primary Wharton’s jelly mesenchymal stem/stromal cells cultured in vitro. Stem Cell Res Ther. 2020;11(1):149.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  73. Klinker MW, Marklein RA, Lo Surdo JL, Wei CH, Bauer SR. Morphological features of IFN-gamma-stimulated mesenchymal stromal cells predict overall immunosuppressive capacity. Proc Natl Acad Sci U S A. 2017;114(13):E2598–607.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  74. Wegmeyer H, Broske AM, Leddin M, Kuentzer K, Nisslbeck AK, Hupfeld J, Wiechmann K, Kuhlen J, von Schwerin C, Stein C, et al. Mesenchymal stromal cell characteristics vary depending on their origin. Stem Cells Dev. 2013;22(19):2606–18.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  75. McLeod CM, Mauck RL. On the origin and impact of mesenchymal stem cell heterogeneity: new insights and emerging tools for single cell analysis. European Cells Materials. 2017;34:217–31.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  76. Weise G, Lorenz M, Posel C, Maria Riegelsberger U, Storbeck V, Kamprad M, Kranz A, Wagner DC, Boltze J. Transplantation of cryopreserved human umbilical cord blood mononuclear cells does not induce sustained recovery after experimental stroke in spontaneously hypertensive rats. J Cereb Blood Flow Metab. 2014;34(1):e1–9.

    Article  PubMed  CAS  Google Scholar 

  77. Yang YK, Ogando CR, Wang See C, Chang TY, Barabino GA. Changes in phenotype and differentiation potential of human mesenchymal stem cells aging in vitro. Stem Cell Res Ther. 2018;9(1):131.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  78. Hladik D, Hofig I, Oestreicher U, Beckers J, Matjanovski M, Bao X, Scherthan H, Atkinson MJ, Rosemann M. Long-term culture of mesenchymal stem cells impairs ATM-dependent recognition of DNA breaks and increases genetic instability. Stem Cell Res Ther. 2019;10(1):218.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  79. Jansen BJ, Gilissen C, Roelofs H, Schaap-Oziemlak A, Veltman JA, Raymakers RA, Jansen JH, Kogler G, Figdor CG, Torensma R, et al. Functional differences between mesenchymal stem cell populations are reflected by their transcriptome. Stem Cells Dev. 2010;19(4):481–90.

    Article  CAS  PubMed  Google Scholar 

  80. Wagner W, Wein F, Seckinger A, Frankhauser M, Wirkner U, Krause U, Blake J, Schwager C, Eckstein V, Ansorge W, et al. Comparative characteristics of mesenchymal stem cells from human bone marrow, adipose tissue, and umbilical cord blood. Exp Hematol. 2005;33(11):1402–16.

    Article  CAS  PubMed  Google Scholar 

  81. Cho KA, Park M, Kim YH, Woo SY, Ryu KH. RNA sequencing reveals a transcriptomic portrait of human mesenchymal stem cells from bone marrow, adipose tissue, and palatine tonsils. Sci Rep. 2017;7(1):17114.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  82. Roson-Burgo B, Sanchez-Guijo F, Del Canizo C, De Las RJ. Transcriptomic portrait of human mesenchymal stromal/stem cells isolated from bone marrow and placenta. BMC Genomics. 2014;15:910.

    Article  PubMed  PubMed Central  Google Scholar 

  83. Gnani D, Crippa S, Della Volpe L, Rossella V, Conti A, Lettera E, Rivis S, Ometti M, Fraschini G, Bernardo ME, et al. An early-senescence state in aged mesenchymal stromal cells contributes to hematopoietic stem and progenitor cell clonogenic impairment through the activation of a pro-inflammatory program. Aging Cell. 2019;18(3):e12933.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  84. Ren G, Zhang L, Zhao X, Xu G, Zhang Y, Roberts AI, Zhao RC, Shi Y. Mesenchymal stem cell-mediated immunosuppression occurs via concerted action of chemokines and nitric oxide. Cell Stem Cell. 2008;2(2):141–50.

    Article  CAS  PubMed  Google Scholar 

  85. Ren G, Su J, Zhang L, Zhao X, Ling W, L’Huillie A, Zhang J, Lu Y, Roberts AI, Ji W, et al. Species variation in the mechanisms of mesenchymal stem cell-mediated immunosuppression. Stem Cells. 2009;27(8):1954–62.

    Article  CAS  PubMed  Google Scholar 

  86. Yuan X, Qin X, Wang D, Zhang Z, Tang X, Gao X, Chen W, Sun L. Mesenchymal stem cell therapy induces FLT3L and CD1c(+) dendritic cells in systemic lupus erythematosus patients. Nat Commun. 2019;10(1):2498.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  87. Diehl R, Ferrara F, Muller C, Dreyer AY, McLeod DD, Fricke S, Boltze J. Immunosuppression for in vivo research: state-of-the-art protocols and experimental approaches. Cell Mol Immunol. 2017;14(2):146–79.

    Article  CAS  PubMed  Google Scholar 

Download references

Acknowledgements

We thank Hailun Wang for his reviews and suggestions. We also thank the members of the China National GeneBank for technical supports.

Funding

This work was funded by Science, Technology and Innovation Commission of Shenzhen Municipality (No. KQJSCX20170322143848413) and Shenzhen Municipal Government of China (No. 20170731162715261). It was also supported by BGI-Shenzhen.

Author information

Author notes

  1. Changbin Sun and Kehua Zhang contributed equally to this work.

Authors and Affiliations

  1. BGI-Shenzhen, Jinsha Road, Dapeng New District, Shenzhen, 518083, China

    Changbin Sun, Jianhui Yue & Xi Zhang

  2. BGI Education Center, University of Chinese Academy of Sciences, Shenzhen, 518083, China

    Changbin Sun

  3. China National GeneBank, BGI-Shenzhen, Shenzhen, 518120, China

    Changbin Sun, Jianhui Yue & Xi Zhang

  4. Cell Collection and Research Center, National Institutes for Food and Drug Control, Beijing, 100050, China

    Kehua Zhang & Shufang Meng

  5. Section of Cell Biology and Physiology, Department of Biology, University of Copenhagen, Copenhagen, Denmark

    Jianhui Yue

Authors
  1. Changbin Sun
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Kehua Zhang
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Jianhui Yue
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Shufang Meng
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Xi Zhang
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

Conceptualization, C.S., K.Z., and X.Z. Methodology and investigation, C.S., K.Z., J.Y., S.M., and X.Z. Bioinformatic analysis, C.S. Writing, C.S., K.Z., and X.Z. All authors read and approved the manuscript.

Corresponding author

Correspondence to Xi Zhang.

Ethics declarations

Ethics approval and consent to participate

This study was approved by the Ethics Committee of BGI-IRB (NO. BGI-IRB 19025-1)

Consent for publication

Written informed consent was obtained from the donors for recording/using their individual details (age, gender, weight, and height) for research purposes.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1: Fig. S1.

Workflow for data search and selection. Fig. S2. Overview of data collected for transcriptome variation analysis. (A) Boxplot showing correlation coefficient of transcriptome expression for each sample to others in the dataset before data selection step 3. 0.9 was selected as a cutoff for filtering out samples with lower median Pearson’s r.(B) Barplot showing number of samples in each study. (C) Histograms showing number of reads (up) and mapping results (down) across samples. Some MSCs derived from different anatomical parts of placenta, AM CM, CP, and CV. AD: Adipose tissue; AM: Amniotic membrane; BM: Bone marrow; CM: Chorionic membrane; CP: Chorionic plate; CV: Chorionic villi; ED: Endometrial; DP: Dental pulp; PL: Placenta; UC: Umbilical cord. Fig. S3. Transcriptome variation across MSC samples. (A) Scatter plot showing DM value and mean expression (CPM) for each expressed gene. (B) Scatter plot showing DM value and gene length for each expressed gene. (C) GESA positive results showing enrichment in GO cellular component (left) and molecular function (right) gene sets database based on ranked genes list in descending order by the DM value. (D) GESA negative results showing enrichment in GO gene sets database based on ranked genes list in descending order by the DM value. Only the top 20 terms with highest NES were presented (p < 0.001). Fig. S4. Differential gene expression and function enrichment analysis of MSCs among G0, G2, G3 and G4. (A) Results of GO biological process enrichment analysis for genes upregulated in G0. (B) Representative genes upregulated in G0. (C) Results of GO biological process enrichment analysis for genes upregulated in G2. (D) Representative genes upregulated in G2. (E) Results of GO biological process enrichment analysis for genes upregulated in G3. (F) Representative genes upregulated in G3. (G) Results of GO biological process enrichment analysis for genes upregulated in G4. (H) Representative genes upregulated in G4. Fig. S5. Expression variability of MSCs in response to IFNγ. (A) Barplot showing genes responded differently to the IFNγ among ssG0, ssG2, and ssG4. (B) Venn diagraming showing overlap of genes responded differently to the IFNγ among ssG0, ssG2, and ssG4. (C-E) Barplot showing enrichment of genes responded differently to the IFNγ for ssG0 (C), ssG2(D), and ssG4 (E). Fig. S6. PCA visualizing distance of samples based on different gene panels. (A) Scatter plot visualizing PCA results using all genes in Genes set1 expressed in the untreated samples. (B) Scatter plot visualizing PCA results using all genes in Genes set2 expressed in the untreated samples. (C) Scatter plot visualizing PCA results using all genes in Genes set3 expressed in the untreated samples. (D) Scatter plot visualizing PCA results using HVGs expressed in the untreated samples. (E) Scatter plot visualizing PCA results using random sampling 100 genes expressed in the untreated samples. (F) Scatter plot visualizing PCA results using top 100 genes with the highest DM values in the Genes set2 expressed in the untreated samples and untreated samples.

Additional file 2: Table S1.

Table S1 Information of MSCs RNA-seq samples downloaded from public GEO. Table S2. DM values for each gene expressed in MSC samples. Table S3. Information of samples clustered into different groups. Table S4 Results of different expression analysis among G0, G2, G3, and G4 groups. Table S5 Results of different expression analysis for INFγ treated vs untreat MSC paried samples. Table S6. List of genes differently respond to IFNγ treatment among ssG0, ssG2, and ssG4. Table S7. Overlap of genes between DRGs and CRGs.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Sun, C., Zhang, K., Yue, J. et al. Deconstructing transcriptional variations and their effects on immunomodulatory function among human mesenchymal stromal cells. Stem Cell Res Ther 12, 53 (2021). https://doi.org/10.1186/s13287-020-02121-8

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s13287-020-02121-8

Keywords

Stem Cell Research & Therapy

ISSN: 1757-6512