Foxm1 regulates neuronal progenitor fate during spinal cord regeneration

Mammals have limited tissue regeneration capabilities, particularly in the case of the central nervous system. Spinal cord injuries are often irreversible and lead to the loss of motor and sensory function below the site of the damage [1]. In contrast, amphibians such as Xenopus tadpoles can regenerate a fully functional tail, including their spinal cord, following amputation [2,3]. A hallmark of spinal cord regeneration is the re-activation of Sox2/3+ progenitor cells to promote regrowth of the spinal cord and the generation of new neurons [4,5]. In axolotls, this increase in proliferation is tightly regulated as progenitors switch from a neurogenic to a proliferative division via the planar polarity pathway (PCP) [6–8]. How the balance between self-renewal and differentiation is controlled during regeneration is not well understood. Here, we took an unbiased approach to identify regulators of the cell cycle expressed specifically in X.tropicalis spinal cord after tail amputation by RNAseq. This led to the identification of Foxm1 as a potential key transcription factor for spinal cord regeneration. Foxm1-/- X.tropicalis tadpoles develop normally but cannot regenerate their spinal cords. Using single cell RNAseq and immunolabelling, we show that foxm1+ cells in the regenerating spinal cord undergo a transient but dramatic change in the relative length of the different phases of the cell cycle, suggesting a change in their ability to differentiate. Indeed, we show that Foxm1 does not regulate the rate of progenitor proliferation but is required for neuronal differentiation leading to successful spinal cord regeneration.


Summary
Mammals have limited tissue regeneration capabilities, particularly in the case of the central nervous system. Spinal cord injuries are often irreversible and lead to the loss of motor and sensory function below the site of the damage [1]. In contrast, amphibians such as Xenopus tadpoles can regenerate a fully functional tail, including their spinal cord, following amputation [2,3]. A hallmark of spinal cord regeneration is the re-activation of Sox2/3+ progenitor cells to promote regrowth of the spinal cord and the generation of new neurons [4,5]. In axolotls, this increase in proliferation is tightly regulated as progenitors switch from a neurogenic to a proliferative division via the planar polarity pathway (PCP) [6][7][8]. How the balance between self-renewal and differentiation is controlled during regeneration is not well understood. Here, we took an unbiased approach to identify regulators of the cell cycle expressed specifically in X.tropicalis spinal cord after tail amputation by RNAseq. This led to the identification of Foxm1 as a potential key transcription factor for spinal cord regeneration. Foxm1-/-X.tropicalis tadpoles develop normally but cannot regenerate their spinal cords. Using single cell RNAseq and immunolabelling, we show that foxm1+ cells in the regenerating spinal cord undergo a transient but dramatic change in the relative length of the different phases of the cell cycle, suggesting a change in their ability to differentiate. Indeed, we show that Foxm1 does not regulate the rate of progenitor proliferation but is required for neuronal differentiation leading to successful spinal cord regeneration.

Foxm1 is specifically expressed in the regenerating spinal cord
We compared the transcriptome of isolated spinal cords at 1day post amputation (1dpa) and 3dpa to spinal cord from intact tails (0dpa, Figure 1A Figure S1E).
To identify the most enriched biological processes by gene ontology (GO), a nonbiased hierarchical cluster for all DE genes was performed ( Figure 1B). We observed three phases: first an increase in expression of genes involved in metabolic processes (cluster I), then a strong upregulation of genes associated with cell cycle regulation (cluster II and III) and finally, a downregulation of expression of genes involved in nervous system development (Cluster IV and V, Figure 1B).
Using Ingenuity Pathway Analysis (IPA), we identified potential upstream regulators that could explain changes in expression of downstream target genes, with Foxm1 showing the highest significance at 3dpa ( Figure 1C). Using published RNAseq of tail regeneration in X.tropicalis [9], we compared changes in expression of known Foxm1 target genes between 0 and 3dpa in whole tail (WT_d0d3), 0 and 1dpa (SC_d0d1) and 0 and 3dpa (SC_d0d3) in spinal cord. Foxm1 and its transcriptional targets are significantly upregulated only in the spinal cord at 3dpa, but not in the whole tail ( Figure 1D).
We wanted to confirm the expression of foxm1 during regeneration by in situ hybridisation (ISH) and RT-qPCR. ISH shows that foxm1 is not expressed in the spinal cord at 0 and 1dpa but is restricted to the regenerating spinal cord at 3dpa ( Figure 1E).
Pelzer et al. Role of Foxm1 during spinal cord regeneration 4 We then performed RT-qPCR for foxm1 over a period of 7 days, its expression peaks at 3dpa and decreases back to baseline levels at 7dpa ( Figure 1F).
We next wanted to identify the upstream signal(s) regulating its expression. As foxm1 expression starts at 3dpa, it is not a direct response to the injury. We tested if signalling pathways required for tail regeneration promote foxm1 expression at 3dpa. A sustained increase of reactive oxygen species (ROS) in the tail is required for its regeneration [10]. ROS levels were decreased following amputation using DPI, an inhibitor of the NADPH Oxidase (NOX). In NF50 tadpoles treated with DPI from 36hpa until 72hpa, foxm1 expression decreases by 69% (p=0.032) compared to DMSO controls (Figures S1F and S1G). ROS are upstream of different signalling pathways, including FGF [10,11].
Furthermore, Sonic hedgehog (Shh) signalling is also required for tail regeneration [12] and induces Foxm1 expression in the developing cerebellar granule neuron precursors [13]. Amputated tails treated with an FGF receptor kinase inhibitor (SU5402, Figure S1H), or a Shh signalling inhibitor (cyclopamine, Figure S1I) failed to regenerate (data not shown) but no changes in foxm1 expression was observed.

