Stem cells, which possess remarkable abilities of self-renewal and differentiation, are classified into embryonic stem cells, adult stem cells, and induced pluripotent stem cells (1). Among adult stem cells, mesenchymal stem cells (MSCs) have gained significant attention due to their ethical advantages over embryonic stem cells, and their ability to be isolated from various tissues, including bone marrow, umbilical cord blood, periodontal ligament, dental pulp, synovial membrane, and adipose tissue (2, 3). In particular, adipose tissue stands out as a rich source of MSCs that can be easily obtained through minimally invasive procedures, such as liposuction (4, 5). Adipose tissue-derived mesenchymal stem cells (AT-MSCs) demonstrate stable proliferation in culture, and exhibit the potential to differentiate into both mesodermal lineages, such as osteocytes, chondrocytes, adipocytes, and muscle cells, and ectodermal and endodermal lineages (3). Consequently, AT-MSCs hold great promise as a valuable cell source for tissue regeneration and regenerative therapies (3).

Since the early 2000s, scientists have been extensively exploring the potential of AT-MSCs in treating various disorders, including articular cartilage lesions (6). In the field of cartilage tissue engineering, both preclinical studies and clinical trials have investigated the use of AT-MSCs (6, 7). However in this area, addressing the inherent heterogeneity of AT-MSCs remains a challenging task. As with other types of MSCs, the heterogeneity of AT-MSCs arises from individual variations among donors, differences in tissue sources and cell subset compositions, and inconsistencies in culture conditions (8, 9). Factors such as donor age and health status both directly impact the availability of cells, and influence their proliferation, colony-forming ability, differentiation potential, and therapeutic effectiveness (10). A comprehensive understanding of cellular heterogeneity is crucial for achieving efficient and reproducible clinical applications of AT-MSCs (11). However, the conventional bulk cell analysis and limited reliance on a few cell surface markers fall short of fully capturing the intricate nature of cellular heterogeneity.

Fortunately, recent advances in single-cell RNA sequencing (scRNA-seq) have resulted in this technique emerging as a powerful tool for dissecting cellular heterogeneity within cell populations at the single-cell level in an unbiased manner (12). The scRNA-seq technique is categorized based on several factors that include single-cell isolation, cDNA amplification, library preparation, and sequencing platform (13). Various approaches to single-cell isolation are currently available, such as limiting dilution, micromanipulation, laser capture microdissection, magnet-activated cell sorting, fluorescence-activated cell sorting, and microfluidics (14). Among these, microfluidics, a novel technology that allows precise control of fluid at the micro-scale, can be further classified into trap-based microfluidics, valve-based microfluidics, and droplet-based microfluidics (14).

In this study, we employed drop-seq, a droplet-based scRNA-seq technique (15), to investigate the response of AT-MSCs to chondrogenic induction at the single-cell level. To specifically focus on the response to chondrogenic induction, and exclude the various factors that were mentioned earlier and that can contribute to the heterogeneity of AT-MSCs, we utilized commercially available AT-MSCs. Our analysis revealed the coexistence of four distinct clusters of cells displaying heterogeneity in translational machinery biogenesis, cell proliferation, and metabolic preferences. These complexities could not be captured by the conventional bulk analyses that rely on averaging response and/or expression levels.

Traditional bulk analysis shows the mean or main responses of an AT-MSC population to chondrogenic induction

Commercially available AT-MSCs were cultured in chondrogenic differentiation media for 1 and 2 weeks, designated as 1W and 2W, respectively. Subsequently, the cells were stained with Alcian Blue (Fig. 1A). In contrast to the undifferentiated cells before induction (designated as the control), glycosaminoglycans, which are important components of cartilage, were distinctly observed in the 1W and 2W cells. Histologically, no apparent differences were detected between the 1W and 2W samples.

We performed bulk RNA-seq, and compared the differentially expressed genes (DEGs) of 1W and 2W cells with those of control cells. As shown in Fig. 1B, 2,482, 2,482 and 2,734,734 genes were up-regulated with at least 2-fold change in 1W and 2W cells, respectively. The number of down-regulated genes in 1W and 2W cells was (1,893 and 2,257), respectively. Gene Ontology of Biological Processes (GO-BP) enrichment analysis was performed to gain insight into the biological significance of these DEGs. The DEGs in 1W and 2W cells were enriched in (109 and 114) GO-BP terms, respectively. Among these, 90 GO-BP terms were common to both 1W and 2W cells. Specific DEGs in 1W and 2W cells were enriched in (19 and 24) GO-BP terms, respectively (File S1 of the Supplementary Information [SI]).

