Cellular Heterogeneity and Lineage Restriction during Mouse Digit Tip Regeneration at Single-Cell Resolution

Abstract

Innate regeneration following digit tip amputation is one of the few examples of epimorphic regeneration in mammals. Digit tip regeneration is mediated by the blastema, the same structure invoked during limb regeneration in some lower vertebrates. By genetic lineage analyses, the digit tip blastema has been defined as a population of heterogeneous, lineage-restricted progenitor cells. These previous studies, however, do not comprehensively evaluate blastema heterogeneity or address lineage restriction of closely related cell types. In this report, we present single-cell RNA sequencing of over 38,000 cells from mouse digit tip blastemas and unamputated control digit tips and generate an atlas of the cell types participating in digit tip regeneration. We computationally define differentiation trajectories of vascular, monocytic, and fibroblastic lineages over regeneration, and while our data confirm broad lineage restriction of progenitors, our analysis reveals 67 genes enriched in blastema fibroblasts including a novel regeneration-specific gene, Mest.

Keywords: blastema; digit tip regeneration; fibroblast heterogeneity; mest; single-cell RNA-seq.

Figure 1. Cellular heterogeneity of the 11dpa blastema

(A) Schematic overview of innate mouse digit tip regeneration following amputation mid-way through the terminal phalanx. Schematic of the experimental design whereby blastemas were dissected, dissociated, and single cells captured. Single cell RNA libraries were prepared and sequenced for computational analysis. (B) Unbiased single cell clustering of 7,610 high quality cells visualized by UMAP. Each dot represents a single cell and cells assigned to the same cluster are similarly colored. Cell type identities are assigned as follows: fibroblasts (clusters 0–2, 4–6, and 8), macrophages (clusters 3, 14, and 16), bone (cluster 7), vascular smooth muscle cells (VSM) (cluster 9), endothelial cells (cluster 10), monocytes (cluster 11), epithelial cells (cluster 12), Schwann cells (SC) (cluster 13), T cells (cluster 15). (C) Gene expression UMAP overlay with examples of highly expressed, cell type specific markers used to assign cluster cell identities: Bglap (bone), Krt14 (epithelial cells), Plp1 (SCs), Lyz2 (macrophages and monocytes), Pecam1 (endothelial cells), Rgs5 (vascular smooth muscle cells), Cd3g (T cells), Prrx1 (fibroblasts). Gray depicts low expression and purple depicts high expression as specified on the scale for each gene.

Figure 2. 11dpa fibroblast heterogeneity

Representative genes differentially expressed among fibroblast clusters. (A-H) Gene expression is shown by violin plot (left) where black dots represent individual cells and the colored curve shows the distribution of cells at a given expression level, and feature plot (right) where purple is high expression and grey is low expression. Cluster numbers and UMAP plots relate to data in Figure 1B. (A) Ccl2 is expressed by cells in clusters 0, 1, 5, and 8; (B) Mmp13 is expressed in clusters 1, 4, and 8; (C) Ndnf is expressed in clusters 0, 1, 2, 4, 5, and 8; (D) Acan is expressed in cluster 4; (E) Aldh1a2 is expressed in cluster 5; (F) Scara5 is expressed in cluster 6; (G) Top2a is expressed in cluster 8. (H-N) HCR RNA FISH for fibroblast cluster markers. Mid-blastemal region of 11dpa section probed for (H) Mmp13, (I) Acan, (J) Scara5, (K) Aldh1a2, or co-probed for (L) Acan (green) and Scara5 (red), (M) Acan (green) and Mmp13 (red), or (N) Scara5 (green) and Aldh1a2 (red). DAPI is shown in white. Scale bars = 20μm. Asterisks (*) denote blood vessel/RBC autofluorescence. Dashed boxes show magnified panels on right.

Figure 3. Integrated analysis of single cell populations through a regenerative time course

All analyses use combined and normalized 11dpa, 12dpa, 14dpa, 17dpa, and unamputated (UA) scRNAseq data sets. (A) UMAP plot of integrated data sets colored by regenerative stage: 11dpa (red), 12dpa (orange), 14dpa (yellow), 17dpa (light blue), and unamputated (blue). (B) UMAP plot of integrated data sets showing clusters and cluster cell type annotations. Assigned cell types are: fibroblasts (FB; clusters 0–5, and 13), vascular smooth muscle cells (VSM; cluster 6), epithelial cells (Epi; clusters 8, 15, and 16), macrophages (Mф; clusters 7, 11, and 19), endothelial cells (Endo; clusters 9 and 18), bone (cluster 10), monocytes (Mono; cluster 12), Schwann cells (SC; clusters 14 and 20), T cells (cluster 17), pre-osteoclasts (PreOC; cluster 21), and neutrophils (N; cluster 22). (C) UMAP plot of integrated data set (gray) showing the cluster distribution of cells from each regenerative stage (pink). (D) The percentage of total cells represented by each cluster for the given regenerative stage. Each stage has been compared to the proportion of cells in UA, and significant changes were determined by differential proportion analysis (marked with asterisk); all p-values reported in Table S4. Clusters are categorized by overarching cell types (fibroblast or bone, immune, vasculature, or neural). Significance values are as follows: * denotes p<0.05, ** denotes p<0.01, *** denotes p<0.001.