Foxm1 is required for spinal cord regeneration
To test the role of Foxm1 during regeneration, we designed a guide RNA (gRNA) targeted at bases 129-152 downstream of the ATG to knockdown and knockout foxm1 expression using CRISPR/Cas9 ( Figure S2A). The efficacy of the gRNA was assessed by Restriction Fragment Length Polymorphism analysis. Co-injection of the gRNA with cas9 mRNA did not lead to the destruction of the NcoI site but co-injection with 0.6ng and 1.5ng of Cas9 protein leads to NcoI-resistant PCR products in dose-dependent fashion ( Figure S2B). We then tested the level of indels by sequencing individual clones leading to 50 to 90% mutation on the foxm1 locus ( Figure S2C). Injection of CRISPR/Cas9 injection at 1-cell stage leads to a 52% reduction of foxm1 expression (p=0.0251; Figure S2D).
We then compared the ability of NF50 tadpoles injected with Cas9 alone (control) and Cas9 + foxm1 gRNA (foxm1 KD) to regenerate their spinal cord and tail. The ratio of regeneration was determined between 3 and 9dpa by dividing the length of the regenerating spinal cord by the length of the amputated spinal cord. No differences were observed at 3dpa, but from 5 to 9dpa the rate of regeneration was on average 40% lower (p=0.04) compared to controls (Figures 2A-C).
Could the impaired regeneration be caused by defective proliferation? To determine the rate of proliferation in the regenerating spinal cord, wt and foxm1-/-tadpoles were injected with EdU at 3dpa, followed by a 2-day chase (Figures 2D-F). As expected, we observed a higher proportion of EdU+ cells in the regenerate than in the non-regenerating spinal cord (~45% and ~20% respectively). However, no difference in EdU+ cells between wt and foxm1-/-was observed, showing that Foxm1 does not affect the overall length of cell cycle.

Characterisation of foxm1+ cells during regeneration
To understand the role of Foxm1 during spinal cord regeneration, we characterised this cell population at the molecular level using single cell RNA sequencing (scRNA-seq) ( Figures 3A and S3A).
As the cellular organisation of the Xenopus spinal cord is not well described, we used the 10XGenomics platform to sequence 2503 cells from uninjured spinal cord.
We then analysed the changes in the transcriptome during regeneration ( Figure   3D). A comparison of the clusters between 0dpa (yellow) and 3dpa (green) on a t-SNE representation shows that whilst some cell types cluster together (Schwann cells, DRGs), neural progenitors display a big shift in their transcriptome. Interestingly, foxm1 is expressed in a specific cluster of neural progenitors only at 3dpa ( Figure 3E). To characterise the foxm1+ cells, we identified DE genes between the foxm1+ cluster and the rest of the dataset ( Figures S3B and S3C). The list of the 20 top DE genes ranked by false discovery rate (FDR) shows that the majority are upregulated and many are linked to the cell cycle (ccna2, pcna, cdk2). We then identified the GO terms that were significantly over-represented with PANTHER and used Revigo to generate a plot representation ( Figure S3D). The majority of the GO terms identified are linked to the cell cycle. We therefore analysed changes in cell cycle dynamics between day0 and day3 ( Figure 3F).
Whilst clusters representing neurons are mainly in G1/G0 both at 0dpa and 3dpa, progenitor clusters have a higher proportion of cells in G2/M and S phases. Surprisingly, almost 50% of the cells in the foxm1+ cluster appears to be in S phase.
To confirm the changes in the cell cycle, we performed anti-PCNA staining on sections at 0, 3 and 5dpa ( Figures 3G, 3H and S3E). The percentage of PCNA+ cells in the spinal cord increases sharply at 3dpa when compared to 0dpa (from 18 to 68%) and remains high until 5dpa (88%; Figure 3G). To estimate the proportion of cells in different phases of the cell cycle, we used the fact that PCNA expression is punctate in S phase and diffuse in G1/G2 phase [14,15] and the chromatin is condensed in M phase. We observed a transient increase of the cells in S phase from 6.6% at 0dpa to 45.7% at 3dpa with a return to baseline by 5dpa (12%; Figure 3H)). These data confirmed our scRNAseq experiments and indicate that proliferative cells in the spinal cord have a long S phase in early regeneration.

