Graphical abstract
graphic file with name fx1.jpg
[37]Open in a new tab
Highlights
* •
Transcriptional landscape of motoneurons (MNs) and its niches upon
axonal injury
* •
High dynamic niches and cell heterogeneity around injured and the
regenerating MNs
* •
Screen critical master regulators that may initiate regeneration
program in MNs
* •
Knockdown of Arid5a significantly reduces neurite outgrowth
in vitro
__________________________________________________________________
Molecular physiology; Neuroscience; Omics
Introduction
Peripheral nerve injury (PNI) is a common and severe clinical problem
([38]Li et al., 2020b). Unlike limited regenerative capacity after CNS
injury, axons of sensory and spinal motor neurons could regenerate and
rebuild neuromuscular connections after PNI. Extensive studies have
explored molecular changes and cellular heterogeneity in sensory
neurons after PNI ([39]Avraham et al., 2020; [40]Renthal et al., 2020;
[41]Wang et al., 2021). However, limited studies have investigated
those in spinal motoneurons (MNs) after PNI. Rodent sciatic nerve crush
(SNC) injury is a well-characterized and robust model for studying
peripheral nerve regeneration. In addition to the activation of the
intrinsic regrowth program, the microenvironment also plays important
roles in the successful axon regeneration, including resident glia,
peripheral infiltrated immune cells (e.g. macrophage), vascular cells,
and extracellular matrix homeostasis ([42]Lei et al., 2012; [43]Scheib
and Hoke, 2013). A pioneering study has explored molecular change under
translating in MNs using RiboTag; ChAT^Cre mouse model at 7 days
post-SNC ([44]Shadrach et al., 2021). However, this only covered
limited mRNA forming their protein products and dynamic gene expression
at the genome-wide level in MNs in the first seven days when undergoing
dramatic change in gliosis and angiogenesis ([45]Milich et al., 2021;
[46]Yu et al., 2019) after injury is still unknown.
Laser capture microscopy (LCM) combined with antibody rapid staining
could capture cell types of interest in situ and coupled with
Smart-seq2 (LCM-seq) could provide the accurate quantification of gene
expression in dissected cells ([47]Nichterwitz et al., 2016). LCM-seq
has been successfully applied in studies of nervous system diseases,
including Parkinson disease and amyotrophic lateral sclerosis (ALS)
([48]Aguila et al., 2021; [49]Nichterwitz et al., 2016, [50]2020).
Unlike other cell isolation methods, such as fluorescence-activated
cell sorting and microfluidics, LCM could harvest target cells from the
microenvironment without dissociation of the tissue and preserved in
situ spatial information ([51]Nichterwitz et al., 2016). Though
inevitable encompassing close its neighboring cells or non-cells around
a target cell, LCM-seq just also provides additional information for
its microecology.
MNs are large and sparse neurons in the ventral horn of the spinal cord
which is difficult and inefficient to obtain sufficient and qualified
MNs. Moreover, dynamic cellular heterogeneity around MNs during injury
and regeneration that nuclei isolation would lose substantial
information. Thus, we employed LCM-seq to investigate the dynamic
expression profile of spinal MNs soma and its minimum niches during
axonal injury and regeneration. We investigated six timepoints post-SNC
and the naive rat with approximately a single MN resolution (a motor
neuron was contained in a sample). Comparison with another study
([52]Shadrach et al., 2021) in gene expression of injured MNs and our
immunostaining indicated the high quality of our data in
reproducibility and accuracy. In the SNC rodent model, not all dorsal
root ganglion (DRG) neurons or motoneurons are axotomized (∼50%)
([53]Holland et al., 2019; [54]Shadrach et al., 2021). Based on
clustering analysis and expression of neuronal injury marker (Atf3), we
could well distinguish “intact” and “injured” MNs at 1, 3, and 7 days
and estimated axotomized MNs ratio were approximated to previous
studies ([55]Holland et al., 2019; [56]Shadrach et al., 2021). We found
numerous robustly differentially expressed genes (DEGs) reported as
regeneration-associated genes (RAGs) in DRG after PNI was robustly
upregulated in “injured and regenerating” MNs by comparison with the
naive and “intact” MNs. We also demonstrated highly dynamic
microenvironments around “injured and regenerating” MNs with dramatic
changes of genes related to the response of MNs, microglia/macrophage,
vascular cells, and extracellular matrix reorganization. We also
demonstrated critical regulators and a potential injury and
regenerative-related co-expression module post-SNC. Specifically, we
found a candidate master regulator (a transcription factor and
RNA-binding protein, ARID5A) which may be associated with axon
regeneration, presenting robust upregulation in injured neurons
(including sensory and motor neurons) and co-expression with other
well-known regulators (Atf3, Klf6, Myc, and Sox11) and significantly
inhibiting neurites outgrowth of primary sensory neurons by knockdown
of Arid5a in vitro. Finally, we present a web-based resource for easily
exploring gene expression of MNs after PNI (MNniches-seq). Our dataset
not only provides the first temporal-spatial profile of individual MNs
and its niche providing dynamic microenvironment information during
axonal injury and regeneration, which are crucial for studying nerve
regeneration and motoneuron degeneration disease, such as ALS.
Results
Two distinct clusters were revealed by LCM-seq
We employed a sciatic nerve crush (SNC) rat model, and investigated
timepoints: 0.5, 3, 12 h, 1 d, 3 d, and 7 d (hours/days), and compared
to that from naive rats to study molecular dynamics for axonal injury
and regeneration ([57]Figure 1A). Motoneurons located in the ventral
anterior horn of the lumbar spinal cord (L4/5) were labeled using rapid
immunofluorescent staining of ChAT antibody and a circular region of
the diameter around 1.5–2 times of a large motoneuron soma was cut
using LCM under direct microscopic visualization, subjected to
Smart-Seq2, and sequenced ([58]Figures 1A and [59]S1). Finally, a total
of 700 samples were sequenced successfully with a depth of 44–134
million clean reads per sample. Samples with low expression of Chat
(FPKM<1) were discarded and left 682 samples (94–104 per condition) and
a total of 16,163 protein-coding genes were expressed, at least three
samples in one group (FPKM ≥1, 12,098–13,074 genes per condition,
[60]Figure 1C). Pearson’s correlation analysis showed a high
correlation among dissected samples in each condition ([61]Figure 1B).
Neuronal activity was estimated with the AUCell package using neuronal
markers ([62]Aibar et al., 2017; [63]Cahoy et al., 2008) and showed a
high neuronal signature with an AUC score of 676 samples (∼99%) greater
than the estimation (0.12, [64]Figure 1D). Besides, high expression of
markers for alpha motoneurons (Rbfox3 and Spp1) and low expression of
the maker for gamma motoneurons (Htr1a) demonstrated that dissected
samples here almost encompassed alpha motoneurons ([65]Figure S2).
Dimension reductions analyses demonstrated a clear two separated
clusters with one consisting of naive samples and early injured samples
(0.5–12 h) and most samples from 1 day (93.9%) and samples
(47.9%–62.9%) from 3 to 7 days, while the other consisting of 37.1% and
52.1% of samples from 3 to 7 days, respectively ([66]Figures 1E and
1F). Identification of the top 20 markers showed that genes with high
expression in cluster 2 were related with the response of neuron and
glia to axonal injury, such as neuronal released factors (Atf3, Ccl2
([67]Shadrach et al., 2021), Reg3b ([68]Namikawa et al., 2005), and
Serpina3n (LOC299282) ([69]Vicuna et al., 2015)) microglial activation
(Reg3b ([70]Li et al., 2020a; [71]Namikawa et al., 2005), C1qa, Csf1r,
Anxa3 ([72]Zhang et al., 2021), and Lyz2), and astroglial activation
(Gfap, Il33 ([73]Gadani et al., 2015b), and Timp1 ([74]Ogier et al.,
2005)). Anxa3 (annexin A3) was most expressed in microglia/macrophages
and endothelial cells in the spinal cord ([75]Blum et al., 2021) and
functions in inflammation, angiogenesis, and rat liver regeneration
([76]Harashima et al., 2008; [77]Park et al., 2005). Immunostaining of
ANXA3 showed obviously increased ANXA3^+ cells in ventral horn were
detected as early as 1 day, but close contact between ANXA3^+ cells and
motoneurons were obviously increased at 3 and 7 days reflected by
adhesion to the neuronal perikaryon, which is consistent with the
result of IBA1 immunostaining ([78]Figures 1I and [79]S3), suggesting
increased crosstalk in neuron-microglia. This was generally consistent
with our sequenced data. The top 20 upregulated genes in cluster 2 also
included known RAGs (Atf3, Gap43, Ngfr, Timp1, and Ccl2), energy
metabolism (Ldha), ECM (Timp1) and Timp1 receptor (Cd63, a
tetraspanin), and pro/anti-inflammation (Reg3b, Ccl2; Il33). Timp1
interacting with Cd63-activating β1 integrin signaling could activate
cellular signaling cascades, promoting cell survival and migration and
regulating cytoskeletal reorganization ([80]Lee et al., 2014). Taken
together, we suspected that samples in cluster 2 contained injured MNs
and underwent regenerating as discussed below.
Figure 1.
[81]Figure 1
[82]Open in a new tab
Two distinct clusters in ventral lumbar spinal motoneurons post-SNC
(A) Overview of the investigated timepoints and workflow of the
experiment.
(B) Pearson’s correlations among samples captured based on expressed
genes (gene with average FPKM >1 was kept).
(C) Gene number identified per condition. The number above the boxplot
represented the sample number.
(D) Histogram of AUC score of neuronal markers reflecting the high
activity of neuronal characteristics.
(E) Principal component analysis (PCA) visualization of samples
distribution.
(F) UMAP clustering of samples. Two distinct clusters were shown.
(G) The sample proportion of each group from different injury
timepoints in two clusters.
(H) Expression profiles of top 20 differential expression genes in
either cluster 1 or cluster 2. (I) Immunostaining of ANXA3 (red) in rat
ventral lumbar spinal cord at different injured timepoints post-SNC.
Neuron and nuclear were labeled with NeuN (green) and DAPI (blue).
Scale bar: 100 μm.
Three response stages were revealed by LCM-seq
Atf3 is a sensitive and specific neuronal marker with injured
peripheral axons ([83]Tsujino et al., 2000). Separate UMAP plots of
Atf3 expression showed high expression in samples from cluster 2 and
also showed high expression in part samples from cluster 1 at 1 day
post-SNC ([84]Figure 2A). We found the expression of Atf3 was greater
than 1 in almost samples from cluster 2 while less than 1 in other
samples except for part samples from cluster 1 at 1 day (26 samples,
26.5%, [85]Figure 2A). Interestingly, this ratio (26.5%) was
approximate to the estimation (∼23.8%) of injured sensory neurons at
1 day post sciatic nerve injury ([86]Renthal et al., 2020). Only 1–2
samples in cluster 2 from 3 to 7 days showed expression of Atf3 was
less than 1 that we directly assigned based on the clustering results
([87]Figure S4). Namely, 37.1% and 52.1% of “injured” MNs were
estimated at 3 and 7 days respectively, which were also consistent with
the estimation in motoneurons using ATF3-CreERT2 knockin mouse post-SNC
(45
[MATH: ± :MATH]
2% for 4 days) ([88]Holland et al., 2019). Thus, we divided these
samples into three response stages and termed these stages with stages
1–3 (namely 0.5–12 h, 1 day, and 3∼7 days) and three states (states
1–3). Atf3 was significantly upregulated with remarkable fold-change in
state 2 and state 3 (81- to 82-fold at 1 day, 122-fold at 3 days, and
131-fold at 7 days) while no significance in state 1 ([89]Figure S5).
Our estimation is approximate to the fold-change (134-fold) of Rpl22^HA
immunoprecipitation (IP) in ipsilateral motoneurons compared to the
contralateral ones from RiboTag; ChAT^Cre mouse model at 7 days
post-SNC ([90]Shadrach et al., 2021), suggesting high reproducibility
and accuracy of our data. Besides, we employed bootstrapping strategy
to identify robust differential expression genes (DEGs) ([91]Aguila
et al., 2021) and showed remarkably robust (≥90%) differential
expression of Atf3 in state 2 and state 3 ([92]Figure S5). The SCENIC
([93]Aibar et al., 2017) analysis showed that the activity of Atf3
regulon was also significantly higher in state 2 and state 3,
indicating Atf3 was mostly activated as early as 1 day (Student’s t
test, ∗p < 0.05, ∗∗p < 0.0001, [94]Figure S5). Immunolabeling of ATF3
in rat spinal cord upon axonal injury supported our sequencing results
that ATF3 was expressed as early as 1 day and co-localized with the
nucleus of part motoneurons, indicating injured MNs ([95]Figure 2B).
ATF3-positive injured MNs most located in ventral posterolateral
nucleus were motor pools for innervation muscles of a limb.
Figure 2.
[96]Figure 2
[97]Open in a new tab
Three response stages in axonal injury and regeneration post-SNC
(A) Separate UMAP plots of Atf3 expression in samples collected from
different post-injury timepoints. The ratio of Atf3 expression ≥1 of
cluster 1 and 2 at 1 day or ratio of cluster 2 at 3 and 7 days post-SNC
was shown.
(B) Immunostaining of ATF3 in the lumbar spinal cord from different
post-injury timepoints.
(C) The radar map depicts the gene number of DEGs (P-value≤0.05 and
|log[2]^fold-change|≥log[2]^1.5) at different injured time points in
different states compared to the naive.
(D) Significantly and consistently upregulated genes at 1 day post-SNC
in samples from state 2 and state 3. DEGs with top rank or critical
were labeled.
(E) Functional enrichment analysis of consistently upregulated DEGs in
1 day post-SNC at state 2 and state 3.
(F) Gap43 expression at different states in different injured time
points post-SNC.
(G) Frequency of detected differentially expressed Gap43 with 100
bootstrap DESeq2 with varied sample sizes.
(H) Expression of well-established regeneration-associated genes
collected from the literature. ∗ indicated transcription factors.
(I) Geneset enrichment analysis of well-known pathways related to
regeneration. Terms with false discovery rates (FDR, FDR<0.25) were
highlighted with ∗ (NOM p < 0.05) and ∗∗ (NOM p < 0.01).
Next, we analyzed DEGs by comparison with the naive to assess whether
separation based on a single marker could well obtain a distinct
expression profile that reflects cell states. The result showed a
significant increase in DEG number after separation with approximately
2-fold of DEG numbers without separation ([98]Figure S6). A smaller
number of DEGs were detected in in state 1 from timepoints (0.5
h∼1 day, <500 DEGs) than that detected from 3 (874 DEGs) and 7 days
(2,049 DEGs) ([99]Figure 2C). A larger DEG number was detected in state
2 and/or state 3 than that detected in the corresponding state 1. A
similar DEG number in state 2 and state 3 was detected at 1 day, while
a larger DEG number was detected in state 2 at 3 and 7 days, suggesting
distinct expression profiles between 1 and 3 days/7 days. 189 and 104
genes showed consistent upregulation and downregulation in state 2 and
state 3 at 1 day, respectively ([100]Figure 2D), and upregulated genes
were significantly enriched in terms of “response to wounding”, “cell
proliferation”, “angiogenesis”, and several signal pathways, such as
MAPK and Jak-STAT signaling pathways ([101]Figure 2E). These evidences
suggested that our sample separation could well represent distinct cell
states by supporting with distinct expression profiles.
Gap43 is a proven marker labeling for developing and regenerating
neurons. Our analysis of Gap43 expression showed significant
upregulation only in state 2 (1.86-fold at 1 day, 4.94-fold at 3 days,
and 8.94-fold at 7 days, [102]Figure 2F). The upregulated fold-change
of Gap43 is consistent with the estimation in MNs (1.96-fold at 1 day,
5.82-fold at 3 days, and 9.93-fold at 7 days) ([103]Shadrach et al.,
2021), confirming our data accuracy and suggesting the regenerating
axon in samples from the cluster 2 population. Gap43 was robustly
upregulated in state 2 neither state 1 nor state 3 (frequency ≥50%,
[104]Figure 2G). Numerous RAGs (Atf3, Gap43, Sprr1a, Jun, Sox11, etc)
with literature evidence showed significant upregulation in state 2 at
3 and 7 days, and some in state 2 and/or state 3 at 1 day
([105]Figure 2H). Notably, important pro-growth transcription factors
(Atf3, Sox11, Stat3, Smad1, Jun, and Myc) showed significant
upregulation as early as 1 day in state 3 and maintain upregulation in
state 2, while terminal RAGs, such as Gap43, Hspb1, Sprr1a, and Ccl2,
showed only significant upregulation in state 2. This result indicated
that state 3 from 1 day was the pre-regeneration phase that activated
pro-regenerative TFs for the initiation of the regenerative program.
Furthermore, most RAGs showed high expression in state 2 from 1 to
3 days, while part RAGs showed high expression in state 2 from 1 day,
suggesting a gradual initiation of regenerating program and 1 day was a
key point for successful regeneration.
In addition, we also investigated the enrichment of several important
pathways related to regeneration, such as “Neurotrophic signaling
pathway”, “MAPK signaling pathway”, “JAK-STAT signaling pathway”, and
reflecting axon regeneration “Regulation of actin cytoskeleton”. The
result showed pro-regenerative pathways were significantly enriched in
state 2 and state 3, but not in state 1 compared to the naive
([106]Figure 2I). Actin cytoskeleton-related term was significantly
enriched in state 2. This result suggested that state 3 from 1 day is a
pre-regeneration phase while state 2 was regenerating phase.
Collectively, dissected samples could be divided into two clusters and
three response phases which could be categorized as “early stress
response phase” (stage 1, 0.5–12 h), “pre-regeneration phase” (1 day),
and “regeneration phase” (3∼7 days) in our investigated injured
timepoints.
Robustly differentially expressed genes of MN niches in response to axonal
injury
We next performed differential expression analysis by comparison with
the naive. We compared our results ( state 2 from 7 days) with DEGs
(301 out of 342 with ortholog in rat were considered) identified by
comparing expression of Rpl22^HA IP in injured and regenerating MNs
with the uninjured in mouse spinal cord at 7 days post-SNC
([107]Shadrach et al., 2021). Significant correlation between these two
datasets (Pearson’s correlation: 0.58, p-value < 2.2e-16) was observed
based on the fold-change of common DEGs compared to the corresponding
control. 189 out of 196 (96.2%) overlapped DEGs showed consistent
either upregulation or downregulation, including Atf3, Bst2, and Cyp1b1
([108]Figure 3A and [109]Table S1). Of the common DEGs, we found
dramatical upregulation of genes was released or expressed by injured
neurons, such as Atf3, Ccl7, Ccl2, Ngfr, Vip, Crh, Gal, Sprr1a, Ankrd1,
Csf1, Adcyap1, and Serpina3n. Besides, we also found genes related to
other cellular (endothelial cells, microglia/macrophage) or
non-cellular (ECM) were present, suggesting their important roles in
axon regeneration. Our result showed a small number of DEGs (<500) at
early injury (0.5–12 h) and state 1 from 1 day ([110]Figure 2C) under
the relaxing threshold. Though a relatively large number of DEGs
(847–2,049) were detected in state 1 at 3 and 7 days ([111]Figure 2C).
However, GO enrichment showed that DEGs in state 1 at 3 and 7 days were
almost enriched in terms related to inflammation and immune response,
such as “positive regulation of cytokine production” and “regulation of
inflammatory response” ([112]Figure S7 and [113]Table S2). Axonal
injury triggers sterile inflammation with peripheral immune cell
infiltration, migration, proliferation and differentiation, and CNS
microglia activation, proliferation ([114]Gadani et al., 2015a),
forming an inflammatory milieu. Thus, DEGs detected in state 1 from 3
to 7 days may be explained by the physiological response to detrimental
stimuli in “intact” motoneurons. Moreover, captured and discriminated
“intact” motoneurons and its minimal niches at 1, 3, and 7 days provide
a valuable control for further screening crucial molecules involved in
axon regeneration.
Figure 3.
[115]Figure 3
[116]Open in a new tab
Robust DE-Gs in response to axon injury and regeneration
(A) Significant correlation of fold-change of DEGs between the study
and our data at 7 days post-SNC.
(B) Venn diagram depicts the overlapped gene numbers by comparison with
the naive or corresponding state 1 following SNC.
(C) Nine robust upregulated DEGs in our study.
(D) Expression of nine robust upregulated DEGs in spinal cord injury
dataset.
(E) GO enrichment of upregulated genes at 1, 3, and 7 days compared to
the naive or corresponding state 1. DEG number annotated a term was
shown within the dot.
(F) The relative fraction of macrophage subtypes estimated in our
collected samples based on the CIBERSORT method.
(G) Top-ranked 100 robustly upregulated DEGs by considering both
expression level and fold-change.
(H) Expression profiles of key genes related to important processes
during axon injury and regeneration.
We next employed bootstrap strategy ([117]Aguila et al., 2021) to
identify robust and reliable DEGs (average frequency ≥50%). 80 robust
DEGs presented in early injury (61 up- and 19 downregulated genes,
0.5–12 h, [118]Figure 3B and [119]Table S3). We found genes, Zfp36,
Egr1, Klf4, Junb, Apold1, Dusp1, and Nr4a1, were also significantly
upregulated at early injury response in spinal cord injury ([120]Yu
et al., 2019) (0.5∼3 h, [121]Figures 3C, 3D, and [122]S8). Most of
these seven genes were expressed in endothelial and pericyte cells
([123]Milich et al., 2021) ([124]Figure S9), including genes related
with angiogenesis (Apold1, endothelial cell-specific) and reported as
early stress response (Egr1 and Nr4a1 ([125]Shaked et al., 2015)) and
related with inflammation (Zfp36, Dusp1, Nr4a1, Junb, Klf4, and Egr1)
([126]Shaked et al., 2015; [127]Trizzino et al., 2021), suggesting the
immediate response of vascular network and inflammation around
motoneurons after injury.
Robust DEGs at 1, 3, and 7 days were further screened with either up-
or downregulation in comparison with the naive or corresponding state 1
group (upregulation was shown in [128]Figure 3B). A total of 4,916
robustly up- or downregulated genes were identified ([129]Table S4). GO
enrichment of upregulated genes ([130]Figure 3E and [131]Table S5)
showed that terms related to inflammation and immune response were the
most enriched at 3 and 7 days, such as “positive regulation of cytokine
production” and “acute inflammatory response”. We also found terms
related to myeloid cells (differentiation, migration, and homeostasis)
were significantly enriched ([132]Figure 3E). To estimate the abundance
of immune cells around a motor neuron upon injury, we employed the
CIBERSORT ([133]Newman et al., 2015) deconvolution analysis, a tool
estimated the abundances of member cell types in a mixed cell
population. The result showed a significantly increased cell fraction
of M2 macrophage at 3 and 7 days, specifically in state 2, based on a
leukocyte gene signature matrix (LM22), suggesting accumulation of
anti-inflammation macrophage around injured and regenerating
motoneurons ([134]Figures 3F and [135]S10). We found
microglia/macrophage markers showed robust upregulation at 1, 3, and
7 days, specifically at 3 and 7 days, such as Tmem119, P2ry12, Aif1
(IBA1), Csf1r, Cx3cr1 ([136]Figure 3G) and several (Tmem119, Trem2,
Aif1, Itgam (Cd11b), Adgre1 (F4/80), P2ry6, and Tyrobp (DAP12)), were
also present in another study ([137]Shadrach et al., 2021)
([138]Figure 3A). Significant robust upregulation expression of
microglia related in comparison with the naive and “intact”, suggesting
the importance of microglia-motoneuron crosstalk in axonal injury and
repair, such as genes involved in activation (Csf1, Spi1, Irf8, Laptm5,
Reg3b, etc), migration (Trem2, Csf1, and P2ry12), proliferation (Il33),
and transition from M1 to M2 (Arg1 and Igf1, [139]Figure 3H).
We also found terms related to axon regrowth were the most enriched in
“injured and regenerating” at 3 and 7 days, such as actin cytoskeleton
organization (Pdlim1, 4, and Rnd1) and microtubule polymerization
(Tubb2b, [140]Figures 3F and 3H). We found 23 RAGs with literature
evidence present at the top 100 of ranked upregulated genes, including
injured MNs released peptide/hormone (Vip, Adcyap1), cytokine (Ccl2,
Ccl7, and Tslp), growth factors/receptor (Ngfr and Grn), and
regenerating axon (Sprr1a, [141]Figure 3G). Of which, injury-induced
cytokines, Ccl2 (322-fold at 3 days and 243-fold at 7 days) and Ccl7
(788-fold at 3 days and 624-fold at 7 days), proven mostly released by
injured neurons to activate or recruit microglia/macrophage for
enhancing axon regeneration were presented at the top 10 ranked
upregulation, which was consistent with the study ([142]Kwon et al.,
2015; [143]Shadrach et al., 2021). Vip (immunoreactive vasoactive
intestinal peptide) and Adcyap1 (also known as PACAP), released by
neurons, could promote and accelerate early myelination and axon
regrowth after PNS injury ([144]Woodley et al., 2019; [145]Zhang
et al., 2002). We found low expression of Vip and Adcyap1 in the naive,
early injury (0.5 h∼1 day)), and “intact” neurons at 3 and 7 days (<1
FPKM) but dramatic increased at state 2 from 3 to 7 days (17–39 FPKM).
Adcyap1 exerts pleiotropic roles with neurotrophic factor,
neuromodulator, and neurotransmitter and could enhance neurite growth
by binding its predominant receptor, Adcyap1r1 (PAC1), in human sensory
neurons ([146]Baskozos et al., 2020).
Except for the response from glia and neurons, we also found genes
related to angiogenesis, such as Procr, Hif1a, Cd40, Grn, and Angptl4
([147]Figure 4H). Procr (Endothelial protein C receptor precursor), a
marker of endothelial stem cells, was only significantly upregulated at
state 2 (20-fold at 1 day, 42- to 47-fold at 3 and 7 days), and
upregulation was also observed in the study (9-fold at 7 days,
[148]Figure 3A) ([149]Shadrach et al., 2021). Neurovascular units are
critical to nerve regeneration and angiogenesis plays vital roles in
wound healing by neurotrophic support, such as released growth factors
(Vegfa and Vegfb) contributing a permissive microenvironment for axonal
regrowth ([150]Guaiquil et al., 2014). Injury induces hypoxia condition
that activation of Hif1a (hypoxia-inducible factor, 3.6- to 3.7-fold at
3 and 7 days), a critical regulator in regenerating neurons,
orchestrate various angiogenic factors, such as VEGF, bFGF,
angiopoietins, and angiopoietin-like (ANGPTL) proteins to enhance nerve
regeneration ([151]Cho et al., 2015; [152]Inoue et al., 2014). Thus,
robust upregulation of proangiogenic genes in endothelial cells may
provide a permissive microenvironment to promote axon regeneration.
Figure 4.
[153]Figure 4
[154]Open in a new tab
Vital TFs and core regulons in axonal injury and regeneration
(A) Three of five clusters of TF expression patterns.
(B) Heatmap plot depicts clustering among regulons inferred from SCENIC
based on connection specificity index (CSI) based on AUC score of
regulons.
(C) Frequency and log2-transformed fold-change of key regulators in M4
from [155]Figure 4A across different injured time points compared to
the naive. Regulators of regulons from M4 were shown at the bottom of
the heatmap and key regulators were shown at the right of the heatmap.
Only samples of state 3 from 1 day and state 2 from 3 to 7 days were
considered.
(D) The top 25 regulons at 1 day were identified by the master
regulator analysis algorithm (MARINa) compared to the naive. Vertical
bars indicated positions of the targets belonging to each TF’s regulon
and colors indicated the Spearman’s correlation between the expression
profiles of the TF and its targets (blue: negative correlation; red:
positive correlation).
(E) Venn diagrams depict the overlap of candidate master regulators
inferred from MARINa and SCENIC. Regulators with FDR <0.05 inferred by
MARINa at 1 day were selected. Co-operated regulators at M4 inferred by
SCENIC were selected.
(F) Expression of Arid5a at different injured time points at distinct
states. Each line depicts a sample.
(G) Binding sites between ATF3 and Arid5a were identified by the FIMO
program. Red blocks indicted ATF3-binding sites. Numbers on the x axis
represent the base number 5′ upstream from Arid5a transcription start
sites.
(H) Immunostaining of ARID5A in the ventral horn of the rat lumbar
spinal cord at different timepoints post-injury.
(I) Knockdown of Arid5a by siRNA significantly reduced the total
neurite length in primary DRG neurons. Up: representative images of
Arid5a knockdown (siRNA-Arid5a) and negative control (NC) in culture
neurons. Down: total neurite length per DRG neuron after Arid5a
knockdown was quantified. The significant difference was performed
using Student’s t-test.
(J) Enrichment score distribution of geneset (873 ortholog genes in rat
with significantly downregulated in Arid5a knockout mice independent of
IL6 stimulation) at different injured time points at distinct states
against the naive. Line colors depict injured time points and shapes
depict distinct states. Line width depicts enrichment with significant
(thick) or not (thin).
We also found terms related to ECM, such as cell-matrix adhesion and
cell adhesion mediated by integrin, were significantly enriched at 3
and 7 days ([156]Figure 3E). ECM is composed of fibrous proteins (e.g.
collagen and glycosaminoglycans) and other proteins secreted by
neighboring cells or exocytosis, which is important for guiding
regenerating axons, regulating cell behavior, and enhancing cell-cell
interaction ([157]Roumazeilles et al., 2018; [158]Yu et al., 2008).
Cell adhesion molecules, such as interacting between immune
cells-matrix (Itgam and Itgb2), neuron-matrix (Csf1 and Adcyap1), and
endothelial-matrix (Nrp1) showed robust upregulation following SNC at 3
and/or 7 days.
Dynamic TFs and inference of co-operated key TFs related to axon regrowth
A total of 345 differentially expressed TFs (47.4% of 727 TFs with
average FPKM ≥1) were identified by comparison with the naive. Of
which, 176 TFs showed robust differentially expressed after axonal
injury (frequency ≥50%, [159]Figure 4A and [160]Table S6) and were
clustered into five distinct clusters. Three clusters (k1∼3 in
[161]Figure 4A)) showed significant upregulation in the
“pre-regeneration phase” and/or “regeneration phase”. We also found
nine TFs (Arid5a, Nfil3, Atf3, Crem, Sox11, etc, [162]Figure 4A and
[163]Table S1) present in DEGs of the study ([164]Shadrach et al.,
2021). GO category of the biological process showed these TFs mainly
involved in myeloid cell differentiation and homeostasis (Spi1, Irf8,
etc), regulation of angiogenesis (12 TFs, Runx1, Hif1a, Xbp1, Stat3,
etc), neurogenesis (Cux1, Nkx6-2, Hif1a, Myc, Sox11, Etv5, etc), and
gliogenesis (Myc, Sox11, Etv5, Nkx6-2, etc), and ER unfolded protein
response (Atf3, Creb3, Xbp1, etc) ([165]Table S7).
To screen co-operated key TFs or other regulators (e.g. epigenetic
modification-related genes), we performed regulatory network analysis
using the SCENIC pipeline. Given a relatively high correlation in
expression of DEGs between rat and mouse post-SNC ([166]Figure 3A) and
lacking cis-regulatory motifs for the rat, we inferred regulon (a
regulator and its targets) activity based on mouse motif dataset
(−500–100 bp) and then performed clustering analysis of connection
specificity index (CSI) ([167]Suo et al., 2018). The regulons could be
divided into seven distinct modules and two modules (M4 and M7) were
focused with pro-growth TFs hits, specifically M4 (Atf3, Stat3, Myc,
Jun, Smad1, Klf6, and Klf7, [168]Figure 4B). Master regulators of 16
regulons (M4) were robust DEGs, respectively ([169]Figure 4C). GO
enrichment showed M4 regulons were significantly enriched in terms
related to cytoskeleton, neurogenesis, signaling transduction,
inflammatory and immune response, and response to wounding
([170]Table S8). Except for well-known regulators (Atf3, Stat3, Jun,
etc) in M4, we also found TFs related to unfolded protein response
(UPR) were enriched (Atf3, Creb3, and Xbp1). ER stress could be
provoked by tissue injury. UPR is important for preserving ER
proteostasis and consists of three branches (IRE1α-XBP1, PERK-CHOP, and
ATF6-SREBP) depending on their activation in either promoting protein
folding or apoptotic process or sterol biosynthesis. We observed
upregulation of Xbp1 at 3 and 7 days post-SNC, but Ern1 (also known as
IRE1α) showed no significant difference. Creb3, a homolog with Atf6,
could be retrograded transported to the nucleus to regulate the UPR and
cholesterol biosynthesis by binding with Srebf1 or Srebf2 ([171]Ying
et al., 2014). We found Creb3 and Srebf2 were co-occurred in M4 and
robust upregulation at 3 and/or 7 days post-SNC. Besides, enrichment
analysis of UPR-related terms showed PERK-mediated UPR was the most
enriched from 1 to 7 days ([172]Figure S11), but its downstream
regulator (CHOP) has no significant difference. These results suggested
Creb3-Srebf2 pathway may be the predominant UPR in spinal motor neurons
after PNI. We also found epigenetic modification-related regulators
occurred in M4 (Hdac2 and Gadd45a). “Chromatin remodeling” was
significantly enriched at 3 and 7 days post-SNC ([173]Figure S12).
Strikingly, we found a chromatin-remodeling enzyme, Apobec1, a known
DNA demethylation mediator, showed a dramatic increase at 3 (74.2-fold)
and 7 days (60.5-fold) post-SNC. We also found other seven known DNA
demethylation mediators ([174]Weng et al., 2017) (Tet3, Tet2, Gadd45a,
b, g, and Apobec1 and 3, [175]Figure S12) showed significant
upregulation, specifically Apobec1, 3 and Gadd45g. Apobec1 was solely
expressed in microglia and knockdown of Apobec1 led to aberrant
myelination, increased inflammation, and lysosomal anomalies in mouse
brain ([176]Cole et al., 2017). Thus, dramatic upregulation of Apobec1
may be beneficial for axon regeneration by maintaining a balance
between the homeostatic and activated immune functions of microglia.
We also employed the Algorithm for the Reconstruction of Accurate
Cellular Network (ARACNe) which could identify transcriptional
interaction based on information-theoretic method on a large gene
expression data ([177]Margolin et al., 2006). We analyzed 1,595 TFs or
genes with “DNA-binding transcription factor activity” and then
employed ssMARINa to infer candidate master regulators (MRs) post-SNC
([178]Aytes et al., 2014). We focused on activated MRs at 1 day which
is critical on axon regeneration as we extensively discussed above.
MARINa identified 41 candidate MRs (FDR <0.05, [179]Table S9) at 1 day
and top 10 ranked MRs including Atf3, Klf6, Myc, and Sox11, which were
well-known pro-regenerative MRs in axon growth ([180]Venkatesh et al.,
2021) ([181]Figure 4D). We found 13 out of 20 co-operated regulators in
M4 were also present in MRs identified by MARINa ([182]Figure 4E).
Specifically, we noticed the top MR, Arid5a (AT-Rich interaction domain
5A), has dual roles as TF and RNA-binding protein, that previously
reported induction by inflammation stimulation to translocate from
nucleus to the cytoplasm, stabilizes the expression of inflammatory
genes, such as Il6 and Stat3 ([183]Masuda et al., 2013; [184]Nyati
et al., 2019). We found the expression of Arid5a was only significantly
upregulated at state 2 and/or state 3 (1 day: 7.69-fold at state 3 and
9.61-fold at state 2; 3 days: 8.58-fold; 7 days: 10-fold,
[185]Figure 4F), showing the activation at 1 day. Additionally, we
found Arid5a was one target of the MR, Atf3, which presented a binding
site of Atf3 motif within 1k bp from 5′ upstream of transcription start
site ([186]Figure 4G) and was also supported by ChIP data
([187]Gilchrist et al., 2006). We found ARID5A mostly expressed in
spinal neurons and obviously upregulated in injured MNs at 1 day and
gradually increased at 3 and 7 days post-injury, specific in the
cytoplasm ([188]Figure 4H). This trend was consistent with our
sequencing data. Furthermore, reanalysis of a single-cell sequencing
([189]Hu et al., 2016) of sensory neurons and translatome datasets
([190]Rozenbaum et al., 2018) of sensory neurons and glia in mouse DRG
after sciatic nerve injury also showed robust upregulation of Arid5a in
injured sensory neurons ([191]Figure S14), which was consistent with
the findings in MNs. These evidences indicated upregulation of Arid5a
was conserved in injured neurons against sciatic nerve injury.
Knockdown of Arid5a in primary sensory neurons showed a significant
reduction of total neurite length in vitro ([192]Figures 4I and
[193]S15, p < 0.001), suggesting its new role in axon regeneration.
Reanalysis of recently released RNA-seq dataset of mouse embryonic
fibroblast cells (MEFs, could be induced into neuronal cells
([194]Treutlein et al., 2016)) from wild type or Arid5a^−/− with or
without IL6 stimulation ([195]Nyati et al., 2021) showed that 987 genes
(873 genes were orthologs in rat) were significantly downregulated in
Arid5a^−/− while less change induced by IL6 compared to the wild type
([196]Figure S13). We termed these genes as Arid5a-dependent and the
top 20 downregulated genes included Stmn2 (also known as SCG10, a
marker of regenerating axons and regulation of microtubule dynamics).
GO enrichment of Arid5a-dependent genes showed the term “regulation of
neuron projection development” was significantly enriched, including
genes related to axon regeneration (Gap43, Stmn2, and Stmn4,
[197]Figure S13). GSEA enrichment of Arid5a-dependent genes in our data
showed significant enrichment in state 2 (injury and regenerating
state) from 3 to 7 days post-injury ([198]Figure 4J).
An axon regeneration-associated co-expressed module
Next, we performed weighted co-expression analysis (WGCNA) to detect
regeneration-associated modules and identified 19 distinct modules
([199]Figure 5A). We found the “turquoise” (1,705 genes,
[200]Table S10) showed a significant correlation with state 2 from 3 to
7 days based on the module eigengene (adjust p-value < 0.05,
[201]Figure 5A). Genes in “turquoise” showed a high expression in state
2 ([202]Figure 5B). We compared this module with the “regenerative”
module (413 genes) in DRG identified by [203]Chandran et al. (2016) and
36.9% were also present in our “turquoise” module, such as Ccl2, Ccl7,
Atf3, Cd44, Il6, and Jun. Next, we used a ternary plot-based method to
explore expression patterns of genes in the “turquoise” module. Genes
in the “turquoise” module were classified into seven major clusters and
then focused on two clusters with upregulation after SNC (cluster 1 and
6, [204]Figure 5C). We found most RAG genes were present in cluster 1
(Ngfr, Pdlim1, Gap43, Mmp2, Tspo, etc) and cluster 6 (Atf3, Il6, Myc,
Jun, Socs3, Sox11, Hif1a, Klf6, etc). GO enrichment analysis of genes
in these two clusters showed inflammation and immune response were the
most enriched categories, such as “microglial cell activation”,
“neuroinflammatory response”, and “positive regulation of cytokine
production” ([205]Tables S11 and [206]S12), suggesting that sterile
inflammation induced by injury is closely related with axon
regeneration. In addition, we found TF-related terms (regulation of DNA
binding) were significant enrichment in cluste r6, including Bcl3,
Jak2, Jun, Myc, and Sox11, supporting our above hypothesis that
pro-growth TFs were activated as early as 1 day post-SNC. In addition,
pro-regenerative pathways, such as MAPK cascade was also enriched at
1 day post-SNC. Genes related with “actin filament organization” and
“axon regeneration” were enriched at stage 3 (namely 3ays and 7days),
such as Gap43, Tspo, and Rac2, suggesting that axon regeneration was
mainly initiated as early as 3 days.
Figure 5.
[207]Figure 5
[208]Open in a new tab
A module potentially related to axonal injury and regeneration was
revealed by co-expression analysis
(A) Relationship between modules and traits (injured time points with
different states here). Only correlation with significant (adjust
p-value ≤ 0.05) was shown (∗).
(B) Gene expression pattern of the “turquoise” module based on the
module eigengene values.
(C) Ternary plots showed relative expression abundance of 1705 triads.
Each circle represents a gene triad with Stage 1, 2, and 3 coordinates
consisting of the relative contribution of each gene to the overall
triad expression. Seven categories were assigned based on triad values.
(D and E) Gene expression trends of cluster 1 (D, left) and cluster 6
(E, left) and GO enrichment of biological processes. N indicates gene
number and p indicates p-value.
(F) Highly connected protein-protein interaction of genes from clusters
1 and 6. Nodes with degrees
[MATH: ≥20 :MATH]
were shown.
Next, we inferred the potential protein-protein interaction (PPI)
network of genes in clusters 1 and 6. A highly connected PPI network
showed the importance of proteins related injured
neuron-microglia/macrophage crosstalk (Il6, Csf1, Ccl2, and Csf1r) with
high connection and dramatic transcriptional changes (degree
[MATH: ≥20 :MATH]
, [209]Figure 5F). We found genes involved in microglia response
(activation and proliferation) were present in the network, such as
Trem2, Tyrobp, Spi1, Irf8, Cx3cr1, and Aif1. We also found the
pro-regenerative TFs, Jun, Atf3, Myc, and Hif1a, were present in the
network. Potential genes identified here will provide important
candidates for further drug development for PNS treatment.
Discussion
Intrinsic regeneration capacity and exogenic supportive
microenvironment greatly contribute to successful axon regeneration,
specific immune system ([210]Scheib and Hoke, 2013). Taking advantage
of LCM-seq featured with collecting the target cell and its
microecology in situ and in-depth sequencing, enabled us to obtain
temporal-spatial profiles of MNs and the in vivo interactions between
injured MNs and non-neuronal cells within the context of axon injury
and regeneration through thorough molecule analysis. To our current
knowledge, this new dataset was the first, massive and in-depth
sequencing at approximate single-MN resolution. To facilitate
exploration by other researchers, our data are easily searchable at
MNniches-seq.
In this study, we established dynamic expression of MN niches at hours
(0.5, 3, and 12 h) or days (1, 3, and 7 days) following SNC and the
naive and divided into three response stages via thorough molecules
analysis (0.5–12 h: early stress response phase; 1 day:
pre-regeneration phase; 3∼7 days: regeneration phase, [211]Figure 6).
Comparison of DEGs at 7 days post-SNC between our study and the study
([212]Shadrach et al., 2021) demonstrated significant correlations and
together with immunohistochemical validation indicating reproducibility
and accuracy of our data. Based on clustering analysis and Atf3, we
could well discriminate samples with injured MNs or “intact” MNs at 1,
3, and 7 days, which were supported by immunostaining results that
ATF3-positive MNs present as early as 1 day of part MNs. Besides, the
estimated injured ratios were in accordance with the estimation in
previous studies ([213]Holland et al., 2019; [214]Renthal et al.,
2020). Two cell states (state 1: “intact” and state 2: “injury”) were
distinguished at 3 and 7 days, which was perfectly supported by both
Atf3 expression and clustering results. Strikingly, in addition to
“intact” and “injury” states, we also found an additional state
(defined as state 3, 26.5% of cells at 1 day, [215]Figure 2A) that
transient dramatic upregulation of Atf3 indicating injury which
represents predominate state at 1 day, while “injury” states were
rarely captured (defined as state 2, 6.1%), which suggested that state
3 at 1 day was intermediate transition states and rarely dissected MNs
in state 2 may be explained by individual fast growth MNs. Considering
the small number of MNs in state 2 at 1 day, we majorly focused on
state 3 at 1 day. DEG analysis compared with the naive showed clear and
distinct gene profiles after sample separation suggested that this
could well reflect cell states at 1, 3, and 7 days post-SNC. Evidences
of regenerating axon marker, RAGs, enrichment of multiple
pro-regenerative pathways showed early pro-regenerative and master
regulators (e.g., Atf3, Jun, Sox11, Smad1, and Stat3) were activated as
early as 1 day while regenerating RAGs (e.g. Gap43 and Sprr1a) were
activated in samples with injured MNs at 3 and 7 days, suggesting
“pre-regenerative” at 1 day and “regenerative” at 3 and 7 days. The
small number of DEGs at 0.5–12 h and significant increase of DEGs at
1 day (state 2 and state 3) and dramatic increase of DEGs at 3 and
7 days (state 2) suggest 1 day post-injury is a key point. Though the
increase of DEG number in state 1 at 3 and 7 days, inflammation and
immune response-related genes were the most enriched which provide a
vital control to explore crucial genes that contributed to the injury
response and regeneration in injured MNs.
Figure 6.
[216]Figure 6
[217]Open in a new tab
Schema of three response phases of injured motoneuron soma post-SNC
1d was a key point in axon regeneration that activated essential
pro-regrowth-related TFs and non-neuronal cells. High dynamic of
injured motoneurons niches occurred at regeneration phase (3 and
7 days) with high cellular heterogeneity and molecules involved in
regeneration.
By comparison with samples containing the naive and corresponding
“intact” MNs, we found genes involved in crosstalk between injured
neuron-microglia/infiltrated macrophage was ranked at the top of DEGs,
such as neuronal Il6, Ccl2, Ccl7, Csf1, and Reg3b, which was not
different between samples with “intact” MNs and the naive. Significant
upregulation of microglia/macrophage markers, e.g., Aif1 (IBA1), Lyz2,
Clec7a, Trem2, Tyrobp, Adgre1, Itgam, Ly86, and C1qa, in samples with
“injury” MNs at 3 and 7 days compared to that with the naive and
“intact” MNs, suggested closely crosstalk between injured
neuron-microglia, which was also present in the DEGs of the study
([218]Shadrach et al., 2021). Our immunostaining of ANXA3 and IBA1
demonstrated that microglia almost adhered around injured MNs at 7 days
post-SNC, which was in accordance with our sequenced result and other
study ([219]Rotterman and Alvarez, 2020). Single-cell sequencing of the
injured spinal cord showed the significant microglia subtypes
transition from homeostatic microglia subtype to disease-associated
microglia subtype at 1 day then revert back to a more homeostatic state
(3 and 7 days) ([220]Milich et al., 2021), suggesting
microglia/macrophage contribute the most to the dynamic
microenvironment. Activated microglia/macrophage released neurotrophic
factors but also produced pro-inflammatory factors contributing to
phagocytose cell debris for promoting regeneration ([221]Scheib and
Hoke, 2013). M2 phenotype of microglia/macrophage could be
pro-regenerative by releasing anti-inflammation factors ([222]Scheib
and Hoke, 2013). We found significant upregulation of genes promoting
microglia/macrophage phenotype transition, e.g., Arg1. In addition to
microglia/macrophage, we also found robust upregulation markers of
endothelial cells (Procr), astrocyte (Gfap). We found enriched genes
involved in angiogenesis included Hif1a, Cd40, Grn, and Angptl4. In
non-cellular, we also found the dynamic extracellular matrix
organization reflected by DEGs involved cells-matrix (integrins, such
as Itgam and Itgb2). In addition to releasing factors mediating the
interaction of neuron-cells, neuron-matrix, and cells-cells for
permissive microenvironment, injured MNs also activate intrinsic
pro-regenerative TFs, such as Atf3, Sox11, Jun, Smad1, and Myc.
Besides, injured MNs released pro-regenerative peptides/hormones and
growth factors/receptors, such as Vip, Adcyap1, Ngfr, Grn, and
anti-apoptotic proteins (e.g., Gadd45a) and upregulate axon regrowth
factor (e.g., Gap43, Sprr1a, Rac2, and Pdlim1).
We also deeply explored potential regulators, including TF and
epigenetic modification (e.g., demethylation and acetylation). We
identified potential co-operated regulators including well-known TFs
(Klf6, Klf7, Atf3, Myc, Stat3, and Smad1) and unfolded protein response
(Atf3, Creb3, and Xbp1) and epigenetic modification (Hdac2 and
Gadd45a). Except for the well-known TFs supported by MARAINa analysis,
we also highlighted a candidate master regulator, Arid5a. Its roles in
the regulation of chondrocyte-specific transcription (acting as a
transcription factor) and inflammatory processes (acting as an
RNA-binding protein) by stabilizing inflammatory mRNA transcripts were
described ([223]Amano et al., 2011; [224]Masuda et al., 2013). However,
studies of Arid5a expression and potential roles during axonal injury
and regeneration were blank. We found significant upregulation of
Arid5a in samples with injured MNs both at transcriptional and protein
level post-SNC ([225]Figures 4G and 4I). This was also supported by
significant upregulation of translating Arid5a level in injured MNs
compared to the control ([226]Shadrach et al., 2021). Knockdown of
Arid5a in primary sensory neurons suggested its potential new role in
axon growth ([227]Figure 4I). Additionally, Arid5a is also a potential
target of the pro-regenerative regulator, Atf3, and showed a high
correlation in expression ([228]Figure 4A, Pearson’s cor: 0.98).
Reanalysis of Arid5a knockout (KO) dataset of mouse embryonic
fibroblast cell showed axon growth-related genes were largely
downregulated in Arid5a KO group and geneset enrichment of these
downregulated genes showed significantly enriched at samples with
injured and regenerating MNs ([229]Figure 4J), indicating the
importance of Arid5a in axon regeneration. Recently, a study of brain
ischemic injury showed that overexpression of ARID5A in primary
cortical neurons elicited neuroprotection against oxygen-glucose
deprivation injury ([230]Cai et al., 2019).Taken together, we suspected
that upregulation of Arid5a in samples with injured neurons may promote
axon regeneration and might be an undescribed role of the less studied
function of Arid5a in PNI, and its mechanism required further
investigation.
In conclusion, we firstly present the landscape of cellular and
molecular interactions between injured and regenerating MNs and its
permissive microenvironment in situ during the complex first seven days
post-injury using LCM-seq. Specifically, distinguishing and comparison
of “intact” and “injury” MN niches at the similar environment yielding
numerous injury-induced and pro-regenerative genes, though most of
their function in promoting axon regeneration of MNs remains to be
explored. Through thorough molecular analysis, we resolved three
response stages and molecular characteristics of the injured MN and its
microenvironment and highlighted critical regulators, for instance,
Arid5a. This study not only provides gene regulation patterns for
regeneration of the injured spinal motoneuron and its dynamic
microenvironment but also gains insights into potential new targets for
the treatment of motoneuron dysfunction and gives the genetic basis for
development of tissue engineering and regenerative medicine.
Limitations of the study
Two major limitations exist in this study. First, acting as a method
for single-cell isolation, LCM has a major disadvantage that inevitably
contains neighboring cells or non-cells of the target cell. Thus, the
expression profile presented here not only contained the motoneuron but
also mixed with its microecology, including other cells or non-cells.
Second, for hypothetical roles in axonal injury and regeneration and
underlying mechanisms of candidate genes, for example, Arid5a, further
functional experiments in vitro and in vivo are required.
STAR★Methods
Key resources table
REAGENT or RESOURCE SOURCE IDENTIFIER
Antibodies
__________________________________________________________________
Anti-Choline Acetyltransferase antibody Abcam Cat# ab137349,
RRID:[231]AB_2754985
Anti-IBA1 antibody Abcam Cat# ab153696, RRID:[232]AB_2889406
Anti-ARID5A antibody Sigma-Aldrich Cat# SAB2100148,
RRID:[233]AB_10598838
Annexin A3 antibody Proteintech Cat# 11804-1-AP, RRID:[234]AB_2057455
ChAT Polyclonal antibody Thermo Fisher Scientific Cat# PA5-29653,
RRID:[235]AB_2547128
ATF3 antibody Abcam Cat# ab58668, RRID:[236]AB_879578
NeuN antibody Sigma-Aldrich Cat# MAB377
Anti-Beta-III Tubulin antibody R & D Systems Cat# MAB1195,
RRID:[237]AB_357520
__________________________________________________________________
Chemicals, peptides, and recombinant proteins
__________________________________________________________________
Diethyl pyrocarbonate Sigma-Aldrich Cat# D5758-100ML
Collagenase type I Sigma-Aldrich Cat# SCR103-250MG
Neurobasal-A medium Thermo Fisher Scientific Cat# 10888022
B-27 Supplement Thermo Fisher Scientific Cat# 17504044
L-Glutamine Thermo Fisher Scientific Cat# 25030149
Penicillin-Streptomycin Beyotime Cat# C0222
poly-L-lysine Sigma-Aldrich Cat# P4832
Trypsin-EDTA Gibco Cat# 25200-072
__________________________________________________________________
Critical commercial assays
__________________________________________________________________
Lipofectamine RNAiMAX reverse transfection reagent Thermo Fisher
Scientific Cat# 13778075
__________________________________________________________________
Deposited data
__________________________________________________________________
Raw and analyzed data This paper CNGBdb: CNP0002605
Single-cell sequencing of mouse spinal cord after injury [238]Milich
et al. (2021) GEO: [239]GSE162610
Bulk RNA-seq of rat spinal cord after injury [240]Yu et al. (2019)
NODE: OEP000369
Mouse embryonic fibroblast cells (MEFs) in wild-type or Arid5a knockout
[241]Nyati et al. (2021) GEO: [242]GSE178250
Single-cell sequencing of mouse DRG neurons [243]Hu et al. (2016) GEO:
[244]GSE71453
RiboTag data from DRG ganglia from Isl1-Cre and Dhh-Cre RiboTag mice
[245]Rozenbaum et al. (2018) GEO: [246]GSE102316
RiboTag data of motoneurons from RiboTag; ChATCre mice [247]Shadrach
et al. (2021) GEO: [248]GSE162028
__________________________________________________________________
Experimental models: Organisms/strains
__________________________________________________________________
Sprague-Dawley rat (SPF) Experimental Animal Center of Nantong
University SCXK(Su)2019-0001
__________________________________________________________________
Software and algorithms
__________________________________________________________________
ImageJ [249]Schneider et al. (2012) [250]https://imagej.nih.gov/ij
NeuronJ-1.4.3 [251]Meijering et al. (2004)
[252]http://www.imagescience.org/meijering/software/neuronj/
R for statistical computing-4.1.0 [253]R Core Team (2018)
[254]https://www.r-project.org
Seurat-4.1.0 [255]Butler et al. (2018)
[256]https://satijalab.org/seurat/
clusterProfiler-4.0.5 [257]Yu et al. (2012)
[258]https://bioconductor.org/packages/release/bioc/html/clusterProfile
r.html
DESeq2-1.32.0 [259]Love et al. (2014)
[260]https://bioconductor.org/packages/release/bioc/html/DESeq2.html
WGCNA-1.70 [261]Langfelder and Horvath (2008)
[262]https://horvath.genetics.ucla.edu/html/CoexpressionNetwork/Rpackag
es/WGCNA/
ARACNe2 [263]Margolin et al. (2006)
[264]http://califano.c2b2.columbia.edu/aracne
ssMARINa-1.01 [265]Aytes et al. (2014)
[266]https://figshare.com/articles/dataset/ssmarina_R_system_package/78
5718
GSEA-4.1.0 [267]Subramanian et al. (2005)
[268]https://www.gsea-msigdb.org/gsea/index.jsp
FIMO-5.4.1 [269]Grant et al. (2011)
[270]https://meme-suite.org/meme/tools/fimo
pySCENIC-0.11.2 [271]Aibar et al. (2017)
[272]https://github.com/aertslab/pySCENIC
AUCell-1.14.0 [273]Aibar et al. (2017)
[274]https://bioconductor.org/packages/release/bioc/html/AUCell.html
HISAT2-2.0.4 [275]Kim et al. (2015)
[276]http://daehwankimlab.github.io/hisat2/
HTseq-0.6.0 [277]Anders et al. (2015)
[278]https://github.com/htseq/htseq
CIBERSORT [279]Newman et al. (2015)
[280]https://cibersort.stanford.edu/
[281]Open in a new tab
Resource availability
Lead contact
Further information and requests for resources and reagents should be
directed to and will be fulfilled by the lead contact, Jian Yang
([282]dna2009@ntu.edu.cn).
Materials availability
This study did not generate new unique reagents.
Experimental model and subject details
Animals
Adult Sprague-Dawley (SD) rats (male, weighted ∼200g, provided by the
Experimental Animal Center of Nantong University) were housed under
standard conditions with a 12-h dark/light cycle and had free access to
food and water. The rats were randomly divided into seven groups
(naive, injured timepoints: 0.5, 3, 12h, 1, 3ays, and 7days) to be
deeply anesthetized by an intraperitoneal injection of a compound
anesthetic (chloral hydrate 4.25 g, magnesium sulfate 2.12 g, sodium
pentobarbital 886 mg, ethanol 14.25 mL, and propylene glycol 33.8 mL in
100 mL) at a dose of 0.2–0.3 mL/100 g.
Method details
Nerve crush
The bilateral sciatic nerve crushed models were conducted as follows:
the sciatic nerve was lifted from the right hindlimb on the lateral
aspect of the mid-thigh and a 1-cm segment was clamp-injured with a
sterile forceps for 30 s. After crush injury, the muscle was then
sutured. One group was euthanized immediately and the other groups at
hours, 0.5, 3, 12h, and days 1, 3, 7 were euthanized after surgery. All
animal experiments in this study were performed in accordance with the
guidelines for animal care and were approved by the Administration
Committee of Experimental Animals of Jiangsu.
Tissue preparation and rapid staining of MNs before laser capture microscopy
All dissection tools and work surfaces were cleaned with RNase ZAP
wipes. Animals were anesthetized by the compound mentioned above and
sacrificed by decapitation. The lumbar intumescence part (L4 and L5) of
the spinal cord was rapidly removed, immediately embedded in a
pre-cooled OCT compound on dry ice to avoid thawing of the tissue. Thin
coronal sections (30 μm) were cut at −20°C using a cryostat (Leica
CM3050 S) and placed onto PEN Membrane Slide (Leica). Slides were
allowed for complete dehydration in the dry box for 10 min at −20°C.
Slides containing tissue sections were incubated with Alexa Fluor® 488
Anti-Choline Acetyltransferase antibody (ChAT,1:50 dilution, Abcam) for
20 min at 4°C, in dark. After being washed with DEPC water, keep the
slides into the dry box at 4°C for 20 min until they were dried. All
steps were performed on ice and all aqueous reagents were freshly
prepared with DEPC water for each experimental day.
Laser capture microscopy
Preheat the Leica LMD6500 Laser Microdissection System and shut on the
CTR6500 controller. Placed slides into the slide holder and captured
neurons of anterior crus in the spinal cord which were labeled green
signal. To avoid damage MNs and ensure a single MN was contained as
possible, we cut circular outlines with 1.5–2 times the diameter of
obvious MN morphology, and MNs with a visible nucleus in the ventral
horn of the spinal cord were dissected. Cells were cut at 40
[MATH: × :MATH]
magnification while keeping laser power to a minimum. The Leica system
uses gravity to collect cells, which fall into the cap of the
collection tube (0.2 mL RNase free PCR tubes) that is placed directly
under the slide. These tubes were added 5 μL lysis buffer acquired from
Annoroad Gene Technology Co. Ltd. company in advance. Samples were
centrifugated in a table centrifuge for 10 s. The tubes were carefully
sealed with parafilm (Pechiney Plastic Packaging), labeled, and
snap-frozen on dry ice. The whole process from thawing the slides until
sample freezing never took longer than 3 h in any of the experiments.
Samples were kept in collection tubes at −80°C until further
processing.
Smart-seq2 procedure for cDNA library preparation and sequencing
The cDNA sequencing library was prepared as follows the Smart-seq2
protocol ([283]Picelli et al., 2013), including cell lysis, cDNA
synthesis, PCR amplification and purification, and quality control.
Qualified libraries were then sequenced using the Illumina HiSeq Xten
platform. Library construction and sequencing in this study were
conducted by Annoroad Gene Technology Co. Ltd. Company.
Immunofluorescence histochemistry
Under anesthesia, the animals were perfused with normal saline and 4%
paraformaldehyde. Lumbar intumescence parts of the spinal cords were
excised. The samples were post-fixed in 4% paraformaldehyde for 8 h,
dewatered in 10%–20%–30% sucrose, and cut on a cryostat into 12 μm
thick sections. The immunofluorescent staining was performed using
rabbit anti-IBA1 antibody (1:1000 dilution, Abcam), mouse anti-ANXA3
antibody (1:1000 dilution, 11804-1-AP, Proteintech), rabbit anti-ARID5A
antibody (1:100 dilution, SAB2100148-100UL, Sigma-Aldrich), mouse
anti-ATF3 antibody (1:100 dilution, ab58668, Abcam), rabbit anti-ChAT
antibody (1:250 dilution, PA5-29653, Invitrogen), and mouse anti-NeuN
antibody (1:250 dilution, MAB377, Sigma) and Hoechst 33342 (1:3000
dilution, Abcam), respectively. The samples were incubated with the
primary antibody at 4°C for 12 h, followed by further reaction with the
secondary antibody (Goat anti-Mouse IgG-Alex-488, 1:500; Goat
anti-Rabbit IgG-Cy3, 1: 100) at room temperature for 1 h. At last,
images were taken under fluorescence microscopy (Leica, Germany).
DRG neuron culture and siRNA transfection
DRGs were dissociated from adult male Sprague-Dawley rats (200-220g,
Provided by Experimental Animal Center of Nantong University, Nantong,
China). DRGs were digested with 0.1% collagenase type I (Sigma, St
Louis, MO, USA) in the oxygenated DMEM (Corning, New York, USA). After
120 min incubation at 37°C, DRGs were transferred into 0.25% trypsin
(Thermo Fisher Scientific, Waltham, MA, USA) for more 15 min at 37°C.
Then, single-cell suspension was centrifuged for 5 min at 1000 r/min
with 15% BSA. After chucking the supernatant, the cells were
resuspended and plated with Neurobasal-A medium (Thermo Fisher
Scientific, Waltham, MA, USA) enriched with L-Glutamine (1%. Thermo
Fisher Scientific, Waltham, MA, USA), B27 supplement (2%, Thermo Fisher
Scientific, Waltham, MA, USA), and penicillin-streptomycin (1%,
Beyotime, Shanghai, China). Experimental plates were pre-coated with
poly-L-lysine (Sigma, St Louis, MO, USA). During the plating, DRG
neurons were transfected with siRNAs (OBiO, Shanghai, China) using
Lipofectamine RNAiMAX reverse transfection reagent (Thermo Fisher
Scientific, Waltham, MA, USA) according to product instructions. After
transfection, the culture medium was changed 16 h later. The sequence
of siRNA-Arid5a was sense 5′-CCACUGAGAAGCUGAAGAATT-3′, antisense
5′-UUCUUCAGCUUCUCAGUGGTT-3′. The sequence of scramble siRNA was sense
5′-UUCUCCGAACGUGUCACGUTT-3′, antisense 5′-ACGUGACACGUUCGGAGAATT-3′.
After cultured purified DRG neurons were transfected with siRNA for 16
h, the medium was changed to remove the redundant siRNA. Then, after
another 48 h of culture, the DRG neurons were resuspended with 0.05%
Trypsin-EDTA (Gibco, Carlsbad, CA, USA) for 1.5 min at 37°C and
replanted to experimental plate. 18 h after replant, DRG neurons were
fixed with 4% PFA. Neurites were labeled by immunostaining of β-tubulin
(Tuj1, 1:500, R&D, Minneapolis, MN, USA), and its immunofluorescence
was acquired by Zeiss microscopy (Zeiss, Axio Imager M2m). Total
neurite length (only neurite length greater than soma diameter was
considered) of each neuron was measured using ImageJ ([284]Schneider
et al., 2012) (Neuron J ([285]Meijering et al., 2004), v1.4.3).
Approximately 100 cells were analyzed for each condition.
Quantification and data analysis
RNA-seq mapping, quantification, and clustering analysis
Rat genome and annotation were retrieved from the Ensembl database
(v89). After trimming low-quality reads, clean reads were mapped
against the reference genome with default setting using HISAT2
([286]Kim et al., 2015). HTseq package ([287]Anders et al., 2015) was
used to obtain gene count and normalized into FPKM by considering
sequencing depth and gene length. Samples with low expression of Chat
(FPKM<1) were discarded in the downstream analysis and 682 samples were
kept. Sample correlation was estimated based on Pearson’s correlation
of log2-transformed FPKM. To estimate neuron signature activity, we
used the AUCell package to calculate the AUC score of neuronal markers.
To investigate sample distribution, we performed Principal Component
Analysis (PCA) based on log-transformed FPKM using the R program
([288]R Core Team, 2018). UMAP clustering was performed using the
Seurat.
Differential expression analysis
DESeq2 package ([289]Love et al., 2014) was employed to detect
differential expression genes (DE-Gs). Due to the small DEG number in
the early injury stage, we adopt a relatively relaxed threshold that
considers both significance and fold-change. Genes with p-value
[MATH: ≤0.05 and |log2fold−change|<
/mo>≥log21
.5 :MATH]
were defined as DE-Gs compared to the control (naive or corresponding
state1). Considering varied sample sizes among different timepoints
after sample separation, we further used a bootstrap strategy to detect
robust and reliable DE-Gs according to the study ([290]Aguila et al.,
2021). Simply as follows.
* a)
Sample sizes were assigned as 3,4,5,6,10, 15/20, …n, where n is the
minimal sample size of the investigated group and the control
group. Internal of 5 was used at times points (1, 3ays, and 7days)
and internal of 10 was used at timepoints (0.5, 3, 12h).
* b)
For each sample size, we randomly chose samples from two groups
respectively. Then performed differential expression analysis with
DESeq2 package and screened DE-Gs.
* c)
Repeated b) for 100 times
* d)
The frequency of DE-Gs was calculated.
To define robust DE-Gs, we adopt the threshold, average frequency≥ 50%.
Transcriptional factors (TFs) analysis
SCENIC pipeline ([291]Aibar et al., 2017) was used to perform
regulatory network inference and mouse RisTarget motif database
(−500–100 bp) was employed. Critical regulons were screened by
clustering analysis of connection specificity index ([292]Suo et al.,
2018) (CSI) based on area under the recovery curve (AUC) score of
regulons (AUC score ≥0.1 were considered). ARACNe ([293]Margolin
et al., 2006) was also used to reconstruct gene regulatory networks,
which was based on an information-theoretic method on large gene
expression data. Then, master regulons were inferred by the ssMARINa
method ([294]Aytes et al., 2014). ATF3 motif (accession: MA0605) was
retrieved from JASPAR database ([295]Fornes et al., 2020) and FIMO
([296]Grant et al., 2011) was used to identify binding sites of Atf3 in
Arid5a within 1k bp upstream from the transcription start sites.
Gene set enrichment analysis
Rat geneset were retrieved using the msigdbr package. GSEA
([297]Subramanian et al., 2005) was used to perform enrichment
analysis. Terms with FDR q-value ≤ 0.25 and NOM p value ≤0.05 was
defined as a significant enrichment.
WGCNA co-expression analysis
WGCNA package ([298]Langfelder and Horvath, 2008) was applied to
construct a single co-expression network and inferred important modules
related to axonal injury and regeneration. Firstly, low expressed genes
were filtered with row average FPKM <1. An automatic, one-step network
construction method was employed using the default parameters except
for the power of β = 4. The expression profile of each module was
estimated based on the first principal component (module eigengene,
ME). The correlation and significance between modules and cell states
at each injured timepoints were calculated using the function “cor” and
“corPvalueStudent” and adjusted p-value was calculated using the
function “p.adjust”. The ternary plot-based method could be used for
defining expression patterns of triads ([299]Ramirez-Gonzalez et al.,
2018). We employed the ternary to explore expression patterns of genes
from “turquoise” module along three injured stages (stage1: 0.5, 3,
12h; stage2: 1day from state3; stage3: 3ays and 7days from state2). For
each gene from “turquoise”, the average expression were normalized
against the total expression within the triad to represent the relative
expression. The R package “ggtern” was used for the ternary plot
([300]Hamilton and Ferry, 2018). Next, we divided these triads into
seven clusters by “20, 20 and 60” rule. For example, if a triad with
relative expression in stage1 and stage2 less than 20 while the
relative expression in stage3 large than 60, then this triad belongs
“Stage3 dominant”.
Functional enrichment analysis
ClusterProfiler package ([301]Yu et al., 2012) was used to perform Gene
Ontology (GO) and KEGG pathway enrichment analysis.
Public dataset analysis
Time-series course bulk RNA-seq of a long-term investigation of gene
expression after rat spinal cord injury (hemi-transection) were
retrieved from our previous study ([302]Yu et al., 2019). For the
single-cell RNA-seq of spinal cord injury dataset, we retrieved the
expression matrix from the GEO database ([303]GSE162610) and performed
standard analysis (normalization and standardization) with Seurat
([304]Butler et al., 2018) (V3.0). Original cell type assignment
([305]Milich et al., 2021) was used and expression of key genes was
visualized using Seurat. For stranded RNA-seq of wild-type and
Arid5a^−/− mouse MEFs with or without IL6 stimulation were retrieved
from the GEO database ([306]GSE178250) and were reanalyzed. DE-Gs
identified by comparison expression of Rpl22^HA IP in injured and
regenerating MNs with the uninjured in mouse spinal cord at 7days
post-SNC were retrieved from the supplement ([307]Shadrach et al.,
2021). Datasets of single-cell sequencing of DRG sensory
([308]GSE71453) and translatome ([309]GSE102316) of Ils1-ribotag and
Dhh-ribotag were retrieved from the GEO database respectively.
Acknowledgments