Figure 4. Vasculature differentiation trajectory of integrated data set

Computationally defined SPRING trajectory analysis of cells from the integrated data set vascular clusters 6, 9, and 18. (A) Force-directed plot of cells showing clusters of VSMs, vascular endothelial cells, and lymphatic endothelial cells. (B) SPRING plot as in (A) with regenerative stages of each cell colored coded: 11dpa (orange), 12dpa (light blue), 14dpa (purple), 17dpa (dark blue), unamputated (yellow). Differential clustering of blastema cells and unamputated cells suggests tissue specific differentiation of VSMs and the vascular endothelium (curved arrows). (C) Gene expression overlay on VSMs. Rgs5 is expressed in all cells, Gadd45b is more highly expressed in UA cells, and Lgals1 is more highly expressed in blastema cells. High expression is in green and low expression is black. (C’) UMAP plots of RNA velocity and expression data for individual genes Gadd45b and Lgals1 at 14dpa and 17dpa; total RNA velocity stream plot for each stage in gray (right). (D) Gene expression overlay on vascular endothelial cells. Pecam1 is expressed in all cells, Rnd1 is more highly expressed in UA cells, and Egfl7 is more highly expressed in blastema cells. (D’) UMAP plots of RNA velocity and expression data for individual genes Rnd1 and Egfl7 at 14dpa and 17dpa; total RNA velocity stream plot for each stage in gray (right).

Figure 5. Fibroblast differentiation trajectory of integrated data set

(A) UMAP plot of unbiased re-clustering of fibroblast and bone cells from integrated data set (Figure 3: clusters 0–5, 10, and 13), reveals 14 refined clusters. (B) Computationally defined SPRING trajectory analysis of cells from the integrated data set showing fibroblasts (FB) or bone (B) are not predicted to transdifferentiate into Schwann cells (SC), monocytes (M), macrophages (MΦ), pre-osteoclasts (pOC), endothelium (Endo), epithelium (Epi), T cells (T), or vascular smooth muscle (VSM). Fibroblast SPRING lineage trajectory overlaid with (C) predicted bone lineage from cluster 1 to cluster 9 (curved arrow). Marker gene expression for each cluster shown with Bglap, Ibsp, and Postn. High expression is in green and low expression is black. (D) Computationally proposed mesenchymal stem cell lineage from cluster 2 to clusters 4, 6, 8, 12, and 13 (curved arrows), with distinct lineages marked by Tnmd, S100a4, and Smoc2. Curved line depicts proposed terminally differentiated cells. (E) UMAP plots of total RNA velocity stream data for 11dpa, 12dpa, 14dpa, and 17dpa individual stage datasets; the position of bone cells in each plot is marked with a ‘B’. (F) Cluster 10 marks mitotic cells (majority highlighted by black circle) and (G) clusters 0, 3, and 5 may not contribute to a lineage, but are (H) enriched for early stage blastema cells (arrow pointing to orange).

Figure 6. Analysis of blastema fibroblast population dynamics

Differential proportion analysis of fibroblast clusters parsed by regeneration profile where clusters in (A) have no significant population dynamics between blastema stages and unamputated. (B) Cells in clusters 2, 4, 6, and 12 are significantly depleted as compared to unamputated and (C) cells in clusters 0, 3, 10, and 11 are enriched during regeneration as compared to unamputated. Significance values are as follows: * denotes p<0.05, ** denotes p<0.01, *** denotes p<0.001. All p-values are presented in Table S6. (D) Subset of genes enriched in blastema stages as compared to unamputated. Gray depicts low expression and dark purple is high expression; small circles depict a low percentage of cells and large circles depict a high percentage. (E) Violin plots of genes in (D) at four blastema stages as compared to unamputated. Black points represent individual cells and the colored curve shows the distribution of cells at a given expression level.

Figure 7. Mest expression during digit tip regeneration

(A-E) UMAP feature plots of Mest expression in all digit tip cells in UA, 11dpa, 12dpa, 14dpa, and 17dpa, respectively. Purple marks high expression and gray marks low expression; plots refer to Figures S5, S4, S3, S2, and 1, respectively. (F-J) DIG-labeled RNA section in situ hybridization for Mest on unamputated and regenerating digit tips. (F) Unamputated digit tip with orientation shown by schematic above. (G) 11dpa with region of the blastema depicted in the schematic above. Additional regenerative stages include (H) 12dpa, (I) 14dpa, and (J) 17dpa. Asterisks (*) denote artifacts from coverslipping. Abbreviations: (N) nail, (CT) connective tissue, (E) epithelium, (B) bone, (BL) blastema. Scale bar = 100μm. (K-O) HCR RNA FISH for Mest. (K) Mesenchymal region in the unamputated digit tip. (L-O) Mid blastemal region in (L) 11dpa, (M) 12dpa, (N) 14dpa, and (O) 17dpa blastemas. Arrows denote examples of Mest positive cells. Asterisks (*) denote blood vessel/RBC autofluorescence. Scale bar = 20μm.