Fig. 1C illustrates the representative GO-BP terms commonly enriched in 1W and 2W cells, which include the up-regulation of GO-BP terms related to cell adhesion, extracellular matrix organization, and collagen fibril organization, as well as the down-regulation of GO-BP terms related to cell division, DNA replication, and DNA repair. Representative GO-BP terms specifically up-regulated and down-regulated in 1W cells include calcium-dependent cell-cell adhesion via plasma membrane cell adhesion molecules and mRNA splicing via spliceosome, respectively (File S1 of the SI). Representative GO-BP terms specifically up-regulated and down-regulated in 2W cells were identified as the transmembrane receptor protein tyrosine kinase signaling pathway and actin cytoskeleton organization, respectively (File S1 of the SI).

We investigated the statistical significance (−log10 [adjusted P-value]) of expression changes in 14 genes known as chondrogenic markers (16, 17), based on bulk RNA-seq data, as illustrated in Fig. 1D. Most of the analyzed genes exhibited a consistent expression pattern in 1W and 2W. Most, though not all, of the stage 1 markers (COL1A1, SOX4), stage II markers (COMP, COL11A1, SOX9), stage III markers (MATN3, WNT11), and stage IV markers (COL10A1, FMOD) displayed significant up-regulation at either 1W, 2W, or both time points. Conversely, COL6A1 and PTHLH, classified as stage I and IV markers respectively, demonstrated down-regulation. The expression level of HOXA7 (stage III marker) exhibited no significant change compared to the control cells, while transcripts of CHAD and COL1A1 (classified as stage III and IV markers, respectively) were expressed at very low levels.

Quality control of the scRNA-seq data

To investigate the response of AT-MSCs to chondrogenic induction at the single-cell level, we performed scRNA-seq. To assess the data quality, we examined various parameters, including the number of genes per cell, the number of transcripts, and the percentage of expressed mitochondrial genes (File S2 of the SI). The number of detected genes can decrease if cells are missed during capture, or if broken beads are present in the droplets. Conversely, if two or more cells are accidentally captured in a single droplet, the number of detected genes increases accordingly. Therefore, we excluded cells with too few (less than 200) or too many (more than 4,000) detected genes. Additionally, cells with a mitochondrial gene expression rate of 10% or higher were excluded from the analysis, as this indicates potential issues, such as cell death during the isolation process, or improper dissolution of cells within the encapsulated droplets (18, 19). After applying these quality control measures, we retained 37,219 transcripts from 378 control cells, 433 1W cells, and 510 2W cells.

scRNA-seq reveals the dynamics of cell cluster proportions

The cell cycle is widely recognized as biologic noise that contributes to cellular heterogeneity (11, 20, 21). Building upon the study by Tirosh et al. (21), we performed principal component analysis (PCA) that considered or excluded genes significantly influenced by the cell cycle (referred to as cell-cycle genes), and compared the results. As depicted in Fig. 2A (upper), the inclusion of these genes in the analysis revealed clustering patterns that were indicative of associations with cell cycle phases. Conversely, when these genes were excluded from the analysis, cells inferred to be in the same cell cycle phase were observed to be randomly dispersed throughout the population (Fig. 2A, lower). To focus specifically on chondrogenic differentiation and minimize potential biological noise introduced by the cell cycle, we excluded the cell-cycle genes, and proceeded with further analysis.

