Single-cell RNA-seq highlights heterogeneity in human primary Wharton's jelly mesenchymal stem/stromal cells cultured in vitro

Abstract

Background: Mesenchymal stem/stromal cells (MSCs) are multipotent cells with a promising application potential in regenerative medicine and immunomodulation. However, MSCs cultured in vitro exhibit functional heterogeneity. The underlying molecular mechanisms that define MSC heterogeneity remain unclear.

Methods: We investigated the gene expression profile via single-cell RNA sequencing (scRNA-seq) of human primary Wharton's jelly-derived MSCs (WJMSCs) cultured in vitro from three donors. We also isolated CD142⁺ and CD142^- WJMSCs based on scRNA-seq data and compared their proliferation capacity and "wound healing" potential in vitro. Meanwhile, we analyzed publicly available adipose-derived MSC (ADMSCs) scRNA-seq data and performed transcriptome comparison between WJMSCs and ADMSCs at the single-cell level.

Results: GO enrichment analysis of highly variable genes (HVGs) obtained from WJMSCs revealed that these genes are significantly enriched in extracellular region with binding function, involved in developmental process, signal transduction, cell proliferation, etc. Pathway analysis showed that these HVGs are associated with functional characteristics of classic MSCs, such as inflammation mediated by chemokine and cytokine signaling, integrin signaling, and angiogenesis. After regressing out the batch and cell cycle effects, these HVGs were used for dimension reduction and clustering analysis to identify candidate subpopulations. Differentially expressed gene analysis revealed the existence of several distinct subpopulations of MSCs that exhibit diverse functional characteristics related to proliferation, development, and inflammation response. In line with our data, sorted CD142⁺ and CD142^- WJMSCs showed distinct proliferation capacity as well as "wound healing" potential. Although WJMSCs and ADMSCs were derived from different tissues and were displaying different differentiation potencies, their HVGs were largely overlapped and had similar functional enrichment.

Conclusion: HVGs identified in MSCs are associated with classic MSC function. Regarding therapeutic potential, these genes are associated with functional characteristics, on which the MSC clinical application were theoretically based, such as development and inflammation response. Altogether, these HVGs hold the potential to be used as candidate markers for further potency association studies.

Keywords: Functional enrichment; Highly variable genes; MSCs; Potency; scRNA-seq.

Fig. 1

Overview of WJMSCs single-cell RNA-seq data. a Expression of marker genes in the three samples. Number on the top showing percentage of cells with at least one UMI. b Boxplot showing top 50 cluster of differentiation (CD) genes ranked by average normalized expression. c Distribution of UMI cross cells after pre-processing to filter out low-quality cells. d Distribution of expressed genes after pre-processing to filter out low-abundance genes with mean-based method (genes with means more than 0.1 were retained)

Fig. 2

Heterogeneity and highly variable genes in WJMSCs. a Phases of cell cycle assigned for each of the three samples. b, c Cell cycle effects (b) and batch effects (c) represent the dominant source of heterogeneity in primary cultured WJMSC population. d Results of pathway enrichment analysis for highly variable genes identified in WJMSCs. e Results of protein class enrichment analysis for highly variable genes identified in WJMSCs. IMMC, inflammation mediated by chemokine and cytokine; GRHR, gonadotropin-releasing hormone receptor

Fig. 3

Candidate subpopulations with different functional characteristics. a UMAP visualizing the results of cell clustering. b Bean plots showing expression of several representative DEGs among the six subpopulations. c Pathways significantly enriched for the genes differentially expressed in one subpopulation compared to others. d GO-slim biological process enriched for the genes differentially expressed in one subpopulation compared to others. For c and d, only the top 10 terms with lowest FDR (FDR ≤ 0.05) were present. e Boxplots showing expression of classic MSC marker genes in subpopulations. f Example of candidate markers showing different expression pattern among the five subpopulations (C0–C4). C0, red; C1, olive; C2, yellow; C3, green; C4, blue; C5, purple. IMCC, inflammation mediated by chemokine and cytokine; DNPD, de novo pyrimidine deoxyribonucleotide

Fig. 4

Candidate subpopulations showing different predicted potencies on differentiation and immunosuppression. a Boxplots showing expression of genes correlated with PBMC suppression across the five candidate subpopulations (C0–C4). b–f Bean plots showing distribution of log (norm_exprs) values of osteogenic score (b), chondrogenic score (c), adipogenic score (d), myogenic score (e), and neurogenic score (f) across the five candidate subpopulations (C0–C4). Wilcoxon rank sum test were performed for significant test, ** P < 0.001

Fig. 5

Transcriptome comparison at the single-cell level between WJMSCs and ADMSCs. a Boxplot showing the top 50 CD genes ranked by average normalized expression in ADMSCs. b Venn diagram showing top 50 ADMSC CD genes overlap with the top 50 WJMSC CD genes, unshared genes were highlighted. c Pathway enrichment of top 50 CD genes expressed in ADMSCs and WJMSCs. d Example of CD genes showing different expression percentage between ADMSCs and WJMSCs. e Violin plots showing expression of markers for reported MSC subpopulations in ADMSCs and WJMSCs. f, g Density distribution showing osteogenic score (f) and adipogenic score (g) between ADMSCs and WJMSCs. Percentage indicating proportion of cells assigned to the right side of the line