Spatiotemporal transcriptomic mapping of regenerative inflammation in skeletal muscle reveals a dynamic multilayered tissue architecture

Abstract

Tissue regeneration is orchestrated by macrophages that clear damaged cells and promote regenerative inflammation. How macrophages spatially adapt and diversify their functions to support the architectural requirements of actively regenerating tissue remains unknown. In this study, we reconstructed the dynamic trajectories of myeloid cells isolated from acutely injured and early stage dystrophic muscles. We identified divergent subsets of monocytes/macrophages and DCs and validated markers (e.g., glycoprotein NMB [GPNMB]) and transcriptional regulators associated with defined functional states. In dystrophic muscle, specialized repair-associated subsets exhibited distinct macrophage diversity and reduced DC heterogeneity. Integrating spatial transcriptomics analyses with immunofluorescence uncovered the ordered distribution of subpopulations and multilayered regenerative inflammation zones (RIZs) where distinct macrophage subsets are organized in functional zones around damaged myofibers supporting all phases of regeneration. Importantly, intermittent glucocorticoid treatment disrupted the RIZs. Our findings suggest that macrophage subtypes mediated the development of the highly ordered architecture of regenerative tissues, unveiling the principles of the structured yet dynamic nature of regenerative inflammation supporting effective tissue repair.

Keywords: Expression profiling; Inflammation; Macrophages; Skeletal muscle.

Figure 1. Spatial and single-cell transcriptomics integration and enhanced-resolution clustering resolve the cellular distribution of myeloid subtypes during regenerative inflammation.

(A) H&E of TAs at 4 days after CTX injury used for ST. Insets indicate histopathological annotations (red-C1: regenerative muscle; blue-C2: necrotic/inflammatory lesions; green-C3: healthy muscle). C2 shows segmental necrosis of pale fibers with loss of cytoplasmic structures, active phagocytosis, C or delta lesions, and membrane damage (95), compared with C1, which includes inflammatory cells and regenerating myocytes. Scale bars: 50 μm. (B) Spatial clustering (Leiden algorithm; cluster resolution 0.3) identifies 3 discrete regions overlapping with histopathology annotations of A. The number of spots (n = 767) is indicated. (C) Enhanced subspot resolution clustering (BayesSpace) identified 7 spatial domains, not resolved at spot-level clustering. (D) Spatial expression patterns at subspot resolution of genes defining the myeloid subsets characterized by scRNA-Seq (7). Color scale shows log-normalized counts for each subspot. Gene label color corresponds to the classification in left panel of H. (E) Heatmap of the spatial expression of top predicted and curated markers, highlighting the specificity of the spatial BayesSpace clusters. Gpnmb is highlighted. (F) Identification of tissue compartments using NMF-based decomposition and day 4 after CTX reference immune subtype signatures (7). Spatial plots show cell abundance. (G) Dot plot of the estimated NMF weights of subtypes across 7 predicted NMF components. Note the differential abundance of MacII and MacIII subtypes and the overlap of DC subsets with MacIV. MacI, resolution-related MFs; MacII, GFEMs; MacIII, infiltrating monocytes/proinflammatory MFs; MacIV, antigen-presenting MFs (7). (H) Distribution and estimated cell abundance of MF and DC subtypes associated with specific NMF cellular compartments. Insets: histological area on NMF3, predicting MacIII and the formation of LGCs (encircled). The local spatial expression of known markers of LGCs (Ccl2, Ccl7) (58, 86) overlaps with the histological features and other MacIII markers (Ly6c2; D, Plac8; Supplemental Figure 1E). Scale bars: 500 μm. (I) IF detection of LGCs (MacIII) and GFEMs (MacII) by GPNMB, Ly6C, and F4/80 (green) in C57BL/6J animals at day 4 after CTX injury. Split channels are shown. White boxes indicate 2 LGC-like structures. Scale bars: 100 μm. (J) Upper: IF of GPNMB⁺ MFs and eMyHC⁺ fibers at day 4 after CTX injury. Lower: high-resolution volume projection confocal image of GPNMB⁺ MFs and eMyHC⁺ regenerating fibers preferential spatial proximity (3D reconstruction distances are indicated). Scale bar: 100 μm.

Figure 2. The sequential appearance of specialized MF subtypes orchestrates skeletal muscle regeneration.