Foxm1 regulates the fate of dividing progenitors during regeneration
Because of the changes in cell cycle and as Foxm1 promotes neuronal differentiation in early Xenopus development [16], we analysed the relative proportions of progenitors and neurons in the regenerating spinal cord when foxm1 expression is impaired. Expression analysis by RT-qPCR shows an 33% increase (p=0.02) in sox2 and a 18% (p=0.01) decrease in ntub expression in tadpoles with impaired foxm1 expression ( Figure 4A).
Interestingly, reducing ROS levels with DPI also impairs expression of ntub and ccnb3, a well-characterised Foxm1 transcriptional target. However, DPI treatment had no effect on the expression of ami, a gene expressed in endothelial cells ( Figure S4A). We next analysed Sox3 expression by immunofluorescence using anti-Sox3 antibodies in the non-regenerating and regenerating spinal cord of wt and foxm1-/-NF50 tadpoles at 5dpa. As expected, in the non-regenerating spinal cord Sox3 is expressed in the cells lining the ventricle. We also observed Sox3+ extensions into the mantle zone of the spinal cord ( Figure 4B, white arrowheads). In the wt regenerate, the spinal cord appears as an almost monolayer of cells around the central canal. Sox3 is expressed only in cells of the lateral spinal cord. These data reveal that the regenerating spinal cord conserved some cellular organisation. In contrast, in foxm1-/-tadpoles, we observed multiple cell layers and Sox3 expression is expressed more broadly ( Figure 4B). To assess the distribution of cells in wt and foxm1-/-spinal cords, the angle of all nuclei (DAPI, Figures S4B and S4C) and Sox3+ nuclei ( Figures S4D and S4E) relative to the centre of the central canal was calculated. Quantification of nuclei distribution either using cumulative number ( Figure S4B) or normalised distribution ( Figure S4C) show that whilst in the wt regenerating spinal cords, there are less cells laterally, this is not the case in foxm1-/-tadpoles. Furthermore, we observed a greater proportion of Sox3+ cells both dorsally and ventrally in the mutant compared to wt ( Figures S4D and S4E).
The proportion of Sox3+ cells in the non-regenerating spinal cord is about 50% in wt and foxm1-/-tadpoles. In contrast, the proportion of progenitors increased in the regenerating spinal cord from 62.1+/-3.5% in wt to 69.4+/-2.9% in foxm1-/-(p<0.001, Figure 4C). Next, we labelled cycling cells with EdU at 3dpa and determined their fate at 5dpa using immunofluorescence. We first quantified the proportion of dividing progenitors (Edu+Sox3+) over the total number of cells (S3+E+/DAPI) or over the total number of progenitors (S3+E+/S3+, Figure 4D) tadpoles following amputation. Both ratios are similar in wt and foxm1-/-tadpoles, indicating that knocking out foxm1 does not alter the rate of progenitor division. In contrast, the proportion of self-renewal (S3+E+/E+) increases significantly in the foxm1-/-compared to wt tadpoles (from 65% to 75%, p=0.006, Figure   4D). The stable rate of proliferation of progenitors combined with the increased rate in self renewal suggests that there is a shift in the fate of dividing progenitors from differentiation towards self-renewal.
To confirm these data, we analysed the expression of Myt1 as a neuronal marker in wt and foxm1-/-NF50 tadpoles at 5dpa ( Figure 4E). Knocking out foxm1 does not affect the percentage of Myt1+ cells in the non-regenerating spinal cord (~43% of spinal cord cells). However, in the regenerate we observed a sharp reduction of Myt1+ cells in the mutant compared to wt (30% vs 14%, p<0.01). Similar data were obtained when we analysed F0 injected with gRNA targeting the foxm1 locus compared to Cas9-injected controls ( Figure S4F). Thus, in foxm1 mutants, there is an increase of progenitors at the expense of differentiated neurons. Rather than promoting proliferation Foxm1 affects cell fate by promoting differentiation.