Based on the significance levels (P-values) and standard deviations of each component (File S3 of the SI), we selected seven PCs out of a total of twenty. Utilizing these components, we generated a uniform manifold approximation and projection (UMAP) (Fig. 2B), which revealed the presence of four distinct cell clusters (clusters (0-3). The proportions occupied by these clusters varied across samples.

Fig. 2C presents the count and percentage of cells allocated to each cluster. Cluster 1 exhibited prominence in both control cells and 1W cells, constituting approximately 43.9 and 34.9%, respectively. In contrast, clusters 3, 2, and 0 were exclusively dominant in control, 1W, and 2W cells, accounting for approximately 37.0, 47.3, and 81.4%, respectively.

scRNA-seq reveals heterogenous chondrogenic marker expression

Fig. 3 presents violin plots illustrating the expression levels of chondrogenic markers based on the scRNA-seq results. Stage III and IV markers (MATN3, WNT11, HOX7, COL10A1, FMOD, and PTHLH) consistently exhibited low expression in the majority of cells. In contrast, stage I and II makers (COL1A1, SOX4, COL6A1, COMP, and COL11A1), except for SOX9, displayed a higher proportion of cells with relatively elevated expression levels following chondrogenic stimulation (1W or 2W cells), compared to control cells (left panels). Notably, the expression levels of these markers showed heterogeneity, even within the same cluster (right panels).

scRNA-seq unveils transcriptional heterogeneity in key biological processes among clusters

We compared each cluster with the others to investigate DEGs. Additionally, we conducted GO-BP enrichment analysis to gain insights into the biological significance of the DEGs within each cluster. Fig. 4A depicts the number of DEGs and enriched GO-BP terms, with the complete information provided in File S4 of the SI. Fig. 4B-E illustrate noteworthy GO-BP terms that varied across clusters: (1) Cluster 1 displayed up-regulation in ribosome production and translation (Fig. 4B); (2) Furthermore, cluster 1 exhibited up-regulation of GO-BP terms associated with mitochondrial oxidative metabolism, such as aerobic respiration, oxidative phosphorylation, and mitochondrial electron transport (Fig. 4C, File S4 of the SI); (3) Cluster 3 showed up-regulation of GO-BP terms associated with cell proliferation, such as DNA replication and cell division (Fig. 4D); (4) Clusters 0 and 2 shared similarities, both demonstrating the up-regulation of collagen fibril organization and the down-regulation of GO-BP terms related to points (2) and (3) (Fig. 4C-E). However, ribosome production was significantly down-regulated only in cluster 0 (Fig. 4B).

Heterogeneous responses to chondrogenic induction stimulation have been well-documented among MCSs derived from various tissues (22). Our study aimed to comprehensively investigate the chondrogenic response of AT-MSCs. To accomplish this, we employed both traditional bulk analyses and scRNA-seq. In bulk analyses, which provide average responses across all cells, we observed a seemingly uniform cellular response in each sample. Expression levels of most chondrogenic differentiation markers from stages 1 to 4 were up-regulated in 1W and 2W cells, compared to control cells (Fig. 1).

However, scRNA-seq, which examines individual single-cell transcriptomes, revealed primarily elevated expression of stage 1 and 2 markers. As shown in Fig. 3, the expression of stage 3 and 4 markers was elevated in a small subset of 1W and 2W cells, indicating that the average values derived from the bulk RNA-seq can lead to misinterpretation by representing the responses of these few cells as representative of the entire cell population. Therefore, a reasonable interpretation of our results suggests that AT-MSCs exposed to chondrogenic induction stimuli (1W and 2W cells) reached stage 1 and 2 of differentiation. This interpretation is further supported by the up-regulation of GO-BP terms related to collagen fibril organization and supramolecular fiber organization in clusters 2 and 0, which dominate in 1W and 2W cells (Fig. 4E). Our findings underscore that relying solely on the conventional bulk RNA-seq for chondrogenic marker analysis may inaccurately assess differentiation steps.

Furthermore, scRNA-seq unveiled marker expression heterogeneity at the single-cell level, even within the same cell group (Fig. 3). For example, in 2W cells, COL1A1 expression displayed a varied range, with some cells exhibiting low levels. This variability could pose a challenge for bulk RNA-seq in detecting significant differences in COL1A1 expression between control and 2W cells (Fig. 1). Additionally, the narrow and elongated violin plots for 1W and 2W cells indicate significant cell-to-cell variation in COL6A1 expression, potentially leading to a misinterpretation of mean expression in the bulk analysis (Figs. 1 and 3). Unlike other stage 2 markers, SOX9 consistently maintained low expression in most cells. This suggests a staggered up-regulation of stage markers, and is consistent with the low expression of COL2A1, which is influenced by SOX9 (23). Similar to stage 3 and 4 markers, SOX9 expression changed noticeably in only a small number of cells within each cell group, which seems to be incorrectly reflected as increased or decreased expression in overall cells in the bulk RNA-seq.

We conducted a comprehensive investigation into the characteristics of each cell cluster, classified based on the similarity of single-cell transcripts. Clusters 1 and 3, predominant among undifferentiated AT-MSCs (control cells), exhibited distinct traits. While cluster 3 exhibited the up-regulation of GO-BP terms related to cell proliferation (Fig. 4D), cluster 1 displayed the significant up-regulation of GO-BP terms associated with ribosome production and oxidative mitochondrial metabolism (Fig. 4B, C). This indicates a distinct metabolic preference for cluster 1, compared to the other clusters. Importantly, even after 1 week of chondrogenic induction, cluster 1 maintained a relatively high proportion, whereas cluster 3 did not (Fig. 2C), implying differing responses or rates of the clusters to chondrogenic induction stimulation. The decrease in the proportion of cluster 3 in 1W and 2W cells (Fig. 2C) may be the main reason for the bulk analysis results showing the down-regulation of cell division and DNA replication in both 1W and 2W cells (Fig. 1C). Furthermore, cluster 0, predominant in 2W cells, exhibited a significant reduction in GO-BP terms associated with ribosome production, in contrast to cluster 2, which dominated in 1W cells (Fig. 4B). Given that cell proliferation and differentiation exhibit an inverse relationship (24), this result suggests a potential correlation between the ribosome regulation, proliferation, and differentiation status of AT-MSCs.

Stem cells possess an impressive ability to sustain self-renewal while retaining the potential for diverse cell lineage differentiation (25). An elevated presence of inactive ribosomes and reduced translational efficiency, compared to differentiated progenies, has been considered to be a shared characteristic among stem cells, playing a pivotal role in maintaining stem cell homeostasis (25). Maintaining a substantial reservoir of accessible ribosomes is likely critical for the rapid remodeling and facilitation of gene expression toward differentiation programs, as observed in embryonic stem cells (25). Since this process incurs significant energy costs, it is particularly interesting, in terms of energy metabolism, that the concurrent up-regulation of GO-BP terms related to oxidative mitochondrial metabolism and ribosome production was prominent only in cluster 1.

Mitochondria are responsible for energy generators in cells, primarily through glycolysis, the TCA cycle, and oxidative phosphorylation pathways (26). The generation of cellular ATP via oxidative phosphorylation involves a series of electron transfer reactions within the mitochondrial electron transport chain (27). Remarkably, in contrast to the present study utilizing AT-MSCs, several studies utilizing bone marrow MSCs reported that during osteogenic and adipogenic differentiation, the energy acquisition route shifts from glycolysis to mitochondrial oxidative metabolism (26, 28, 29). However, it is essential to note that these previous investigations were based on bulk analyses, i.e., an ‘average values for whole cells’ approach, and did not reveal changes in energy metabolism preference at the single-cell level. In contrast to these bulk analyses, our findings unequivocally demonstrate the existence of a subset of AT-MSCs with a greater reliance on oxidative phosphorylation than other cells, which is presumed to be linked to meeting the energy requirements necessary to maintain accessible ribosome reserves. To determine whether these traits, i.e., a subset showing increased ribosome levels, and dependence on oxidative phosphorylation, are exclusive to AT-MSCs, or are prevalent among MSCs in general, single-cell level studies focusing on bioenergy conversion throughout the differentiation process should be conducted.

We acknowledge the limitation of our study in not encompassing the entire chondrogenic differentiation stage. Instead, we focused on analyzing 37,219 transcripts derived from 1,321 cells spanning the spectrum from an undifferentiated state to a 2-week chondrogenic differentiation period. Despite this constraint, our findings vividly illustrate the heterogeneity of AT-MSCs before and after chondrogenic induction, including metabolic preferences and ribosome production. Future studies are needed to explore whether these fundamental biological processes contributing to cellular heterogeneity are directly sensed by factors governing the fate of AT-MSCs, or whether they indirectly influence the translational capacity, ultimately affecting the expression of cell fate determinants. A comprehensive understanding of AT-MSC heterogeneity, attained through single-cell approaches, holds paramount importance for the efficient utilization of this exceptional cell source in practical, reproducible, and effective clinical applications.

See supplementary file for detailed procedures.

Cell culture and chondrogenic induction

Human AT-MSCs (Promocell, C-12977) were cultured using Mesenchymal Stem Cell Growth Medium 2 at 37°C with 5% CO₂. Cells from passages 3 to 5 were cultured in chondrogenic differentiation medium supplemented with 10% fetal bovine serum for 7 and 14 days.

Bulk RNA-seq

RNA-Seq libraries were prepared using TruSeq Stranded mRNA LT Sample Prep Kit. Fragments (∼300 bp) were gel-isolated, PCR amplified, and sequenced as paired-end reads (2 × 101 bp reads) on Illumina HiSeq 2500. HISAT2 aligned reads to human genome hg38. Gene annotation data were obtained from Ensembl (release 98) biomart. DEGs were identified using DESeq2 R package (≥ 2-fold change, adjusted P-value < 0.05, Benjamini-Hochberg correction). GO term enrichment analysis was performed using DAVID.

scRNA-seq

The single-cell isolation process entailed seeding cells at a concentration of 5.0 × 10⁵ cells/flask in T-25 flasks. After chondrogenic differentiation, cells were isolated, enzymatically treated, and passed through 40 μm cell strainers. Drop-seq was conducted as described previously with modifications. Reads were processed using Drop-seq_Tools pipeline v2.3.0 and mapped to hg38 with STAR v2.5.2b.

scRNA-seq data analysis

We utilized R version 4.1.2, Seurat v4.0, and MergeSeurat. Cells with gene counts below 200 or above 4,000, or with a mitochondrial gene percentage exceeding 10%, were excluded. Data were scaled (default = 10,000) and log-transformed. PCA and clustering were performed based on JackStraw plot and Elbow plot. DEGs between cell clusters were identified using FindMarkers with significant criteria (adjusted P-value < 0.05, Bonferroni correction, ≥ 1.3-fold change). Functional enrichment analysis of DEGs utilized GO terms and DAVID.

This research was supported by the National Research Foundation of Korea (NRF) funded by the Ministry of Science & ICT (NRF-2022R1F1A1071248) and the Chung-Ang University Research Scholarship Grants in 2021.

The authors have no conflicting interests.