(A) Analysis workflow for CD45⁺ cell scRNA-Seq and ST (Visium) of regenerating and dystrophic muscle. Cell suspensions were collected from digested TAs of adult mice at 1, 2, and 4 days after CTX injury and steady-state GAST of 2-mo D2.mdx. PBMC datasets from noninjured C57BL/6J mice were from 10x Genomics. Enhanced spatial resolution, deconvolution, and cooccurrence of myeloid subtypes are achieved by single-cell and spatial dataset integration (BayesSpace and Cell2location). (B) Single-cell transcriptomes from CD45⁺ cells at days 1, 2, and 4 after CTX injury, and PBMCs were harmony integrated and batch-effect corrected. Data (24,382 cells) are presented as a PaCMAP projection and color coded by origin. (C) Integrated transcriptomic atlas of 10 major populations (SingleR automated cell type annotation ImmGen database) (96). Cell types are color coded. Right: cell-type proportions and compositional dynamics. MFs account for 41.2% of all immune cells. (D) Cells in the macro-clusters of interest (monocytes, MFs, and DCs) were reanalyzed in isolation. t-distributed stochastic neighbor embeddingg (t-SNE) visualization reveals local differences. Cells are colored by major cell-type classification. (E) Clustering of the isolated cell types from D resolved 10 subtypes of monocytes, MFs, and DCs. Subcluster composition (absolute numbers) is presented as an alluvial plot. (F) t-SNE visualization of subtypes of monocytes, MFs, and DCs. (G) Dot plot of top DEGs distinguishing the monocyte/MF/DC clusters (3 DC: clusters 7, 9 and 10; 2 monocytic: clusters 3 and 6; 2 MF subtypes: clusters 1, 4; and 3 MF transitional states: clusters 2, 5, and 8; A). Dot size represents the percentage of cells expressing each marker within a cluster. Gpnmb, the top GFEM marker (cluster 1) is highlighted in red. (H) t-SNE colored by cluster and inferred pseudotime (Slingshot; principal curves are smoothed representations of each lineage) with 4 predicted cell fates: 1 monocyte (patrolling monocytes), 1 MF (GFEMs), and 2 DC lineages. Origin determines the circulating Ly6C^hi monocyte population, projected at the start of all trajectories. Trajectory 1 predicts the patrolling monocyte differentiation (not relevant in injury) (37). (I) Gene expression dynamics of monocyte/MF/DC subpopulations resolved along latent time. Cells were subjected to trajectory inference using Monocle’s (97) differential expression analysis to identify lineages. Top likelihood-ranked genes by branch and pseudotime are shown. Gpnmb is highlighted.

Figure 3. Infiltrating myeloid cell transcriptomic profile comparison from healthy, acutely injured, and early stage dystrophic muscle.

(A) Module score of gene sets representing functional markers of each myeloid subset visualized with a t-SNE in the CTX injury time-course scRNA-Seq dataset. Genes for each module are indicated. (B) Module score of gene sets representing functional markers of each myeloid subset visualized with a t-SNE in the 2-mo D2.mdx scRNA-Seq dataset. Genes for each module are indicated. (C) Composition bar plot of the major immune cell types in the Harmony integrated dataset (2-mo D2.mdx + CTX injury). (D) Heatmap of top genes in 2-mo D2.mdx vs. CTX-injury monocytes, MF, and DC subsets. (E) Pairwise Spearman’s correlation plot of monocytes, MF, and DC subsets identified in 2-mo D2.mdx, CTX-injury, and resident muscle MFs from healthy quadriceps (36) (1,300 total; presto Wilcox AUC; logFC > 0.5, P-adj < 0.1, AUC > 0.5). The resident MFs are represented according to the following nomenclature chosen by the authors (36): cluster 0 (cluster 0), Cd209 cluster (cluster 1), Ccr2 cluster (cluster 2), and proliferating cluster (cluster 3). Color intensity and circle size are proportional to the correlation coefficients. Note (a) the relative uniqueness of resident MFs, (b) the high correlation of proinflammatory monocytes (cluster 3) in both CTX and D2.mdx datasets, and (c) the high correlation of cycling MFs in D2.mdx (cluster 6) with the resident MFs (proliferating cluster). (F) Heatmap of known regulators of regeneration in 2-mo D2.mdx versus CTX-injury monocyte, MF, and DC subsets.

Figure 4. GPNMB is a marker and component of GFEMs, and its deficiency impairs regeneration.