Discussion
Using a combination of bulk and single-cell RNAseq experiments on isolated spinal cord during tail regeneration in X.tropicalis, we identified a new population of cells present exclusively in the regenerating spinal cord. It is characterised by the expression of foxm1 and GO analysis reveals that genes involved in cell cycle and metabolism are overrepresented. Impairing foxm1 expression blocks spinal cord regeneration and alters the fate of the dividing progenitors with an increase in self-renewal division and a decrease in the production of neurons.
To characterise the Xenopus spinal cord at the molecular level, we undertook a single cell RNAseq approach. The number of cells and the depth of sequencing did not allow us to unambiguously determine the identity of progenitors and neurons subtypes but we identified the main cell types present in a vertebrate spinal cord, such as roof plate, dorsal and ventral progenitors, neurons and Schwann cells. Comparison of their transcriptomes at 0 and 3dpa reveals that differentiated cells did not respond much to injury. In contrast, the transcriptome progenitors displayed large changes as shown by their distinct locations on a t-SNE plot. Interestingly, at 3dpa, we isolate both the regenerate and the stump proximal to the amputation plan. However, no overlap of progenitors at 0dpa and 3dpa was observed, indicating that even cells positioned anteriorly from the amputation plan do respond to injury. Furthermore, foxm1 expressing cells are a subset of the 3dpa progenitors, suggesting that this cluster represents the cells present in the regenerate.
Foxm1 has been shown to have a role in neuronal differentiation during primary neurogenesis in Xenopus [16] and in the mouse telencephalon [17]. However, in both cases, the role of Foxm1 seems to be dependent on its ability to control the overall length of the cell cycle. This is not the case in the regenerating spinal cord, where we do not observe a difference in proliferation in the regenerating spinal cord in foxm1-/-compared to controls. This suggests a cell-cycle independent role for Foxm1 during this process. Cell cycle regulators such as Cdc25b and Ccnd1 have been shown to play a role during neuronal development independently of their primary function [18][19][20] suggesting that it might be a general principle of spinal cord regeneration [8,21]. It has been suggested that long S phase might be necessary for progenitors undergoing selfrenewal division to ensure genome integrity, especially in response to high level of ROS present in the nervous system [22][23][24]. The regenerating tail is an oxidative environment and we show here that foxm1 expression requires ROS. Furthermore, Foxm1 has been shown to ensure chromosomal stability and genome integrity in U2OS and aged fibroblasts [25,26], raising the intriguing possibility that Foxm1 might ensure that the expansion of the progenitor pool does not lead to genetic instability during regeneration.
The expansion of the neural stem cell pool is required for spinal cord regeneration in axolotl, zebrafish and Xenopus [4,8,27]. Here we show that regeneration also requires the precise control of neuronal differentiation. In mammals, ependymal cells also re-enter the cell cycle upon spinal cord injury (SCI), are able to self-renew but differentiate mainly into astrocytes [28]. Understanding how Sox2/3+ cells are re-activated and able to differentiate into neurons may open new opportunity to enhance spinal cord regeneration in species that have limited regenerative capabilities.

Xenopus tropicalis growth and tail amputation
X.tropicalis embryos were obtained and raised as previously described [29]. Tail amputation was performed at Nieuwkoop and Faber stages (NF)42-50 using a scalpel [30]. Tadpoles were anesthetized with 0.1% MS222 in 0.01X Marc's Modified Ringer (MMR) solution (10mM NaCl, 0.2mM KCl, 0.1mM MgSO4, 0.2mM CaCl2, 0.5mM HEPES, pH 7.4) for the procedure followed by recovery in 0.01X MMR. All animal procedures complied with the UK Animal (Scientific Procedures) Act 1986 and were conducted with UK Home Office approval.

