Overview of investigate the spatial gene expression in T10 spinal cord sections
To ensure the uniformity of the wound area, we adapt the lateral hemisection model, instead of contusion (41). Previous study revealed that the region of glial scar approximately distributed 0.5 mm away from the site of injury after SCI (42). To decipher the cell type and region of glial scar as much as possible, we harvested the spinal cord tissue 1.5 mm away from the site of T10 right lateral hemisection (figure 1A). The Visium Spatial Gene Expression Slides were used to capture and profile the mRNAs in the spots of spinal cord sectons at 3d, 7d, 14d and 28d after T10 right lateral hemisection, which intended to cover the dynamic process of glial scar formation after injury (43). Because each section contained up to 4,992 spots in its capture area and the size of the Visium Spatial Gene Expression Slides is 6.5 mm×6.5 mm (figure 1B), approximately four spinal cord tissues can be placed on the same chip to capture more spots as shown in figure S1A. Supplementary Table S1 contains a summary of the data evaluation. Overall, 15449 spots within the 32 sections were analyzed (Table S2, figure S1A and B). Data showed the median sequencing depth of single spot at approximately 219,442 Unique Molecular Identifiers (UMIs) and 4,875 genes in our study (Table S2, figure S1C). The H&E staining images served as a reference for subsequent unsupervised analysis. We focused on the extent of scar foci and boundary that conduct the scar microenvironment, identification cell type, the reconstruction of cell trajectories, cell-cell interactions within the scar formation (figure 1B).
Defining the spatiotemporal boundary of scar in the injury site
Integrated clustering analysis of ST from 32 samples using Seurat revealed 4 distinct clusters when visualized on a UMAP plot (figure 1C). These 4 clusters covered all cell types in the SCI site and adjacent uninjured spinal cord. The numbers of UMIs and genes in cluster 1, cluster 2 and cluster 3 were larger than that in clusters 0 (figures 1D and 1E). To better understand the spatial distribution of these cell types, we compared the H&E staining images (figure 1F) with their ST data (figure 1G) and mapped the clusters back to their spatial location (figure 1H). These identified spatially patterns could match the changes in pathological morphology. Specifically, cluster 2 and cluster 3 characteristically distribute in the scar regions showed in H&E staining images, consistent with the larger UMIs numbers. Based on these findings, we inferred that cluster 2 and cluster 3 represented multiple cellular components that are known to comprise the glial scar. On the contrary, cluster 0 and cluster 1 represented cell types in the uninjured spinal cord (figure 1E). Together, we firstly identify the approximate scope of glial scar at three stages.
Meanwhile, we listed the top 5 highest differentially expressed genes (DEGs) in each cluster during this process from ST (figure 1I, figure S1D). The highest DEGs provided a unique molecular signature for each cell type. For example, cluster 0 cells highly expressed typical neuron marker Snap25 and mainly represent the neurons in the gray matter area of the uninjured spinal cord. Cluster 1 cells highly expressed myelin associated protein (such as Plp1, Mbp and Mal) and mainly represent the oligodendrocytes in the white matter area of the uninjured spinal cord.
The highest DEGs in cluster 2 were secreted phosphoprotein-1 (Spp1), glycoprotein nonmetastatic melanoma protein B (Gpnmb), galactose-specific lectin 3 (Lgals3) and cathepsin D (Ctsd). Spp1, also known as OPN (Osteopontin), is a broadly expressed pleiotropic protein (44). It has multifunction in the pathophysiology of several inflammatory, degenerative, autoimmune and oncologic diseases (45). Spp1 was firstly discovered as a component of bone matrix, but it is also expressed by subsets of neurons including alpha-retinal ganglion cells (αRGC), alpha motor neuron, and several classes of glial cells including astrocyte, Schwann cells, Müller cell and disease-associated microglia (46); OPN is an alpha motor neuron marker in the mouse spinal cord (47). Recent studies have shown that the levels of Spp1 are affected by neural injury. Spp1 exert both pro- and anti-inflammatory roles and contributes to tissue damage (the Yin) not only by recruiting harmful inflammatory cells to the site of lesion, but also increasing their survival (the Yang) (45, 48). For example, Spp1 is capable of stimulating mTOR activity (49) and promotes RGC regeneration when introduced in combination with either insulin-like growth factor 1 (IGF-1) or brain-derived neurotrophic factor (BDNF) (46). Additionally, Spp1 was reported to stimulate the outgrowth of injured motor axons but not sensory neurons (50). Our data showed the spatiotemporal distribution of Spp1 (figure S1E). Spp1 indeed appear on some motor neurons and significantly increased in injury site, consistent with previous studies (45, 47, 48). Lgals3 plays a role in numerous cellular functions including cell growth, cell adhesion, apoptosis, pre‑mRNA splicing, differentiation, transformation, angiogenesis, inflammation, T-cell regulation, host defense and fibrosis (51, 52). Based on the high expressions of stress response and inflammatory regulator, cluster 2 cells were considered as the glial and immune cells of scar. The identification of scar-associated cluster 2 serves important functions in defining the boundary of scar accurately from the cellular level.
The highest DEGs in Cluster 3 were collagen type I alpha 1 chain (Col1a1), collagen type I alpha 2 chain (Col1a2), insulin like growth factor binding protein 6 (Igfbp6) and matrix gla protein (Mgp). Col1a1 and Col1a2 belong to the type I collagen (Col I). To date, pericytes and fibroblasts have been reported to be the cell types that produce Col I after SCI (53, 54). Among the extracellular matrix (ECM), Col1a1 and Col1a2 were reported to be the most highly expressed in injured spinal cord of contusion SCI model at 14 post-SCI (dpi).
In order to further explore the function of these cell types, the Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis data are enriched for each cluster. The Go annotation matched well with the anatomical annotation (figure 1F). Specially, some significant GO such as synaptic vesicle cycle, regulation of neurotransmitter levels and ATP metabolic process are highly related to nomal activities of neuron in Cluster 0. Some significant GO such as myelination, axon ensheathment and positive regulation of neuron projection development are highly related to electrical signal transmission in Cluster 1. On the contrary, cluster 2 cells and cluster 3 cells exhibited damage associated patterns of gene expression, with a significant enrichment for genes that have functions in the extracellular structure organization, angiogenesis, wound healing, regulation of neuron death, regulation of inflammatory response and regeneration (figure 1J).
Based on the symmetric distribution pattern of cluster 0 cells and cluster 1 cells in normal spinal cord, we divided the T10 right lateral spinal cord into four layers (figure S1F) to further explain the relations between the subgroup and spatiotemporal distribution. After countering all samples, we found at least 1.2 mm distance from the scar center (approximately 12 spots) can cover the scar boundary. Subsequently, GSVA score of apoptosis, ferroptosis, GABAergic synapse and PI3K-Akt signaling pathway were respectively counted (figure S1G). We found the GSVA scores of ferroptosis and GABAergic synapse signaling pathway peaked in the scar center and were symmetrically distributed on the left and right side at three stages in Layer 3 (figures 1K-L). Layer 3 enriched by cluster 0 neurons suffered from ferroptosis after spinal cord injury, consistent with previous studies (55). Together, we firstly present a spatiotemporal atlas that systematically describes the spatial archetypes and cellular heterogeneity of glial scar at three stages in the early acute stage (3 dpi), subacute stage (7 dpi) and intermediate stage (14 dpi, 28 dpi) after spinal cord injury.
Deciphering spatiotemporal distribution of resident cells in the scar
To determine the heterogeneity within the cluster 2 and 3 which represented multiple cellular components that comprise the glial scar, re-clustering analysis was performed and visualized on UMAPs (figure 2A), which revealed twelve clusters identified by annotated lineage markers (figure 2B-C). These 12 clusters represented approximately all major cell types that are known to comprise the scar including microglia, macrophages, astrocytes, oligodendrocytes, fibroblasts and endothelial cells. From the re-clustering UMAPs at the different stages after spinal cord injury (figure 2D), we found the cell spots distribution at 14 dpi and 28 dpi were similar, and the cell spots distribution at 3 dpi were significant different from that at 14 dpi and 28 dpi. By the contrast, the cell spots distribution at 7 dpi were observed between cell 3 dpi spots and 14 dpi spots. The cell spots distribution provided compelling evidence about the division of time window-the early acute stage (3 dpi), subacute stage (7 dpi) and intermediate stage (14 dpi, 28 dpi) after spinal cord injury (1, 56-58).
To gain molecular insights into the injury responses that are mediated by these cell types, we combine the subgroups and all the re-clustered celltypes into six distinct cell clusters, microglia, macrophage, astrocyte, oligodendrocyte, fibroblast and endothelial cell (figure 2E). The transcriptional signatures were shown in figure 2F. We listed the highly expressed genes in each group in detail. For example, Decorin (Dcn) was expecially expressed on fibroblast. Plp1 was expecially expressed on oligodendrocyte. Metallothionein 3 (Mt3) was expecially expressed on astrocyte. By the contrast, Spp1 was expecially expressed on multiple cellular components (such as macrophage, fibroblast and microglia). Subsequently, these six cell clusters were are mapped bact to slices (figure 2G), and we found macrophage firstly appeared in the scar area at the early acute stage (3 dpi), and then fibroblast occupied the center of the scar at the subacute stage (7 dpi). From the center of scar to the boundary, fibroblast, microglia, astrocyte and oligodendrocyte appeared in proper order at the subacute stage (7 dpi) and intermediate stage (14 dpi, 28 dpi) after spinal cord injury. We used a novel gene set analysis method (methods) infer activity maps by countering the GSVA score of cell spots in scar area. Wound healing (59, 60), a selected GO term, was mapped the in the scar region at three stages, and was significant active (figure 2H). Taken together, our data revealed the spatiotemporal distribution of multiple celltypes in the scar region and clearly depicted the boundary of each celltype.
Identifying macrophage infiltrating into the scar by ST analysis
To determine the heterogeneity within the macrophage population, re-clustering analysis was performed and visualized on spatiotemporal maps (figures 3A and 3B). We identified three macrophage subtypes, which were labeled cluster 0 to 2. Cluster 0 was identified by high expression of lysozyme (Lyz2), a specific marker for macrophage, and thyrotropin releasing hormone (Trf). Cluster 0 comprised 45.3% of macrophage present at 3dpi (figures 3C and 3D). Cluster 1 expressed higher levels of platelet factor 4 (Pf4) (a specific marker for macrophage) and the anti-inflammatory marker Arginase 1 (Arg1) than Cluster 0 and 2 (figure 3E, figure S2A), comprised 51.8% of macrophage present at 3dpi. Cluster 2 expressed moderate level of Lyz2 and low level of Pf4, and comprised 2.9% of macrophage present at 3dpi (figure 3D, figure S2A).
Although these three subtypes did not precisely correspond to the M1/M2 division of macrophage subsets, Cluster 1 cells expressed higher levels of the Arg1 and CD14 than Cluster 0 and 2 (figure S2A). It indicates Cluster 1 cells were likely to be M2 macrophages. Interesting, ST maps show the macrophages of Cluster 1 cells were always in the central party of injury site and surrounded by the macrophages of Cluster 0. The macrophages of Cluster 2 were always in the periphery of injury site at 3 dpi (figure 3A). Subsequently, cluster 1 cells decreased dramatically and cluster 2 cells significantly increased. Interesting, cluster 2 cells interlaced with cluster 0 cells and formed a circular band of cells at 7 dpi. These three clusters cell become rare at 14 dpi and 28 dpi. In addition, ST maps show the expression of macrophage marker Pf4 was significantly upregulated at 3 dpi and 7 dpi, and then back to baseline level at intermediate stage (figure S2A). There was an obvious shift in the expression of M2 macrophage marker Arg1, which peaked at 3 dpi and then fell rapidly. This was consistent with the change of macrophage phenotype after injury (61). The Pf4+ cells resided in the scar at intermediate stage (14 dpi, 28 dpi) mainly represent M1 macrophages, which resulted in an unfriend regenerate microenvironment in scar and were also characterized by CD86 and CD32 maps (figure S2A). This is the first identification of macrophage subsets spatiotemporal distribution in the scar.
It is traditionally believed that M2 macrophages could promote CNS repair while limiting secondary inflammatory-mediated injury (61). We found the macrophages of cluster 1 cells not only decreased dramatically, but also surrounded by the macrophages of cluster 0 and cluster 2. To investigate the functions of these clusters, we listed top 10 marker genes of these three clusters (figure 3E), and enriched the GO terms base on highly expressed genes (figure S2B). The top terms showed cluster 1 cells indeed have high enrichment factor in wound healing. However, cluster 1 cells also exhibited unique patterns of gene expression, with significant enrichments for chemotaxis, regulation of inflammatory response, response to interferon-gamma and cellular response to interleukin-1. Interesting, cluster 0 cells were associated with border and exhibited significant enrichments for axonogenesis, axon guidance, axon extension and axon ensheathment. It seems they can promote axon regeneration into the scar from both sides. The impressive functions matched the unique spatiotemporal distributions of cluster 0 cells resided in the scar periphery, which allowed cluster 0 cells to have better access to axon terminal. The GSVA score of axon guidance pathway analysis also supported this view (figure 3F). In addition, the top 4 GO terms exhibited by cluster 2 cells were extracellular structure organization, extracellular matrix organization, positive regulation of cell migration and angiogenesis. Particularly, the GO term of negative regulation of immune system process matched the traditional M2 macrophages functions limiting secondary inflammatory-mediated injury (61). So, we hypothesized that cluster 2 cells represent the continuators of the cluster 1 cells.
To test the hypothesis, we conducted the trajectory analysis of these three subtypes (figure 3G-J, figure S2C). Two branches were builded (figure 3G) and then the pseudo-time trajectory algorithm was applied (figure 3H). We found the evolation trajectory and transcriptional connections among the subtypes. The results indicated trajectory direction started from cluster 1 cells to a separate branch consisting of cluster 0 cells and cluster 2 cells (figure 3H-I). We observed the cell fate 1 branch enriched for genes associated with axon ensheathment, such as Sema7a, Vim, Mal, Cntn2 and Fgf1. While the cell fate 2 branch showed increased expression of extracellular matrix organization and negative regulation of immune system process related genes, such as Gpnmb, Mmp12, Col3a1, Col1a1 and Col1a2 (figure 3J).
The differentially expressed genes (DEGs) provided a molecular signature for each cell type, which were different from canonical markers (Table S3). For example, the highest DEGs in Cluster 1 were non-canonical genes such as Thbs1 (Thrombospondin 1) and Plin2 (Perilipin 2). Microvascular dysfunction is a critical pathology that underlies the evolution of secondary injury mechanisms after traumatic SCI. Thbs1 signaling was implicated in the acute neuropathological events that occur in spinal cord microvascular ECs (smvECs) following SCI (62). Additionally, unsatisfied intrinsic neurite growth capacity results in significant obstacles for injured spinal cord repair. Thbs1 was reported to be a neurite outgrowth-promoting molecule. Overexpression of Thbs1 in bone marrow mesenchymal stem cells (BMSCs) can promote neurite outgrowth in vitro and in vivo of Oxygen-Glucose Deprivation (OGD) treated motor neurons and SCI rat models. Specially, the transplantation of BMSCs+ Thbs1 could promote neurite outgrowth and functional recovery after SCI partly through the TGF-β1/p-Samd2 pathway(63). ST maps show the expression of Thbs1 was significantly upregulated at 3 dpi and 7 dpi, and then back to baseline level (figure 3K, figure S2C). Immunostaining showed the majority of Thbs1+ cells are P4/60 (a macrophage marker) positive cells at 3 dpi and 7 dpi (figure 3M). Although Spp1 was highest DEGs in Cluster 1, Spp1 expressed in multiple cell types and was unfit for non-canonical macrophage markers (figure S1E).
Another unexpected finding is that the macrophage in Cluster 2 expressed higher levels of Col I, including Col1a2, Col1a1, Col3a1, which is different from the previous study that the fibroblasts and pericytes produced Col I after SCI. Col I was highly expressed in the spinal cord, which is consistent with previous findings. Specially, Col I induced astrocytic scar formation via the integrin-N-cadherin pathway during the scar-forming phase (64). This was consistent with GO terms associated with extracellular matrix organization. ST maps show that the expression of Col1a2 was rare in the intact spinal cord tissues, and obviously upregulated after injury. It peaked at 7 dpi and was maintained at high level for 28 days (figure 3L). Immunostaining showed the majority of Col1a2+ cells are CD68 (a macrophage marker) positive cells at 3 dpi (figure 3N). Additionally, Msr1 (macrophage scavenger receptor 1) has been implicated in many macrophage-associated physiological and pathological processes including atherosclerosis, Alzheimer's disease, and host defense. In the scar, Msr1 promote the formation of foam macrophages and neuronal apoptosis after spinal cord injury (65). ST maps showed that Msr1+ cells significantly increased at 3 dpi and 7 dpi in the injured site (figure S2D). Activation of the innate immune system promotes regenerative neurogenesis in zebrafish. TNFa from pro-regenerative macrophages induces Tnfrsf1a-mediated AP-1 activity in progenitors to increase regeneration-promoting expression of hdac1 and neurogenesis. But it is unknown whether the mammals retain the mechanism. Our ST maps showed that at least Tnfrsf1a was widespread at 3 dpi and 7 dpi in the scar (figure S2E) (66).
In terms of macrophage origin, the traditional view believed that bone marrow-derived mononuclear macrophages (also known as recruited macrophages) are recruited into damage tissue via chemokine (C-C motif) receptor 2 (Ccr2) in response to chemokines and inflammatory signals and affect wound healing by secreting inflammatory factors such as TNFα (61) (67). Ccr2 encodes a receptor for monocyte chemoattractant protein-1, a chemokine that specifically mediates monocyte chemotaxis, and is expressed on monocytes and macrophages, but not in microglia or tissue resident macrophage at the resting state (68). Recent study found the mouse meninges hosted a rich repertoire of immune cells mediating CNS immune surveillance and supplied not from the blood but by adjacent skull and vertebral bone marrow. Under spinal cord injury, local bone marrow–derived monocytes can infiltrate spinal cord from the dural meninges (69). We want to distinguish position and proportion of the macrophages from the blood and tissue-resident macrophages from adjacent the dural meninges. ST maps showed that Ccr2+ monocyte-derived macrophages were present in the scar, but their proportion in the injured site was low (figures 3O and 3P), consistent with previous study. Violin Plots showed that the expression of Ccr2 was low in these three clusters (figure 3Q). Last, we profiled the combination of the activated macrophage markers to distinguish these clusters, and found cluster 1 cells expressed higher levels of heme oxygenase 1 (Hmox1) at 3 dpi. In summary, our data reveals the spatiotemporal distribution of macrophage subtypes at the injury site after SCI.
ST analysis reveals spatiotemporal changes in fibroblast subtypes
Fibroblasts are the main producers of the matrix (including ECM) and constitute the basic framework of tissues and organs (70). To assess the cellular heterogeneity among fibroblast at the injury site, re-clustering analysis was performed and visualized on spatiotemporal maps (figure 4A) and a UMAP (figure 4B-C). We identified three fibroblast subtypes, which were labeled cluster 0 to 2. Base on the spatiotemporal maps, we found cluster 2 mainly appeared in the scar at 7 dpi. Cluster 1 appeared in the scar at 7 dpi and locate in the center of scar. On the contrary, cluster 0 located in the meninges of spinal cord. The fibroblast number in Cluster 1 reached the peak at 14 dpi (figure 4D). Cluster 0 were identified by high expression of Prostaglandin D2 Synthase (Ptgds) and insulin like growth factor binding protein 6 (Igfbp6) (figure 4E). Ptgds preferentially expressed in brain and functions as a neuromodulator as well as a trophic factor in CNS. In the chronic constriction injury (CCI) of sciatic nerve model, Ptgds was overexpressed in the lumbar spinal cord, but not in the striatum (71). Igfbp6, a member of the insulin-like growth factor-binding proteins family, modulates insulin-like growth factor (IGF) activity. In spinal cord, meningeal cells, interneurons in the deep part of the dorsalhorn and around the central canal, and motoneurons were Igfbp6 positive cells (72). Igfbp6 plays a key role in neuronal apoptosis after SCI. In acute SCI model, Igfbp6 was upregulated significantly and was co-localized with active caspase-3 and p53 in neurons. When Igfbp6 was knocked down, the protein levels of active caspase-3 and Bax as well as the number of apoptotic primary neurons were significantly decreased (73). ST maps showed Igfbp6 and Ptgds were high expression levels in Cluster 0 cells (figure S3A).
To investigate the functions of these clusters, we listed top 5 differential genes between these three clusters (figure 4E), and enriched the GO terms base on highly expressed genes (figure S3B). The top terms showed cluster 2 cells have high enrichment factor in wound healing. In addition, cluster 2 cells also exhibited significant enrichments for angiogenesis, extracellular matrix organization, response to transforming growth factor beta and collagen metabolic process. The GSVA score of ECM-receptor interaction pathway analysis also supported this view (figure 4F). But cluster 2 cells diminished gradually. Interesting, cluster 0 cells were associated with meninges and exhibited significant enrichments for axonogenesis, axon development, axon guidance, axon ensheathment and regulation of epithelial cell proliferation. What is more, the top 10 GO terms exhibited by cluster 1 cells were leukocyte chemotaxis, phagocytosis, regulation of immune effector process, T cell activation and activation of immune response. cluster 1 cells gradually become the main component in the scar and their enrichments were associated with inflammatory hyperactivation in the intermediate stage after spinal cord injury, consistent with previous studies (74). The GSVA score of glycosaminoglycan degradation in top 10 pathway analysis enriched in cluster 1 cells significantly increased from 7 dpi to 28 dpi (figure 4F), which may be related to the maturation and contracture of the scar, resulting in compartmental structures containing multiple cells (figure S3C).
To verify the evolutionary relationship between these three cell populations, we conducted the trajectory analysis (figure 4G-I). Two branches were found (figure 4G) and then the pseudo-time trajectory algorithm was applied (figure 4H). The trajectory direction started from cluster 2 cells and cluster 0 cells to a separate branch consisting of cluster 0 (figure 4H-I). We observed the cell fate 1 branch showed increased expression of extracellular matrix organization, such as Dcn, Col1a2 and Bmp7. While the cell fate 2 branch enriched for genes associated with axon guidance and axonogenesis, such as Sema3b, Sema5a, Periostin (Postn), Cntf, Ogn, Mmp15 and Fgf1 (figure 4I). Interesting, the origin of cluster 1 cells came from three directions. That is undoubtedly true that a certain percentage of cluster 1 cells came from adjacent the dural meninges after spinal cord injury. Our data showed cluster 2 cells appeared at 3 dpi and could transform into cluster 1 cells. Where these cells originate from? Previous studies reported that in mouse models that develop fibrotic tissue, the primary source of scar-forming fibroblasts was type A pericytes. Perivascular cells with a type A pericyte marker profile also exist in the human brain and spinal cord, suggesting it is conserved across diverse CNS lesions (75, 76). Attenuation of pericyte-derived fibrosis represents a promising therapeutic approach to facilitate recovery following spinal cord injury(77). We hypothesized that cluster 2 cells originate from type A pericyte and then tested the expression of type A pericyte marker GLAST (gene name Slc1a3). We verified the expression level of GLAST and GLAST+ cells proportion at the three clusters, and found GLAST+ cells were the primary origin of fibroblast in the scar (figure 4J and K), consistent with previous study(76).
Homeostatic fibroblasts of cluster 0 were identified based on its spatiotemporal distribution and the expression of several annotated markers of steady-state fibroblast, such as P4hb and Gsn. Homeostatic fibroblasts of cluster 0 were the predominant subtype in the uninjured spinal cord, but by 3 dpi they were replaced by the damaged-associated fibroblasts (DAF) of cluster 1 and cluster 2 (figure 4D, figure S3A). DAF were always in the center of scar. DAF were identified by high expression of Gpnmb, Ftl1, Col3a1 and Postn. Myofibroblast differentiation marker Postn is an important extracellular matrix protein and coordinate a variety of cellular processes and functions in tissue development and regeneration, including wound healing, and ventricular remodeling following myocardial infarction (78). Upregulation of Postn was reported in different diseases characterized by oxidative stress and inflammatory response (79) (80). In addition, upregulation of Postn suppresses SLC7A11 expression through the inhibition of p53 in VSMCs and increased the sensitivity of cells to ferroptosis. More importantly, Postn could remodel the injury environment and was identified as a therapeutic agent for traumatic injury of the CNS. Astroglial-derived Postn promotes axonal regeneration through FAK and Akt signaling pathway after spinal cord injury (81). But Postn was also reported to be a key player in scar formation after traumatic SCI. Genetic deletion of Postn in mice reduced scar formation at the lesion site by inhibiting the proliferation of the pericytes. Moreover, the pharmacological blockade of Postn restrained scar formation and improved the long-term functional outcome after SCI (82). Conversely, our ST maps found the expression of Postn predominantly in the DAF of cluster 1 and cluster 2. Although the expression of Postn also peaked at 7 dpi, it maintained a high level of expression until 28 dpi (figure 4L). We also found the non-canonical genes such as Kininogen 2 (Kng2) (figure 4M) and Serum Amyloid A3 (Saa3) displayed better specificity than canonical fibroblast markers such as P4HB which was expressed in multiple cell types (figure S3A). The expression of Saa3 in cluster 1 increase gradually and reached the highest level at 7 dpi (figure 4N). Saa3 is a pseudogene, and acts as an inflammatory ligand and target of long non-coding RNA (lncRNA) metastasis associated lung adenocarcinoma transcript 1 (MALAT1). Saa3 was accompanied by increase of inflammatory mediators, tumour necrosis factor alpha (TNF-α) and interleukin 6 (IL-6) (83). Additionally, Saa3 is a key mediator of the protumorigenic properties of cancer-associated fibroblasts in pancreatic tumors (84). So far, the role of Saa3 in scar formation is unclear.
In spinal cord, fibrotic scar become mature gradually (figure S3C) and limits CNS regeneration in adult mammals. Similarly, skin wounds generally heal by scarring, a fibrotic process mediated by the Engrailed-1 (En1) fibroblast lineage (85). Preventing Engrailed-1 activation in fibroblasts yields wound regeneration without scarring. It is unknown whether fibrotic scar-forming in spinal cord need Engrailed-1 activation. In the current study, ST maps showed En1 almost was rarely expressed in spinal cord before and after injury (figure S3D), suggesting CNS having a different fibrotic scar-forming. Another study found Jun promotes hypertrophic skin scarring via CD36 in vitro and in vivo models (86). Indeed, we observed that injury enhance Jun expression (figure S3D). Particularly, CD36 was obviously induced in the injury and very limited to the area of scar (figure 4O). Previous study has found CD36 deletion improved recovery from spinal cord injury in mouse, but they didn’t identify the cell type expressing CD36 (87). In our data, there's a fair amount of overlap between CD36+ spots and a-SMA+ (gene name Acta2) activated fibroblasts (figure S3A), suggesting CD36+ involved in contractures of the activated fibroblasts for scar formation. Specially, 59.83% percent of CD36+ cells were GLAST positive at 28 dpi. Targeting CD36 may be a direction for spinal scar therapy. Taken together, we identified the spatiotemporal distribution and origin of fibroblast subtypes at the injury site after SCI.
ST analysis of microglia heterogeneity reveals the specific roles during gliosis
Microglia play a crucial part in scar formation (70, 88). To assess the cellular heterogeneity among microglia at the injury site, re-clustering analysis was performed and visualized on spatiotemporal maps (figure 5A) and a UMAP (figure S4A-B). We identified six microglia subtypes, which were labeled cluster 0 to 5 (figure 5B, figure S4C). Base on the spatiotemporal maps, we found cluster 1 cells mainly distributed in the boundary of fibroblast foci in scar at 7dpi and diminished rapidly. On the contrary, cluster 0 cells replaced cluster 1 cells in the boundary of fibroblast foci. The fibroblast number in Cluster 0 increased gradually from 7dpi to 28dpi (figure 5B). ST maps characterized the expression of P2Y12-a marker for homeostatic microglia, CD68 and Iba1-the markers for activated microglia (figure S4D).
To characterize the functions of these clusters, we listed top 5 differential genes for each cluster (figure 5C), and enriched the GO terms base on highly expressed genes (figure S4E). The top terms showed cluster 0 cells have high enrichment factor in axon ensheathment, axon development and ensheathment of neurons. The pathway of ECM-receptor interaction was also significantly enriched in cluster 0 cells. Interesting, ST maps showed the minority of cluster 0 cells crossed the boundary of fibroblast scar. Cluster 1 cells, cluster 2 cells, cluster 3 cells exhibited significant enrichments for angiogenesis, extracellular matrix organization and extracellular structure organization. The GSVA score of cardiac muscle contraction pathway analysis peaked at 7 dpi and might explained the phenomenon of scar contracture (figure 5J). In addition, cluster 5 cells exhibited significant enrichments for negative regulation of cell motility, negative regulation of cell migration, negative regulation of locomotion and r negative regulation of immune system process, which were similar to the function of neonatal microglia (88). Unfortunately, cluster 5 cells are few, and was not enough to change the whole inflammatory environment.
We next investigated the communication and interaction between microglia and fibroblast, the interface regions of clusters were selected with the range of at least 2 spots for each cluster. The top enriched gene-pairs included TNFRSF1B_GRN, TNFRSF1A_GRN, SPP1_CD44, PLXNB2_PTN, GPR37L1_PSAP and GPR37_PSAP (figures 5D and 5E, figures S4F and S4G). For example, PLXNB2 are transmembrane receptors that participate in axon guidance and cell migration in response to semaphorins (89). Plexin B2 was induced in injury-activated microglia and macrophages early after SCI and facilitated the formation of concentric rings of microglia and astrocytes around the necrotic core of the lesion in a mouse model of SCI (90). ST maps showed Plexnb2+ cell enriched in scar area (figure 5K). Additionally, GPR37L1_PSAP and GPR37_PSAP were top 2 enriched gene-pairs between microglia cluster 0/1/3 and fibroblast cluster 1 from 14 dpi to 28 dpi. Prosaposin (gene name PSAP) can be secreted from various cell types in response to cellular stress. Secreted PSAP can initiate endocytose or pro-survival signaling pathways via binding to GPR37 and GPR37L1 (91, 92). ST maps showed PSAP was enriched in scar-resident fibroblasts (data not shown) and GPR37L1 was enriched in scar-resident microglia (figure 5L), but the exact role of fibroblast secreted PSAP s should be explored in future.
Subsequently, we conducted the trajectory analysis (figure 5F-I) and the pseudo-time trajectory algorithm (figure 5G). The trajectory direction started from cluster 1 cells and cluster 0 cells to a separate branch consisting of cluster 2-5 (figure 5H). We observed that the expression level of enriched genes associated with axon ensheathment and axon development, such as Tnc, Fn1, Sirt2, Ugt8a and Klk8, gradually increased following the scar stabilization. On the contrary, Sema3b, Mmp12 and Laminin gradually declined. We also observed that the expression level of G protein-coupled receptor 17 (Gpr17) significantly increased. Gpr17, a P2Y-like receptor, may act as a 'sensor' of damage that is known to be activated by both uracil nucleotides and cysteinyl-leukotriene released in the lesioned area, and could also participate in post-injury responses (93), In non-injured spinal cord parenchyma, Gpr17 is present on a subset of neurons and of oligodendrocytes at different stages of maturation, whereas it is not expressed by astrocytes. Induction of SCI resulted in marked cell death of Gpr17+ neurons and oligodendrocytes inside the lesion followed by the appearance of proliferating Gpr17+ microglia/macrophages migrating to and infiltrating into the lesioned area (94). In vivo pharmacological or biotechnological knock down of Gpr17 markedly prevents brain infarct evolution, suggesting Gpr17 as a mediator of neuronal death at this early ischemic stage (95). Gpr17 also acts as an intrinsic timer of oligodendrocyte differentiation and myelination. Pharmalogical targeting of Gpr17 signaling in OPCs and microglial inhibition of oligodendrocyte maturation together promote robust myelination of regenerated axons after CNS injury (96). The classically activated neuroinflammatory microglia induce A1 astrocytes by secreting Il-1α, TNF and C1q, and that these cytokines together are necessary and sufficient to induce A1 astrocytes (97). ST maps characterized the enhanced expression of C1qa in six microglia subtypes at different time point (figure S4H).
Recent study found microglia-organized scar-free spinal cord repair in neonatal mice. Fibronectin+ microglia in cluster 3 (MG3) formed bridges between the two stumps starting from 3 dpi. Activated microglia were observed inside the lesion in the absence of GFAP+ astrocytes and collagen I+ fibroblasts (88). In our adult mouse model, ST maps showed A few microglia in cluster 0 crossed the fibroblasts scar and formed bridges between the two stumps at 14 dpi. Unfortunately, the bridges eventually disappeared at 28 dpi. We profiled the selected cells that crossed the fibroblasts scar and found Nxpe3, Pirb and Nr4a1 specificity highly expressed in these cells (figure 5M). IHC verified the expression of Nxpe3 in the boundary of fibroblast scar (figure 5N). Finally, we profiled ECM-related genes (for example, Fn1, ECM1 and Thbs1), wound healing genes (for example, Anxa1), positive regulation of cell adhesion (for example, Spp1), negative regulation of hydrolase activity (for example, Cstb, Serpinb6a and Stfa1), negative regulation of immune system process (for example, Arg1), disease-associated genes (for example, Igf1 and Clec7a), embryonic stem cell-related genes (for example, Ms4a7, Lgals and Ms4a6c), phospholipase inhibitor activity (for example, Anxa2, Anxa5 and Apoc1) in six microglia subtypes. These genes were highly expressed in MG3. But, we found that Stfa1 and Clec7a weren’t expressed in six microglia subtypes, the expression level of Arg1 and Thbs1 were very low (figure 5O). Specially, Anxa1, a potent anti-inflammatory protein, may contribute to the ability of MG3 cells to mediate rapid inflammation resolution and were barely expressed in the microglia of cluster 2 and cluster 3, which were the two major types of microglia in the scar at 28 dpi (figure S4I anf figure 5B). These difference between the microglia in six clusters with the MG3 cells in neonatal mice may contribute to decipher the difference in regenerative ability between adult and neonatal mice.
Six different astrocyte subtypes are identified by ST.
Astrocytes are thought to be the main drivers of encirclement. Tightly connected astrocytes continuously reshape the boundary at the edge of injury, wrapping immune cells and fibroblast-like cells with ephrin-mediated cell adhesion and spatially isolating remaining neural tissue from injury and fibrosis (59, 74, 98). We next assessed the cellular heterogeneity among astrocyte at the scar, re-clustering analysis was performed and visualized on spatiotemporal maps (figure 6A) and a UMAP (figure S5A-B). We identified six astrocyte subtypes, which were labeled cluster 0 to 5 (figure 6B, figure S5C). Base on the spatiotemporal maps, we found cluster 2 cells mainly distributed in the grey matter layer of microglia concentric rings in scar at 7dpi and diminished rapidly. Cluster 2 cells might represent protoplasmic astrocytes. On the contrary, cluster 0 cells mainly distributed in the white matter layer of scar, suggesting these cells represented fibrous astrocyte. The astrocyte number in Cluster 0 peaked at 14 dpi (figure 6B). The astrocyte number in cluster 1 continuously increase. ST maps characterized the expression of GFAP and Lcn2-the markers for activated astrocyte (figure S5D-E).
To further explain the biological process of the subgroup, we characterized top 5 differential genes for each cluster (figure 6C), and enriched the GO terms base on highly expressed genes (figure 6D). The top terms showed cluster 0 cells have high enrichment factor in traditional functions (99), such as neurotransmitter transport and synaptic vesicle cycle. In addition, cluster 0 cells also exhibited significant enrichments for myelination, axonogenesis, axon ensheathment, axon development and regeneration. The GSVA score of synaptic vesicle cycle pathway analysis and axon guidance also supported this view (figure 6E). Cluster 1 cells also exhibited significant enrichment for ATP metabolic process. The corresponding KEGG analysis showed cluster 1 cells have high enrichment factor in cardiac muscle contraction pathway. The GSVA score showed it peaked at 7dpi around the microglia scar (figure 6E). The top terms showed cluster 2 cells have high enrichment factor in regulate brain energy metabolism and homeostasis (59), such as ATP metabolic process and generation of precursor metabolites and energy. Interesting, cluster 4 cells exhibited significant enrichments for osteoblast differentiation, regulation of ossification, cellular response to zinc ion and cellular response to copper ion. Astrocytes were known to remove excessive potassium ions in the microenvironment (100). Here, we found cluster 4 cells contributed to the homeostasis in zinc ion /copper ion, which caused neuronal cell death (101).
We next conducted the trajectory analysis (figure 6F) and the pseudo-time trajectory algorithm (figure 6G). The trajectory direction started from cluster 2 cells to a separate branch consisting of cluster 0, 1, 3, 4, 5. We observed the genes with significantly elevated levels of expression were associated with immune response Lectin induced complement pathway, such as complement C1q Chain A (C1qa), complement C1q Chain B (C1qb) and complement C1q Chain C (C1qc) (figure 6H). Slute carrier family 25 member 22 (Slc25a22), solute carrier family 7 member 10 (Slc7a10) and solute carrier family 1 member 2 (Slc1a2) were involved in homeostasis of excitatory neurotransmitter glutamate, and markedly decreased. At the same time, solute carrier family 32 member 1 (Slc32a1), involved in gamma-aminobutyric acid (GABA) and glycine uptake into synaptic vesicles, also markedly declined. The relation between the local excitatory disorder after SCI and these transporters downregulation remain unknow (41). In addition, the pro-regeneration molecular inhibin subunit beta A (Inhba) markedly declined followed by the increase of Semaphorin 3B (Sema3b) (figure 6I), which inhibits axonal extension. These changes suggested the microenvironment of the astrocyte scar switched to inhibit regeneration with the maturation of scar (figure S5F) (60, 102).
Subsequently, we next investigated the communication between astrocyte and microglia. The top enriched gene-pairs included TNFRSF1B_GRN, TNFRSF1A_GRN, SPP1_PTGER4, SPP1_CD44, SPP1_a9b1 complex, SPP1_a4b1 complex, PLXNB2_SEMA4D, PLXNB2_PTN, GPR37L1_PSAP, GPR37_PSAP, CD74_MIF and CD74_APP (figure 6J). Interesting, GPR37L1_PSAP, GPR37_PSAP and PLXNB2_SEMA4D were in the top relationship between microglia and fibroblast. Here, these pairs were still critical. For example, GPR37L1_PSAP and GPR37_PSAP were top 2 enriched gene-pairs between astrocyte cluster 2/4 and microglia cluster 0/1.
Additionally, the expression level of Sema4d gradually increased, and the pair of PLXNB2_SEMA4D between astrocyte cluster 2/4 and microglia cluster 1 has been enriched. These data indicated there were some common mechanisms mediate cell connections in the scar.
Finally, we profiled the single-cell regulatory network inference and clustering (SCENIC) analysis(103) and found the six subgroups were grouped into two major groups (figure 6J). The cell number of groups 1 (cluster 0, cluster 1 and cluster 4) increased gradually. Groups 1 cells might represent the traditional A1 astrocyte (97). In addition, ST maps showed that group 1 cell wrapped microglial scar at intermediate stage (14 dpi, 28 dpi). Interesting, we found 10 transcription factors drived astrocyte transformation to the A1 like-groups 1 cells. Additionally, Eph/ephrin signaling can the activated astrocytes maintain the homeostasis of extracellular glutamate. The EphA4 signaling prevents glutamate excitotoxicity (104). ST maps showed that the expression level of EphA4 was quite low, suggesting the neurotransmitter disorder in scar (figure S5G). Transglutaminase 2 (TG2) plays a key role in regulating the response of astrocytes to damage. Astrocytic-specific TG2 deletion or inhibition results in a significant improvement in functional recovery after SCI (105). ST maps showed TG2 was upregulated in astrocytes in scar at different time point after SCI (figure S5H). Glial BAI3 and BAI1 binding to RTN4/NoGo receptors regulate axonal elongation and synapse formation, thereby controlling neural network activity (106). Our ST showed the expression level of Bai3 in astrocyte scar was low at different time point after SCI (figure S5I). Additionally, astrocyte-derived saturated lipids contained in APOE and APOJ lipoparticles mediate toxicity. Eliminating the formation of long-chain saturated lipids by astrocyte-specific knockout of the saturated lipid synthesis enzyme ELOVL1 mitigates astrocyte-mediated toxicity (107). Our ST showed the activated astrocytes in scar have high expression of Apoe, Apoj and Elovl1 (figure S5J). An unexpected finding was that the activated astrocytes in scar have high expression of Clu (figure S5K), an anti-neuroinflammatory gene. The previous study showed that CLU was secreted by Schwann cells in peripheral sciatic nerve and induced outgrowth of sensory neurons, but not motor neurons (50). Recently, an exciting result showed 'runner plasma', collected from voluntarily running mice and infused into sedentary mice, reduces baseline neuroinflammatory gene expression and experimentally induced brain inflammation. Mechanically, the complement cascade inhibitor CLU in plasma reduced neuroinflammatory gene expression in a mouse model of acute brain inflammation and a mouse model of Alzheimer's disease (108). The exact role of Clu in spinal cord injury remains unknown. In a word, these findings call for a reinterpretation of astrocyte infiltration into the scar and may inform future therapeutic approaches.
ST identified the different oligodendrocyte subtypes in the scar
Oligodendrocyte wrapped around the axons of neurons in the CNS, forming myelin insulators on the surface of neurons that allow electrical signals to travel more efficiently. Interesting, subpopulations of oligodendrocyte and their progenitor cells have much in common with immune cells in mouse models of multiple sclerosis, and they can be involved in clearing myelin sheaths damaged by disease, in a manner similar to immune cells (109). We next assessed the cellular heterogeneity among astrocyte at the scar, re-clustering analysis was performed and visualized on spatiotemporal maps (figure 7A) and a UMAP (figure S6A-B). We identified three oligodendrocyte subtypes, which were labeled cluster 0 to 2 (figure 7B, figure S6C). Base on the spatiotemporal maps, we found cluster 0 cells mainly distributed in the outer layer of astrocyte concentric rings in scar at 3 dpi and diminished rapidly. On the contrary, cluster 2 cells mainly distributed in the white matter of scar (figure 7A). The oligodendrocyte number in cluster 1 continuously increase. ST maps characterized the expression of Mbp, Mog, Mag, and Cldn11-the markers for activated astrocyte (figure S6D-E).
To further explain the biological process of the subgroup, we characterized top 5 differential genes for each cluster (figure 7C), and enriched the GO terms base on highly expressed genes (figure S6F). The top terms showed cluster 0 cells have high enrichment factor in response to wounding. In addition, cluster 0 cells also exhibited significant enrichments for epithelial cell proliferation, regeneration, angiogenesis, axon regeneration and glial cell migration. The GSVA score of p53 signaling pathway analysis and also supported this view (figure 7D). But cluster 0 cells diminished gradually. Additionally, cluster 1 cells exhibited significant enrichments for neutrophil activation involved in immune response, leukocyte mediated immunity, antigen processing and presentation and macrophage activation. The GSVA score of antigen processing and presentation analysis enriched in cluster 1 cells significantly increased from 7 dpi to 28 dpi (figure 7D), suggesting they can work in a manner similar to immune cells (109). Cluster 2 cells exhibited significant enrichments for neurotransmitter transport, regulation of synaptic plasticity and neurotransmitter secretion. The GSVA score of GABAergic synapse in top 3 pathway analysis enriched in cluster 2 cells significantly increased from 7 dpi to 28 dpi (figure 7D), which may be related to the inhibition of neuronal excitability after injury (41).
We next conducted the trajectory analysis (figure 7E) and the pseudo-time trajectory algorithm (figure 7F). The trajectory direction started from cluster 0 cells to a branch consisting of cluster 1 and 2. We observed the expression level of serpin family E member 1 (Serpine1), a component of innate immunity, gradually declined. But beta-2-microglobulin (B2m), a component of the class I major histocompatibility complex (MHC), was involved in the presentation of peptide antigens to the immune system and gradually increased. Galectin 1 (Lgals1), playing a role in regulating apoptosis, gradually declined. Interesting, solute carrier family 6 member 1 (Slc6a1) and solute carrier family 32 member 1 (Slc32a1) mediated the rapid removal of GABA and maintain low extracellular levels. The expression level of them significantly increased, suggesting oligodendrocyte also participated in the regulation of excitability in the scar. In addition, cell adhesion molecules neurotrimin (Ntm) exerts an inhibitory impact on neurite outgrowth (110), whereas its expression level significantly increased during the maturation of scar.
Finally, we profiled SCENIC analysis (figure 7I) and found Foxc1 potentially controlled cluster 0 cells. Foxc1 was as an important regulator of cell viability and resistance to oxidative stress in the eye (111). Irf7 and Ifr8 controlled cell fate of cluster 1 cells. Interesting, Sox10 controlled cell fate of cluster 2 cells. Sox10 were verified to play a central role in oligodendrocyte maturation and CNS myelination. Specifically, Sox10 activates expression of myelin genes, during oligodendrocyte maturation, such as Dusp15 and Myrf (112).
Oligodendrocyte death after SCI contributes to demyelination of spared axons. Bone morphogenic protein-7 (BMP7) could potently inhibit TNF-α-induced oligodendrocyte apoptosis (113). Overexpression of BMP7 reduced oligodendrocytes loss and promoted functional recovery after spinal cord injury (114). In our data, ST maps showed cells secreted BMP7 were mainly in the fibroblast scar, suggesting its beneficial role (figure 7J). Last, ST maps characterized the expression of CD74-a component of the class II major histocompatibility complex (MHC) (figure 7K). CD74 is a critical chaperone that regulates antigen presentation for immune response (115). All the three oligodendrocyte subtypes have a moderate expression of CD74. Taken together, these findings call for a reinterpretation of oligodendrocyte infiltration into the scar and deciphered the potentially functions.
Gene regulatory networks controlling the scar programming
To better understand scar-relevant changes in genes regulation and interactions between cell types, we carried out WGCNA and an unbiased coexpression analysis (38) of our ST data. In WGCNA, 2000 genes (Table S4) from the spatiotemporal transcriptome dataset in cluster 2 and cluster 3 (figure 1H) are clustered. The branches of highly correlating genes were formed, which were cut and assigned a color. Finally, the correlation of genes and modules was calculated (figure 8A), and 14 modules were identified (figure 8B and figure S7A). Among the 14 modules, the MEgreenyellow module had the most significant correlation with the time (figure S7A), which gradually increasing following the scar formation and approximately representing microglia and astrocyte (figure 8B). We next conducted the network of the enriched hub module genes involved in Megreenyellow to identify key controlling genes (figure 8C). Meblue module and Meblack module had the most significant correlations with the early acute stage (3 dpi) and represnting macrophage and oligodendrocyte, and hub module genes in Meblue (figure 8D) and in Meblack (figure S7B) were respectively carried out.
Subsequently, we adapted another unbiased coexpression analysis more suitable for ST data (38). We identified 18 major coexpression modules on the basis of 2,000 highly variable genes from all spatiotemporal transcriptome dataset (Table S5). Among them, module 3, 8, 10, 16, and 17 are highly correlated and used as candidate analysis modules (figure 8E). ST maps showed the spatiotemporal activities of module 3, 8 and 10 (figure 8F). We next grouped the genes of module 3, 8 and 10 on the basis of their expression pattern, resulting in submodules (figure 8G, figure S7C-7F). Col3a1, col1a2 and col1a1 (submodule 8.1) have highest corelation score (figure 8G). KEEG pathways enrichment analysis of the submodules 8.1 show that they are enriched for ECM-receptor interaction (figure 8H). Gpnmb, Lyz2, Ftf1, Lgals3, Ctsd and Ctsb (submodule 10.3) have highest corelation score (figure S7E). KEEG pathways enrichment analysis of the submodule 10.3 show that they are enriched for lysosome (figure S7F). Collectively, the pathway activities encompassed by module 8 and 10 reveal signaling within and between cell types during glial and fibroblastic activation in the scar.
The scar boundary was reinterpreted by ST
The traditional characteristics of lesion site in mice after adult injury were scare structure visualization with immunofluorescence for laminin, fibronectin, chondroitin sulfate proteoglycan (CSPG), GFAP, CD68, CD31 and collagen (88). All these data came from the sections and we lack a global understanding of its distribution throughout the spinal cord scar. Although diffusion tensor imaging (DTI) could accurately provide information of cavity area after SCI, this method is difficult to distinguish cell types (116, 117). With the help of tissue transparency, laminin, fibronectin and collagen were re-characterised in our mice model of T10 right lateral hemisection (figure 9A). Laminin was highly enriched in the wound area at 7 dpi. On the contrary, fibronectin and collagen was relatively low expression, which were different form the expression pattern of neonatal mice in spinal cord crush. We next applied our ST data to exhibit the laminin (figure 9B), fibronectin (figure 9C). Overall, the ST maps showed fibronectin was upregulated in the injury site and roughly characterised the region of scar at 3 dpi, 7 dpi and 14 dpi. On the contrary, the ST maps showed laminin was upregulated in the injury site and clearly characterised the region of scar at 7 dpi, 14 dpi and 28 dpi. The ST maps of CD31 showed the neovascularization failed to penetrate the scar (figure 9D), which was similar to slices staining (88). Collectively, it is a potential breakthrough that combining the ST maps of scar markers to identify the size of scar. To define the distribution of each cell type more accurately, we counted the number and fraction of each cell type maintaining scar architecture at different time after SCI (figure 9E and 9F), and found that fibroblasts are always in the center during scar shrinkage and maturation. Fibrotic scars are surrounded by microglial scars, which are wrapped by the astrocytic scars (118). Specially, microglial scars and astrocytic scars did shrink at the gray matter. But they expanded a larger area along white matter at 28 dpi (figure 9E and 9F). Conventional wisdom holds that secondary injury can be temporally divided into acute, sub-acute, and chronic phases following SCI (1, 57). In our study, ST decipher all or most genes systematically throughout scar formation. Combining cell type, functional diversity, the possible trajectory and the fraction of scar-resident cell, we defined four phases of the scar formation: macrophage infiltration, proliferation and differentiation of scar-resident cell, scar emergence and scar stationary. Taken together, we identified the spatiotemporal dynamics of key genes in specific cell types, thus providing a global atlas of gene expression and molecular regulation driving scar formation following SCI (figure 9G).