(A) Upper: Brightfield images of day 4 after CTX Ly6C^loF4/80^hiGPNMB^– and Ly6C^loF4/80^hiGPNMB⁺ muscle-infiltrating MFs ex vivo after 12 hours in culture (equal number of cells were seeded). Lower: apoptosis was assessed by cleaved caspase-3 immunostaining. Scale bars: 100 μm. (B) Percentage of apoptotic (CASP3⁺) Ly6C^loF4/80^hiGPNMB⁺ and Ly6C^loF4/80^hiGPNMB^– MFs (unpaired t test with a P = 0.0007; n = 3). (C) Effect of GPNMB^– and GPNMB⁺ MF-derived conditioned media on the proliferation and differentiation of C2C12 myoblasts. Scale bars: 100 μm. (D) Effect of GPNMB^– and GPNMB⁺ MF-derived conditioned media on C2C12 myoblasts (n = 3). Proliferation index as percentage of Ki67⁺ cells (n = 10 fields/group; unpaired t test, P = 0.0008). Fusion index as percentage of myotubes (visualized by heavy chain of myosin II) with more than 3 nuclei (n = 10 fields/experiment/group; unpaired t test, P = 0.002). (E) Spatial cell proximity quantification of GFEMs (CD68⁺GPNMB⁺) versus other MF subtypes to regenerating eMyHC⁺ fibers at day 4 after CTX in C57BL/6J animals (n = 3; >200 mm² tissue area/sample; unpaired t test, P = 0.0267). (F) Detection of GFEMs by CD68 and GPNMB, in relation to eMyHC⁺ fibers in C57BL/6J animals at day 4 after CTX injury. Insets indicate the split channels. Scale bars: 100 μm (left); 20 μm (insets). Lower panel indicates high-resolution volume projection confocal images (3D reconstruction distances are shown). (G) Left: H&E images of D2.Gpnmb^– (KO) and D2.Gpnmb⁺ (WT) TAs at day 8 after CTX injury. Note the near complete absence of regenerating fibers (highlighted in white) and extensive inflammation areas (highlighted in red) in the D2.Gpnmb^–. Right: IF detection of newly formed fibers by eMyHC in D2.Gpnmb^– and D2.Gpnmb⁺ TAs at day 8 after CTX injury. Scale bars: 1 mm (left H&E); 500 μm (right H&E); 100 μm (IF panels). (H) IF detection of mature MFs by CD68 in KO and WT at day 8 after CTX injury, correlating to the extent of unresolved inflammation. Scale bars: 100 μm (main); 1 mm (insets). (I) H&E images of regenerating TAs (day 8 after CTX injury) from WT (C57BL/6J) and D2.Gpnmb^– animals used for ST. Insets show magnified H&E areas (green rectangles). Scale bars: 500 μm (middle panels); 50 μm (far left and right panels). (J) Enhanced subspot resolution clustering of regenerating TAs (day 8 after CTX injury) from WT and D2.Gpnmb^– identified 5 spatial clusters (n spots/group are indicated). Scale bars: 500 μm. (K) Top marker gene expression after z score transformation for each spatial cluster. Dot size represents the percentage of subspots expressing the gene. (L) Spatial expression of representative healthy muscle (Myh4), differentiating myoblasts (Myog), and persistent inflammation/mature MF (Cd68) genes in WT and D2.Gpnmb^–. Note the loss of the distinct structure of regenerative zones in the KO. Scale bars: 500 μm. (M) Gpnmb and Myh3 spatial expression patterns in the C57BL/6J day 8 after CTX ST sample confirm the proximity of GPNMB⁺ MFs to early stage regenerating fibers and distinct tissue organization around lesions. Scale bars 500 μm (upper); 50 μm (lower). H&E previously presented in I. This duplication is intended to provide the location and context for the presented magnified feature plots. In all bar graphs, bars represent mean ± SD. *P < 0.05; **P < 0.01; ***P < 0.001.

Figure 5. scRNA-Seq and ST integration with enhanced-resolution clustering resolve the complex dystrophic muscle architecture and cellular distribution and profiles of myeloid subtypes.