RNA sequencing
Twenty spinal cords were isolated and immediately transferred into Trizol (Life Technologies) for each timepoint in triplicate. Following RNA extraction according to manufacturer's instructions, total RNA concentration was quantified using Qubit HS RNA assay kit (Invitrogen) on the Qubit Fluorometer 2.0 (Invitrogen). Integrity was tested with the Agilent RNA 6000 Pico kit on the Agilent Bioanalyser. Samples with an RNA Integrity Number (RIN) ≥7 were considered of acceptable quality. RNAseq was performed with Illumina NextSeq500 using unpaired-end sequencing at the GeneCore facility (EMBL).
Adapter sequences were trimmed using Trimmomatic v0.38. After quality control, the reads were converted to a FASTQ format and mapped on the v9.1 of the X.tropicalis transcriptome using bwa. The number of reads per transcript was determined using HTseq and the idxstats files were used to identify differentially expressed (DE) genes using the general linear model glmQLFit in DESeq2 with R. DE genes with a |log2(FC)|>1 and FDR<0.01 were considered significant. For the clustering, the k-mean was determined at 5 using the Elbow method and enriched gene ontology (GO) terms were identified using Fidea [31]. The full RNAseq dataset was then uploaded onto the Ingenuity Pathway Analysis (IPA, Qiagen) software. Upstream regulatory analysis was performed for DE genes with a |Log2FC|>1 and FDR<0.001 (992 genes for 0dpa/1dpa and 2720 genes for 0dpa/3dpa). This analysis predicts upstream molecules such as transcription factors that may be causing the observed change in gene expression [32]. RNA-seq data have been and viability (above 90%) were estimated on a Countess Cell Counter (Invitrogen). The cells were then loaded onto the 10X Genomics platform for processing.

Single cell RNAseq analysis
-Building custom genome for mapping: The X.tropicalis genome v.9.1 was downloaded from Xenbase.org. A number of genes had duplicated entries where the same gene ID was assigned to different gene names. The majority of these genes had an overlapping transcript, we therefore used the GFFReads merging option to merge these transcripts.
-Data pre-processing: The sequence files from the sequencer were processed using 10x -Gene filtering and Normalization: Genes with average UMI counts below 0.01 were filtered out. After this filtering 11,215 genes were left for downstream analyses. To take into account the effect of variable library size for each cell, raw counts were normalized using the deconvolution-based method [35]. Counts from many cells are pooled together to circumvent the issue of higher number of zeros that are common in scRNA-seq data.
This pool-based size-factors are then deconvoluted to find the size factor of each cell.
These normalized data are then log-transformed with a pseudo-count of 1.
-Visualization and Clustering: The first step for visualization and clustering is to identify the Highly Variable Genes (HVGs). To do this, we first decomposed the variance of each gene expression values into technical and biological components and identified the genes for which biological components were significantly greater than zero. These genes are called HVG. HVG genes were then used to reduce the dimensions of the dataset using PCA. The dimensions of dataset were further reduced to 2D using t -SNE, where 1 to 14 components of the PCA were given as input. The cells were grouped into their putative clusters using the dynamic tree cut method. Dynamic tree cut method identified the branch cutting point of a dendgoram dynamically and combined the advantage of both hierarchical and K-medoid clustering approach. This method identified 15 clusters in the population.
-Identification of marker genes: To identify the marker gene for a cluster, we compared a given cluster with all other clusters. We selected the marker genes for the cluster using manual curation of DE genes with an FDR<0.01. For PANTHER analysis, all DE genes between the foxm1+ cluster were uploaded and tested against the pseudo-bulk of all the genes expressed in the scRNAseq experiment. scRNA-seq data have been deposited in the ArrayExpress database at EMBL-EBI (www.ebi.ac.uk/arrayexpress) under accession number E-MTAB-8839.
For chemical inhibitor treatments, tadpole tails were amputated at NF50 and left to recover for 36h in 0.01XMMR. Then, ROS signalling was inhibited with 4µM diphenyleneiodonium (DPI, Merck), FGF signalling with 20µM SU5402 (Calbiochem) and Shh signalling with 2.5µM of Cyclopamine (Merk). At 3dpa, the regenerates were collected, total RNA extracted and then processed for RT-qPCR.

Genotyping
Genomic DNA was extracted by incubating the embryo or a section of the tail removed by amputation for 3h at 55ºC in 10mM Tris pH8, 1mM EDTA, 80mM KCL, 0.3% NP40, 0.3% Triton X100 and 0.2mg.mL -1 of proteinase K. The samples were subsequently processed for PCR amplication using the following primers: fwd 5'-GTATGTTGCAGAGCAGGGCAT, rev 5'-GTATGTTGCAGAGCAGGGCAT. The PCR product was subsequently digested with NcoI for Restriction Fragment Length Polymorphism (RFLP) analysis.

Isolation of RNA and qPCR
RNA was isolated from total embryos, tail regenerates and isolated spinal cord using Trizol

In situ hybridisation, EdU labelling and immunofluorescence
Whole mount in situ hybridisation on embryos and tails was performed as previously

Statistics
When 2 or more conditions were compared, the normality of the distribution of the data was tested using a Shapiro-Wilk test. For two conditions, a two-tail unpaired t-test was used. If three or more conditions were compared, a one-way ANOVA followed by a Tukey Post-hoc test was used. For the rose plots and the cumulative distribution( Figure S4), a Kolmogorov-Smirnov test was used. All the statistical tests were done using Prism 8.