To identify the transcriptional changes resulting from growth in spheroid vs. monolayer conditions we performed single cell RNA sequencing (scRNA-seq) on cells grown under both conditions. After quality control steps to remove presumed low viability cells (>25% mitochondrial transcripts) and cell doublets, a total of 11,515 cells (8,022 from spheroid culture, and 3,493 from monolayer culture) were retained for analysis. Data from the two libraries were combined using the Seurat integration pipeline to group similar cells from both libraries and the data were visualized using Uniform Manifold Approximation and Projection (UMAP). A nearest neighbor graph was constructed using the integrated assay and clustering was performed using the Louvain algorithm at a resolution of 0.3, identifying eight clusters (
Figure 4A). Coloring cells by source library in the UMAP projection shows that the integration procedure intermixed cells from the two libraries in UMAP space, compared to a UMAP projection generated on the non-integrated data, in which the cells from the two libraries are separate in UMAP space (
Figure S3A). The cells of both libraries can also be split into populations with a high (>2,500) or low (≤2,500) number of detected genes, with an increased number of lower gene count cells from the spheroid culture (
Figure S3B).The low and high gene cells group separately in UMAP space and Louvain clustering, with clusters 0, 1, 2, and 6 containing the low gene cells, and clusters 3, 4, 5, and 7 containing high gene cells. Identification of marker genes for each cluster found between 233 and 1,671 markers per cluster (with an adjusted p-value <0.05), with a total of 3,365 genes identified as markers of at least one cluster. In these clusters, spheroid cells comprise the majority of clusters 0, 1, 2, 4 and 7 (
Figure S3C). The heatmap suggested differential gene expression across all clusters (
Figure S3D), which supports their separate distribution in UMAP (
Figure 4B). As seen in
Figure S3D, clusters 1 and 6 grouping together in the top dendrogram, as do clusters 4 and 7; these cluster pairs were both grouped as a single cluster when the clustering was run at a lower resolution. Analysis of differentially expressed genes (DEGs) revealed that each cluster was characterized by a specific transcriptional profile (Supplementary Excel spreadsheet 1).
Next we examined expression of the MSC surface markers observed in IF staining (
Figure 2C) in each cluster. Consistent with the IF results in
Figure 2C, the expression levels of mouse MSC genes,
Sca-1 (
Ly6a),
CD29 (
Itgb1),
CD44,
CD90.1 (
Thy1) were high in clusters 2, 3, 4, 5 and 7 (
Figure 4C and
Figure S4). On the basis of their gene expression profile, we focused on clusters 4 and 7 for further analysis. Our gene expression analyses revealed that cluster 4, which includes 9% of 8,022 original spheroid cells (
Figure 4D), exhibited significantly increased expression of MSC markers compared to the rest of the cells (
Supplementary Excel spreadsheet 2, and
Figure 4E). We observed that the expression of many extracellular matrix (ECM) genes, including
Col1a1, Col1a2, Col3a1, Fbn1, Lama2, and Lama4, were significantly greater in cells of one of more of clusters 2, 3, 4, 5 and 7 (
Table S1 and
Figure S5A, B), and that expression of some of these genes was significantly higher in spheroid than in attached cells. Further, their expression was significantly increased in spheroid cells of cluster 4 compared with the rest of the cells (
Figure 4F). Since approximately 10% of spheroids cells showed self-renewal property (
Figure 2F), we sought to explore the genes important for stem cell self-renewal. Analysis of differentially expressed genes revealed significantly increased expression of one or more of the genes
Notch2, Sox4, Sox9, Klf2, and
Foxp1 [
22,
23,
24,
25,
26] in clusters 3, 4, 5 and 7 (
Table S2 and
Figure S6A), but only in cluster 4 was expression all of these genes significantly greater when compared to the remainder of the dataset, with
Sox9 and
Foxp1 also significantly higher in spheroid vs monolayer cells within cluster 4 (
Table S2 and
Figure 4G and
Figure S6B). The heatmap shown in
Figure S6C displays expression profile of selected self-renewal and collagen genes across clusters. The DEGs analysis of cluster 7, which had the lowest percentage of spheroid cells among the clusters (
Figure 4D), showed significant upregulation of a family of transcription factors,
Hox genes [
27], including
Hoxb7, Hoxb9, Hoxb13, and
Hoxa11os (
Figure 4H and
Figure S7A) and
Mmp10 and
Mmp13 genes which are implicated in ECM remodeling [
28] (
Figure S7B, C). The heatmaps demonstrating average normalized expression per cluster of genes expressed at a significantly higher level in cluster 7 and all the
Hox and
Mmp genes in the dataset is presented in
Figure S8A, B. Altogether, characterization of transcriptional profiles of the TEFs-derived spheroid cells using scRNA sequencing revealed the presence of specific clusters associated with mesenchymal stem-like cells and self-renewing capacity.