(A) Left: H&E images of mouse GAST from 2-mo D2.mdx used for ST. Histopathological annotation areas are noted: regenerative muscle (yellow), necrotic/inflammatory lesions (green), healthy muscle (blue). Right: Percentage of spots in the annotated areas. (B) Enhanced subspot resolution clustering (BayesSpace) identified 7 spatial clusters (color coded), which were not resolved by pathologist annotations. The white rectangle highlights a lesion with structured inflammation and regeneration zones. (C) Top marker gene expression after z score transformation for each spatial cluster. Dot size represents the percentage of subspots expressing the gene. (D) Spatial expression of representative genes coding for markers of each spatial cluster: Ccl2 and Ccl7 (LGCs: cluster 2), Itgax, Mmp12, Trem2, and Gpnmb (resolution-related MFs and GFEMs: cluster 6), Myog and Myh3 (newly regenerating fibers: cluster 1), Pvalb and Tpm3 (healthy muscle: cluster 7), Col1a1 (ECM: cluster 4), and Esam (endothelial cell/vasculature-enriched areas: cluster 3). Note the differential spatial expression patterns in the highlighted region of B. (E) Single-cell transcriptomes derived from CD45⁺ cells from 2-mo D2.mdx GAST. A total number of 4,811 myeloid cells (MFs, monocytes, DCs; SingleR automated annotation using the ImmGen database) were analyzed. Data are presented as a t-SNE projection to visualize variation in single-cell transcriptomes. The subsampling-based clustering approach (chooseR) resolved 7 myeloid subsets (color coded). (F) Top marker genes for the 7 identified clusters. Dot size represents the percentage of expressing cells within a group, and color scale represents the average expression level (row z score) across all cells within the cluster. (G) Left: 2D embeddings visualizing cell cycle phases of the 2-mo D2.mdx scRNA-Seq dataset using t-SNE. NA indicates the number of unassigned cells. Right: 2D embeddings visualizing 3 subclusters of cycling cells from parent cluster 6 using VeloViz (98) embeddings. Cell numbers for each subcluster are indicated. Arrows show velocity projections (velocyto.R). (H) Identification of tissue compartments using NMF-based decomposition and 2-mo D2.mdx reference immune subtype expression signatures (7). Heatmap of the estimated NMF weights of subtypes (rows) across 6 predicted NMF components (columns), corresponding to the identified cellular compartments. Relative weights normalized across domains for every MF subtype are shown. (I) Spatial plots show cell abundance for each immune cell subtype calculated in H. Scale bars: 500 μm.

Figure 6. ST organization of dystrophic muscle upon intermittent Pred treatment.

(A) Left: H&E images of GAST from 2-mo D2.mdx mice treated weekly (Q.W.) for 4 weeks with Pred, used for ST. Histopathological annotation areas: regenerative muscle (yellow), necrotic/inflammatory lesions (green), healthy muscle (blue). Right: Percentage of spots in annotated areas. Each section is from a different biological replicate, and each library was obtained from a separate Visium experiment followed by bioinformatic integration to remove batch effects. (B) Enhanced subspot resolution clustering (BayesSpace) identified 7 meaningful spatial clusters (color coded) unresolved by pathologist annotations. (C) Comparison of the spatial cluster (color coded) subspot composition in untreated D2.mdx and D2.mdx treated Q.W. with Pred. UNT, untreated. (D) Spatial expression of representative glucocorticoid receptor (GR) targets, regenerative muscle, inflammation, and atrophy marker genes is shown. Myh3 indicates newly regenerating fibers, Tsc22d3 GR target engagement, Lyz2 inflammatory myeloid cells, and Trim63 atrophy-inducing pathways. (E) Representative DEGs in untreated D2.mdx versus D2.mdx+Q.W. Pred comparison, grouped in functional categories. Dot size represents the percentage of spots within a treatment group. (F) Unbiased global tissue-wide DGE of all spots in D2.mdx+Q.W. Pred versus D2.mdx-UNT (red dots indicate significant DEGs; P < 0.05, logFC > 1). Top DEG names are indicated. (G) GO pathway enrichment analysis of the DEGs in D2.mdx+Q.W. Pred versus D2.mdx-UNT ST datasets. Top significant up- and downregulated pathways are shown (P < 0.001, fold enrichment > 2). Gray box, enriched GO terms with Pred treatment; red box, downregulated terms with Pred-treatment. (H) Identification of tissue compartments in D2.mdx+Q.W. Pred-treated ST samples using NMF-based decomposition and 2-mo D2.mdx reference immune subtype expression signatures (7). Dot plot of estimated NMF weights of cell subtypes (rows) across 6 predicted NMF components (columns) corresponding to the identified cellular domains. Relative weights normalized across components for every MF subtype are shown. (I) Spatial plots show cell abundance for each immune cell type calculated in H. Scale bars: 1 mm.

Figure 7. RIZs are disrupted by GC treatment in early stage dystrophy.

(A) Magnified view of representative structures and RIZs in untreated 2-mo D2.mdx muscles. Zone A represents the center of an inflammatory lesion occupied by LGCs (Ccl7⁺); zone B is occupied by a gradient of resolution-related MFs (Mmp12⁺) and GFEMs; and zone C represents the regeneration zone marked by developing myofibers (Myh3⁺). The correlation of observed/expected zone organization is quantified per subspot. H&E data have been previously presented in Figure 4A and provide the location and context for the magnified feature plots. (B) Abnormal tissue zones in 2-mo D2.mdx+Q.W. Pred animals. Note the disintegration/absence of regenerating fibers (Myh3⁺) in zone C. H&E images have been previously presented in A and provide the location and context for the magnified feature plots. (C) Global subspot correlation (Spearman’s) of spatial gene expression in 2-mo D2.mdx-UNT and 2-mo D2.mdx+Q.W. Pred samples. Higher correlation in Pred samples suggests a collapse and rearrangement of inflammatory (zones A, B) and regenerative (zone C) zones (color coded). (D) Example of RIZs formed with alternative markers. Scale bars: 500 μm (left panel,A, H&E); 1 mm (left panel, B, H&E); 100 μm (right panel, A); 50 μm (right panel, B); 100 μm (D). (E) Representative H&E region of RIZs in 2-mo D2.mdx GAST validated by IF. The MF subtypes and zones were visualized with IF (bottom panel shows the absorbed signal) for CCL2 (zone A), MMP12 (zone B), and eMyHC (zone C). Dotted lines indicate the zones and interface layer (red: necrotic lesion) selected for cell density quantification in F. Scale bars: 1 mm (lower left); 100 μm (others). (F) Stacked bar histogram of CCL2⁺, MMP12⁺, and eMyHC⁺ cell density inside (–1 to –100 μm) and outside (+1 to +440 μm) the necrotic boundary in E. (G) Representative IF region with 2 inflammatory lesions in 2-mo D2.mdx+Q.W. Pred GAST samples. MF subtypes and regenerating fibers were visualized as in E. Bottom panel indicates the cell density and distribution. Scale bar: 100 μm.

Figure 8. ATF3 directly regulates a GFEM-like transcriptional program.

(A) De novo motif enrichments around the TSSs of GFEM-associated genes. ATAC-Seq peaks of day 4 after CTX Ly6C^lo repair muscle MFs within 50 kb of TSSs were selected as input. Detected motif matrices, P values, and background are shown. (B) Predicted scores and motif map of 2 distal enhancers (E1, E2) and 1 proximal (P) site around the Gpnmb locus, selected based on ATAC-Seq. Open and closed circles indicate the absence or presence of corresponding TF mRNA, respectively. Motifs of ATF3 and JUN are highlighted. (C) Genome browser view of the Gpnmb locus depicting capture Hi-C (in naive BMDMs), ATAC-Seq (blood monocyte and muscle-infiltrating MFs; normalized scale), and ChIP-Seq (in naive BMDMs and muscle-infiltrating MFs) for indicated TFs, active transcription histone marks (H3K27Ac), and elongating polymerase II (S2P). CTCF and RAD21-defined transcriptional unit boundaries, distal (E1, green; E2, blue) and proximal (P, red) Gpnmb-associated regulatory elements and track scales are indicated. (D) Heatmap of the highest expressed TFs (decile filtered) in the myeloid subtypes of the CTX scRNA-Seq dataset (Figure 2F). Hierarchical clustering and average log-normalized expression values are shown. Atf3 is highlighted. (E) Spatial expression feature plots of top TFs with detected binding in the regulatory elements in the D2.mdx samples. (F) Magnified view of Atf3 spatial expression in representative RIZs in untreated 2-mo D2.mdx muscles. The correlation of observed/expected zone organization is quantified per subspot for each zone and indicates an overlap of Atf3 with zone B. Data (left and middle) have been previously presented in Figure 4A and Figure 7A, respectively, and provide the location and context for the magnified feature plot and expected spatial organization. Scale bars: 500 μm (E, left panel, F). (G) IF region of a lesion in 2-mo D2.mdx GAST muscle. MF subtypes were visualized with CCL2 (red, zone A) and ATF3 (yellow, zone B), and regenerating fibers with eMyHC (green, zone C). Scale bar: 100 μm. (H) Volcano plot showing the DEGs in the Atf3^–/– naive BMDMs (P < 0.01, FDR < 0.01). Number of DEGs and gene labels of GFEM-predicted markers among top DEGs are shown. (I) Atf3 mRNA expression in WT and Atf3^–/– naive BMDMs (n = 3; unpaired t test, P < 0.0001).