Stemness Activity Underlying Whole Brain Regeneration in a Basal Chordate

Understanding how neurons regenerate following injury remains a central challenge in regenerative medicine. Adult mammals have a very limited ability to regenerate new neurons in the central nervous system (CNS). In contrast, the basal chordate Polycarpa mytiligera can regenerate its entire CNS within seven days of complete removal. Transcriptome sequencing, cellular labeling, and proliferation in vivo essays revealed that CNS regeneration is mediated by a newly formed neural progeny and the activation of neurodevelopmental pathways that are associated with enhanced stem-cell activity. Analyzing the expression of 239 activated pathways enabled a quantitative understanding of gene-set enrichment patterns at key regeneration stages. The molecular and cellular mechanisms controlling the regenerative ability that this study reveals can be used to develop innovative approaches to enhancing neurogenesis in closely-related chordate species, including humans.


Introduction
One of the most challenging questions in regenerative medicine is that of how to promote neural regeneration within the adult human brain. In mammals, regeneration is mediated by tissue-specific stem cells capable of self-renewal and multilineage differentiation [1]. While the involvement of neural stem cells (NSC) in neurogenesis of the mammalian brain has been identified [2], the rarity of long term NSCs and the difficulty of monitoring NSCs in vivo, limits our ability to study the mechanisms that direct their activation and differentiation [3,4]. Non-mammalian animals, such as teleost fish, amphibians, and invertebrate species such as tunicates, however, possess expansive regenerative abilities that enable them to replace injured nerve cells, and are used to study the cellular and molecular underpinnings of neural regeneration [5][6][7].
Comparing regeneration mechanisms among a wide array of model systems with varied regenerative capabilities suggests that differences in CNS regeneration ability are due to variations in several key features, including changes in wound-healing responses in neurogenic cell populations and changes in the genetic circuitry of the stem cells that maintains access to embryonic transcriptional programs [8]. Comparing the mouse and zebrafish brain for example, revealed that unlike the mouse, the zebrafish adult CNS contains numerous neurogenic niches with active neuronal progenitors that maintain a life-long ability to produce new neurons in response to injury [9]. In the salamander, CNS regeneration involved the expression of embryonic morphogens to stimulate cell growth, re-patterning, and diversification [8].
In the current study, we developed a new model system to study whole CNS regeneration. By utilizing the solitary tunicate Polycarpa mytiligera's natural ability to regenerate all its body organs and tissues, including a simple brain and CNS, from a small body fragment [6], we have gained insight into neuroregeneration in a species at the base of the chordate evolutionary tree. Tunicates are basal chordates, a sister group of vertebrates, that share structures and cell types considered to be homologous to those in vertebrates [10], and are used as model systems for chordate development and regenerative studies [11,12]. A population of circulatory putative stem cells was suggested to mediate regeneration in this group [13][14][15][16][17][18][19]. Recent transcriptome analyses suggest that these stem cells are tissue specific [10,11]. Following injury, undifferentiated cells accumulate at the regenerating area, forming a proliferation zone capable of producing regenerative tissue [6,16,[20][21][22]. A recent study suggests that following CNS ablation, candidate stem cells differentiate into specialized cell types, including neurons [23].
While considerable progress has been made in characterization of the cellular events that lead to CNS regeneration in tunicates, only a few studies have described the expression of morphogens and their transduction pathways that stimulate cell growth, re-patterning, and diversification to regenerate the CNS [11,12,24]. To overcome this gap, we integrated transcriptome sequencing of major CNS regeneration stages with cell labeling and cell proliferation in vivo essays, to characterize the stepwise events and genetic changes that lead to whole CNS regeneration in P. mytiligera.

Experimental Design
Our main research aim was to provide a detailed description of CNS regeneration in a new chordate model system. To this end we characterized the morphological, functional, and molecular processes underlying P. mytiligera CNS regeneration.
Morphological analysis. We first described the control CNS in juvenile and adult animals in order to better understand this system's anatomy and to determine the end point of the regeneration process. Next, we monitored the reconstruction of the CNS using various methods, including histology, immunolabeling, and in situ hybridization at different time points following CNS removal ( Figure 1). Our findings enabled us to divide the regeneration process into three key stages that reflect similar morphological criteria in both age groups (Figure 2A). (H) Transmission electron microscopy images of P. mytiligera neural complex showing the brain structures; medulla (ME) and cortex (CO). Synapse in the medulla is marked with yellow arrows. Atrial siphon (AS), digestive system (DS) endostyle (EN) and oral siphon (OS).  Enlargements of the square area appear in the following image in the panel. The CNS is absent following amputation. By 2 dpa we see the beginning of the closing of the open wound followed by complete regeneration of the CNS by 21 dpa. Anterior siphon (AS) body wall (BW) cerebral ganglion (CG), cortex (CO), medulla (ME), nerve fiber (N) neural gland (NG), and oral siphon (OS). (C) Whole mount of regenerating animals showing troponin expression in regenerating muscle fibers at 2 and 5 dpa. (D) Transmission electron microscopy images of regenerating muscle fibers at 6 dpa. (E) Whole mount of regenerating animals showing tubulin expression in regenerating nerves fibers at 2 and 5 dpa. (F) Transmission electron microscopy images showing the brain structures and cellular components in control (nonamputated) and 6 dpa animals. Synapse in the medulla is marked with yellow arrows. Immune cells are marked with red arrows. (G) Gene enrichment plot of regeneration-associated and CNS related gene sets during CNS regeneration. Light shaded regions indicate the 50% and 99% confidence intervals under a hypergeometric model. Functional analysis. Cell proliferation underlies the formation of new tissues in diverse regenerative organisms, including P. mytiligera. Our goals were therefore to better characterize the identity of post-amputation dividing cells and study their involvement in nerve regeneration. First, we sought to describe the location of the proliferation events and determine whether CNS removal resulted in local or systemic cell proliferation. Next, we sought to monitor the identity of the proliferating cells using the tissue-specific markers alpha-tubulin (nerves and cilia) [25][26][27] and troponin (muscle) [28]. To achieve these goals, we applied the EdU proliferation assays to regenerating juveniles and monitored the location of cells that entered the S-phase following CNS removal. We then analyzed the expression of tubulin and troponin in these dividing cells to determine the regenerative stage at which they acquire their identity (Figures 2 and 3).  Molecular analysis. We performed de novo transcriptome analysis using CNS samples of adult animals from the three chosen regeneration stages, in order to analyze the molecular landscape of neural regeneration. As adult stem cells have been shown to mediate whole body regeneration in other tunicates, we hypothesized that CNS regeneration would involve the activation, proliferation, and specialization of a similar population in our highly regenerative model. Focusing on regeneration and stem cell markers, we compared the expression patterns of each regeneration stage with the control and with all other time points in order to identify the enriched pathways and key transcription factors activated or inhibited during CNS regeneration (Figures 2, 4, S3 and S4). This comprehensive analysis provided a valuable insight into the factors and pathways that conduce wound healing and neurogenesis, as revealed by our new chordate model system.

Animal Collection and Culturing
During March-April 2018 P. mytiligera adult individuals were collected from the Gulf of Aqaba (Eilat), Red Sea (Israel). The animals were maintained at the Eilat Inter-University Institute (IUI) in aquaria with running seawater for four days of acclimation prior to the onset of the experiments. P. mytiligera juveniles were obtained from a breeding culture established in the IUI, as previously described [29].

Regeneration Experiment
For the morphological description of the regeneration process, and selection of representative time points, adults (n = 24) of similar size range (4 ± 1 cm, from the tip of the oral siphon to the base of the animal) were anesthetized [30] and dissected. The CNS was removed using a scalpel. Both treated and non-amputated (control) animals (n = 4) were maintained in running seawater for the duration of the experiment. Response to touch was used to determine survival. Juveniles of the same age (3 months old, n = 24) were dissected, and their CNS was removed. Due to their small size, the oral siphon and the CNS were amputated together to ensure removal of the entire CNS.

Histology and Immunohistochemistry
Regenerating adult animals (n = 20) were fixed at different time points along the regeneration process (amputation day, 2-, 7-, 14-, and 21 post-amputation). Following relaxation with menthol, the animals were fixed in 4% formaldehyde in seawater, partly separated from the tunic, and dehydrated with ethanol at increasing concentrations prior to paraffin (Paraplast Plus, Lecia) embedding. Sections (7 µm) were mounted on glass slides, deparaffinized, and stained with standard hematoxylin and eosin solutions.
For immunolabeling, chosen sections were rehydrated in xylene and graded series of ethanol to phosphate-buffered saline (PBS). Juveniles (n = 10) were fixed in 4% paraformaldehyde in phosphate-buffered saline (PBS) overnight at 4 • C. Fixed juveniles were then washed three times (5 min each) in PBS followed by a 10 min wash in TritonX-100 at room temperature (RT). Samples were blocked in PBT + 3% BSA for 3 h at RT. Antibodies were added directly to the blocking solution overnight at 4 • C. To visualize nerves and cilia, we used a mouse monoclonal anti-acetylated tubulin antibody (Sigma T7451) diluted 1:1000 in blocking solution (36,3,7). Samples were then washed twice for 10 min each in PBT. Secondary antibodies (Alexa Fluor goat anti-mouse IgG 488, ThermoFisher, Waltham, Massachusetts, USA, A-11001) were added at 1:500 dilution to the blocking solution for 3 h at RT. Samples were then washed in PBS three times for 10 min each. Samples were stained with DAPI (Sigma-Aldrich, MO, USA, 1 µg/mL in PBSTx) and mounted. Images were taken using a fluorescent microscope.
After washing in buffer and post-fixation in 1% OsO4 in 0.2 M cacodylate buffer, the specimens were dehydrated and embedded in epoxy resin (Sigma-Aldrich, MO, USA). Ultra-thin sections (80 nm thick) were stained with uranyl acetate and lead citrate to provide contrast. Photomicrographs were taken with a FEI Tecnai G12 electron microscope operating at 100 kV. Images were captured with a Veleta (Olympus Soft Imaging System) digital camera.

An RNA-Seq Catalog for Central Nervous System Regeneration
To characterize the molecular events that take place during P. mytiligera regeneration, and to establish a reference map of CNS regeneration amenable to inter-species comparison, we profiled 12 samples of CNS tissue from adult animals and surveyed gene expression changes. CNS samples (n = 3) from untreated adult animals were used to reflect the start and end points of regeneration. Based on these morphological features, pools of regenerating CNS tissues were collected at five time-points: 2 days (n = 2), 7 days (n = 2), 14 days (n = 3), and 21 days (n = 2) post-amputation, to provide an informative perspective on the molecular events leading to CNS regeneration.

RNA-Seq, Read Mapping, and Transcriptome Assembly
Samples were prepared at the IUI during September 2017. Fifteen P. mytiligera adult individuals (6 ± 1 cm, from the tip of the oral siphon to the base of the animal) were used for neural complex regeneration experiments. Neural complex tissue samples were prepared for RNA extraction (RNAeasy mini kit from Quiagen) at different time points along the regeneration process. Total RNA was extracted following the manufacturer's protocol. cDNA libraries were then prepared from high-quality samples (RNA integrity number (RIN) > 8) at the Weizmann Institute of Science, Life Sciences Core Facilities, Israel. Barcoded library samples were sequenced on an Illumina NextSeq 500 (2 × 150 bp) at Stanford University, CA, USA.
Only transcripts with p_good >0 score were used for further analysis. TransDecoder [35] (v5.5.0) was used for identifying candidate coding regions within the transcript sequences. Homologous sequences were identified and annotated, based on the candidate coding region transcripts, using BLASTP (NCBI), cut-off 1e-10, to align with (1) the Swis-sProt protein database; and (2) solely on the Mus musculus curated SwissProt database (downloaded on 21 February 2020).
Sets of genes were tested for enrichment of Gene Ontology (GO Biological Process, Molecular Function and Cellular Component) terms. For a set of genes with significant up and down effects, found using the Wald test, an over-representation analysis (ORA) was performed using the enricher function of the clusterProfiler [40] package (v3.16.1) in R with default parameters. A gene set enrichment analysis (GSEA) was also performed on the entire assembly using a scoring based on the log fold change of each gene on its respective time point using the GSEA function of clusterProfiler [40] package (v3.16.1) in R with default parameters. In both cases Significant GO terms were identified with an FDR < 0.05.
The gene ontology terms were obtained from a manually created database based on the SwissProt curated Mus musculus GO annotations, using the makeOrgPackage function of AnnotationForge [41]. Pipeline summary is presented in Figure S1.

Identification of Regeneration-Associated Differentially Expressed Genes
The generated dataset was analyzed using two methods: Identification of differentially expressed genes between each time point and the control homeostatic CNS sample. For each time point, transcript levels were compared between each pool of regenerating fragments and the control sample (Table S1).
Identification of differentially expressed genes between all possible combinations of all time points. The identification of chronologically differentially expressed genes and the formation of binary tables have been described previously [11]. Briefly, we used DEseq2 [39] (REF; FDR < 0.05) to identify the differentially expressed genes between all possible combinations of contiguous and individual time points, resulting in a hierarchy of up and down regulated time points for each gene for each age group. For each gene in each experimental group its specific time point hierarchy was then used to score all the possible binary patterns, with each pattern's score being the number of shared up-anddown regulated time points between it and the hierarchy, while subtracting the number of up-and-down regulated time points for which the pattern and hierarchy disagreed. The pattern with the highest score was used for that gene and regeneration stage pairing. Based on these analyses, a binary gene-time expression matrix for every expressed gene recorded along the time points was produced, with 1 indicating dynamically "high" expression and 0 indicating dynamically "low" expression (Table S2). Using the binary matrix, we identified pathways and Go terms associated with each time point (Figures 2G,  4A and S3B, Table S3).

Gene Enrichment Plots
Gene enrichment plots have been previously described in [11]. Briefly, at each time point, the proportion of genes in a gene set that are active (indicated by a 1 in the gene-time expression binary matrix defined above) is calculated. This gives a value between 0% (no genes in common) and 100% (all genes in the gene set are active at that time). A baseline expectation of the proportion of overlapping genes is calculated using a hypergeometric model that gives the likelihood that the same number of genes as in the selected gene set would be randomly selected from the matrix. In addition, the 68% confidence interval (1 standard deviation) of the proportion of shared genes ('enrichment') from the hypergeometric model is calculated, plotted, and presented as a shaded region in the plot. The baseline is then subtracted from the values calculated, with the confidence interval also subtracted, to present the expected range of values and the extent to which the actual enrichment result differs from a null result. If the baseline expectation is greater than the actual enrichment (i.e., the subtracted value will be negative) a value of 0% is used (as a negative percentage is considered meaningless).

Use of Specific Pathway Gene Sets
For the gene enrichment plots the following gene sets were used: For the pathway and go terms identified by GeneAnalytics within our time data, we focused on the pathways of the same names from PathCards, an integrated database of human biological pathways and their annotations. PathCards clustered Human pathways into SuperPaths based on gene content similarity (https://pathcards.genecards.org/, accessed on 21 July 2021). From each of these gene sets we removed any gene that did not have a putative homolog to a known Polycarpa gene. Consequently, the percentages within the enrichment plots refer to these curated Polycarpa specific gene sets. If a gene name appeared more than once in the Polycarpa gene model annotation, all matching Polycarpa gene ids were included in the gene set.

Gene Cloning and Transformation
Gene-specific primers were designed from the transcriptome sequence and synthesized by Integrated DNA Technologies (IDT). Their oligonucleotide sequences were as follows: Troponin C 5 GGATTTGACGGAAGAGCAGA3 and 3 TCATTGCCACGAACTCTTCA5 . Tubulin alpha-1A chain 5 CTGCAGACGAAACCTTGACA3 and 3 TAAACCGTATCACCG TGCAA5 .
Genes were amplified from P. mytiligera cDNAs using gene-specific primers and cloned into pGEM-t vectors using the manufacturer's protocol (Promega; CAT #A1360). Vectors were transformed into E. coli Top10 by the heat-shock method. Briefly, 100 µL of bacteria were mixed with 5 µL of cloned vector, incubated on ice for 30 min, and then

Fluorescent In Situ Hybridization
Fluorescence in situ hybridization (FISH) was based on the published protocol [42,43] with some modifications (Appendix A). Essentially, fixed juvenile animals were separated from their tunic and opened along the endostyle. Following bleaching and proteinase K treatments, samples were incubated with DIG/ DNP-labeled riboprobes overnight in hybridization solution, at 60 • C. After hybridization, samples were washed twice in: prehyb solution, 1:1 pre-hyb-2X SSC, 2X SSC, 0.2X SSC, PBSTx. For the blocking step prior to antibody incubation, 0.5% Roche Western Blocking reagent and 5% inactivated horse serum in 1xPBSTx were used. Samples were then incubated with antibodies (anti-DIG-POD, Fab fragments) overnight at 4 • C. Antibody solution was washed with 1xPBSTx. Rhodamine or FITC dyes were used in tyramide development. For peroxidase inactivation, samples were washed with 1% sodium azide solution for 90 min at room temperature. Finally, samples were counterstained with DAPI (Sigma-Aldrich, MO, USA, 1 µg/mL in PBSTx) for 1 h, mounted and photographed using a Zeiss LSM 880 scanning laser confocal microscope.

In Vivo Cell Labeling Experiments
Cell proliferation was detected by incorporating 5-ethynyl-2-deoxyuridine (EdU) into replicating DNA. P. mytiligera juveniles (3 months old) were divided into two groups, dissected, and separated from their neural complex and oral siphon using a scalpel (Extended data Figure 3a). One group (n = 3) was exposed to a 16 h EdU pulse following amputation and was fixed immediately afterwards. The second group (n = 9) was exposed to a 16 h EdU pulse 32 h following amputation and left in seawater to regenerate. Fixation was done at 2-, 5-, and 7 days post-amputation (n = 3 per each time point) (Figure 3).
For the pulse experiments animals were incubated with 10 µmol/L EdU (Invitrogen, Carlsbad, CA) in 5 mL of MFSW for 16 h in Petri dishes. Following completion of the labeling, animals were fixed for 12 h in 4% FA, rinsed three times in 1 × phosphate-buffered saline (PBS), and processed for EdU detection using Alexa Fluor azide 488 at room temperature, according to the instructions of the Click-iT EdU Alexa Fluor High Throughput Imaging Assay Kit (Invitrogen). Samples were stained with DAPI (ThermoFisher 33342) (1 µg/mL in PBS) and mounted in VECTASHIELD (Vector Laboratories RK-93952-28) using coverslips.

Statistical Information
All results are expressed as mean ± SE. Statistical analyses were performed using R Studio. Statistical analysis was performed using two-sided Mann-Whitney Wilcoxon test throughout the study. Difference was considered significant as follows: * p < 0.05.

Results and Discussion
To obtain a precise description of the morphological and cellular events underlying P. mytiligera CNS regeneration, a detailed description of the CNS anatomy of juvenile and adult animals was acquired using naive animals ( Figure 1). Our results revealed a similar morphology in both life stages. The CNS is composed of a single ganglion (brain) [44] and connected by a network of nerves to the different body parts. P. mytiligera's brain consists of an outer cortex and a central neuropil ( Figure 1D,H). Nerve cell bodies are present in the cortex, surrounding a medulla where synaptic connections can be seen ( Figure 1H). Four main anterior nerves, two posterior nerves, and one ventral nerve exit from the brain ( Figure 1F,G). The anterior nerves innervate the incurrent siphon and the anterior body wall, while the posterior nerves innervate the excurrent siphon and the remaining body wall. The ventral visceral nerve extends over the branchial basket toward the digestive tract ( Figure 1G). The brain is located adjacent to the neural gland ( Figure 1A,B,E), a nonnervous structure that shares with the brain a common embryonic origin [45]. Together, the brain and the dorsal neural gland comprise the neural complex. P. mytiligera's CNS presents a simple brain composed of tubulin-positive nerve cells that interact via synaptic connections, representing a basic model comparable to that in other chordates. Based on these findings, we were able to monitor the CNS regeneration process, differentiate between distinct developmental stages, and determine its endpoint (Figure 2A). Both juvenile and adult animals were able to survive and regenerate following CNS removal. However, the time needed to complete a successful regeneration strikingly differed between the two age groups: young animals regenerated their entire CNS within 7 days, while older animals needed 3 weeks to complete the same process (Figure 2A). To study the cellular and transcriptomic signatures associated with CNS regeneration, we compared CNS samples taken from control animals (to establish a basal state) with CNS collected at diverse time points along the regeneration process. Based on these morphological and transcriptomic results, CNS regeneration was divided into three milestones: early, middle-, and late-stage regeneration (Figure 2A).
To better characterize the cells that had contributed to the regenerated tissues we labeled endogenous dividing cells through the systemic administration of EdU, which incorporated into the DNA during the S-phase of cell division and stained the nucleus of the cells that underwent division within a few hours post-treatment (Figure 3). Fluorescent in situ hybridization of tissue-specific markers was then used to track these newly dividing cells and identify the ones that had differentiated into neurons ( Figures 2E and 3B) [46]. Our results indicate a possible contribution of the undifferentiated progenitor cells that enter the S-phase cellular state during the early regeneration stages, to the formation of the regenerated CNS. Following amputation of the CNS, only the isolated tubulin-positive nerve fibers in the surrounding tissue remained (Figure 2B). At the initial stage postamputation epidermal tissue started to regenerate and close the open wound ( Figure 2B). This tissue was composed of regenerating troponin-positive muscle fibers ( Figure 2C,D), while nerve fibers had not yet began to re-innervate the tissue ( Figure 2E). During the first 16 h post-amputation (hpa) no significant difference in the number of EdU positive cells was found between non-regenerated and regenerated areas. At 48 hpa a higher number of proliferating cells were located specifically in the regenerated area, indicating a local increase in proliferation rate. At this time point, these dividing cells did not express CNS specific markers (Figures 2E and 3B).
During the mid-regeneration stage, the wound became fully healed, and the CNS regenerated, with newly formed nerve fibers exiting from the brain and continuing along the regenerating tissue to the siphons ( Figure 2E). This stage was also characterized by the highest number of EdU-positive cells. Some of these cells were found embedded in the regenerated brain and nerve fibers and were tubulin-positive, indicating their neuronal identity ( Figure 3B). This finding indicated a recent differentiation event, as these neural progenies had acquired their identity as nerve cells at this time point, following their proliferation during the early regeneration stage.
At the late regeneration stage, the de novo brain morphological structures and cellular distribution pattern resembled those of the control ( Figure 2B). TEM analyses verified that the regenerated brain contained a cortex and a medulla with visible synapses. However, the brain morphology was still not completely identical to the control brain at this point, as the cortex layer was not yet well organized into distinct layers ( Figure 2F). In addition, we observed several immune cells embedded in the cortex layer ( Figure 2F), implying a possible contribution of the immune system to the regeneration process [47].
To characterize the molecular signature of neural regeneration we performed the first transcriptome assembly of a solitary chordate CNS regeneration ( Figure S1).
The generated dataset was then analyzed using over-representation (ORA) and functional enrichment (GSEA), to identify the key pathways that had been enriched during the CNS regeneration in comparison to the control CNS (non-amputated brain) ( Figure S1). Moreover, by combining the DEseq2 [39] and bioinformatic tools that we had developed (all vs. all analysis [11]), we were able to construct a detailed gene expression atlas of the differentially expressed gene sets for each stage of the regeneration process (Table S2). Our analyses revealed the dynamic expression of 239 pathways, demonstrating significant changes in gene expression along the regeneration timeline.
The early regeneration stage was characterized by activation of cell proliferation and chromatin organization pathways ( Figure 2G), consistent with the results of our EdU assay. The main upregulated GO terms comprised cell signaling, nerve proliferation, and injury response pathways, reflecting the animal's fast response to the amputation and damaged tissues ( Figure S2). Interestingly, this stage was enriched in several stemcell-related pathways ( Figure 4A), including the p53 pathway, a stress response pathway required for neurite outgrowth as well as for axonal proliferation and regeneration [48]. In the mouse brain, p53 is involved in maintaining regenerative ability by regulating the proliferation of stem and progenitor cells [49]. p53 has been shown to play several roles, including apoptosis, cell cycle arrest, and proliferation, depending on the different stages of the regeneration process [50]. In P. mytiligera, activation of this pathway occurs at the early regeneration stage and is inhibited at the later stages. These findings indicate a negative regulation of neural stem cell self-renewal and an investment in specialization processes at the early time point, followed by stem cell maintenance as the CNS undergoes reconstruction [51,52]. p53 transcriptional activities are primarily regulated through posttranslational modifications, including sumoylation [53]. The sumoylation pathway was also enriched during the early stage of CNS regeneration ( Figure 4A). This pathway has a role in neural crest development [54] and stem-cell proliferation [55], and was shown to increase in mice following nerve lesion [56]. Another post-transcriptional regulatory pathway that was highly enriched following P. mytiligera's CNS ablation was the Piwi-piRNA pathway ( Figure 4A). This pathway is involved in stem cell maintenance in diverse organisms [57], and was suggested to have an inhibitory role in neuron or Schwann cell responses during peripheral nerve injury in both nematodes and rodents [58].
As regeneration progressed, dynamic changes in the enriched pathways could be observed. Upregulation processes associated with differentiation and tissue formation were expressed, as well as the enrichment of extracellular matrix organization processes ( Figure S2). Enrichment of Vegf signaling indicates a revascularization process, an essential process during regeneration and consistent with the morphological description ( Figure 4A). At this stage the pathways associated with stem-cell activity that had been detected at 2 dpa were replaced by a different set of stem-cell associated pathways, including Notch, Hedgehog, and Wnt ( Figure 4A), indicating a possible shift from stem cell proliferation and differentiation to stem, cell self-renewal and maintenance [50,[59][60][61][62]. Closer observation of the gene sets composing these pathways revealed potential key regulators of CNS regeneration, including wnt7b, wnt5a, β-catenin, low-density lipoprotein receptor-related protein 5 (lrp5), and lrp6, all of which demonstrated a dynamic expression pattern along the regeneration process ( Figures 4B,C, S3 and S4). Components of the Notch signaling pathway included upregulation of the receptor protein gene notch-1 and its ligand, delta, down-regulation of the signal mediator rbpj, and dynamic expression of a downstream transcriptional target hey-2 ( Figures 4B,C, S3 and S4).  The middle and late regeneration stages demonstrated an enrichment of the signaling pathways associated with tumor necrosis receptor factor 1 (TNFR1), fibroblast growth factor receptor (FGFR), and nerve growth factor (NGF) ( Figure S3B). These pathways play key roles in controlling the normal inflammatory and wound response to muscle and neuronal damage, necessary to achieve regeneration [47,63].
The late stage de novo regenerated brain revealed an upregulation of GO terms related to extracellular matrix secretion and axis specification ( Figure S2B), implying a continuous process of tissue assembly and specialization. Additional support for such a shift came from the upregulation of voltage-gated ion channel genes, including kcnh8 and cacnb2 ( Figure  S3A), and the enrichment of gene sets associated with the Nanos pathway ( Figure S4B) and neurotransmitter release cycle ( Figure 2G), reflecting morphogenesis processes and the establishment of functional synapses [64].
Next, we sought to identify key regulatory factors that are known to be involved in CNS regeneration in other model organism, that were differentially expressed during regeneration in P. mytiligera, and to determine precisely when their expression was enhanced or inhibited along the regeneration process in relation to the control (Figures 4B,C, S3 and S4). Our results revealed a dynamic expression of transcription factors related to neuronal differentiation and neural cell fate specification ( Figures 4B,C, S3 and S4). Among this gene list, sox10, tbr1, and pax-6 were differentially expressed following CNS removal. Fluctuation in expression patterns of these markers has been previously shown to mediate CNS regeneration in other chordates [65,66]. These findings indicate a possible conserved role of these factors during nerve regeneration in P. mytiligera, as both sox10 and tbr1 were upregulated throughout the entire process, while pax-6 was consistently down-regulated.
The nervous system development genes, foxa2 and six1, were also upregulated during P. mytiligera CNS regeneration ( Figures 4B,C, S3 and S4). A similar trend in expression patterns of these factors underlies regeneration in mice, where foxa2 is required for the specification and differentiation steps of neurons at progressively higher doses [67] while Six1 expression is required for the proliferation of neuroblast progenitors [68]. Interestingly, the neurogenic transcription factors c-myb and sox2 were upregulated during the late stage of CNS regeneration. The universal neural stem cell marker Sox2 is involved in the proliferation and maintenance of neural stem cells as well as in neurogenesis [69]. The hematopoietic transcription factor and proto-oncogene c-myb regulates neural progenitor proliferation and has an important regulatory role in the neurogenic niche of the adult mouse brain. The expression of these neural stemness markers following CNS removal strongly suggests the involvement of stem cells during regeneration in P. mytiligera.
Our study presents gene sets and pathways known to be associated with neural regeneration and stem cells and describes the exact time point in which their expression was upregulated or down-regulated (Table S3). Our results indicate that CNS regeneration in P. mytiligera shares basic functional stages with other regenerative animals. These include an increase in cell proliferation and reorganization process at early stages of regeneration, followed by the activation of specialization and morphogenesis processes at late regenerative stages ( Figure 5). While we chose to focus on conserved pathways, known to mediate regenerative response in other organisms, our analyses revealed novel relations between the expression dynamics of some of these pathways that were unique to our model. These include the high expression of P53 and piRNA pathways preceding the activation of Notch, Wnt, and Nanos pathways at early stages of the CNS regeneration. As a member of an evolutionary clade considered to be the closest living invertebrate relative of vertebrates, our CNS regeneration atlas may well constitute a fundamental resource in understanding how evolution has shaped CNS regeneration processes across phylogeny, and in obtaining mechanistic insight into this complex process. Utilizing a variety of model organisms contributes greatly to our understanding of regenerative biology, and establishing new models bears the potential to uncover new insight into this complex process. However, working with non-model organisms has its limitations and challenges [70,71]. One of the main challenges of the current study was the low number of available specimens as well as availability of experimental tools and protocols. To overcome these challenges, we established new protocols (Appendix A) and divided our pool of accessible animals into several selected experiments. Juvenile animals were used for FISH and EdU experiments, where their small size and transparency represent a significant advantage, while adult animals were used for RNA extraction and gene expression analysis as their large body size enables accurate removal of the CNS. Another limitation of working with a non-model is the need to establish a broad and somewhat basic foundation for future research. Our regeneration atlas has revealed the main pathways and genes that are activated or inhibited throughout the entire CNS regeneration Utilizing a variety of model organisms contributes greatly to our understanding of regenerative biology, and establishing new models bears the potential to uncover new insight into this complex process. However, working with non-model organisms has its limitations and challenges [70,71]. One of the main challenges of the current study was the low number of available specimens as well as availability of experimental tools and protocols. To overcome these challenges, we established new protocols (Appendix A) and divided our pool of accessible animals into several selected experiments. Juvenile animals were used for FISH and EdU experiments, where their small size and transparency represent a significant advantage, while adult animals were used for RNA extraction and gene expression analysis as their large body size enables accurate removal of the CNS. Another limitation of working with a non-model is the need to establish a broad and somewhat basic foundation for future research. Our regeneration atlas has revealed the main pathways and genes that are activated or inhibited throughout the entire CNS regeneration process in an invertebrate chordate species, and suggests the involvement of candidate stem cells in this process. These findings will pave the way for future studies aimed at functionally isolating neural stem cells in adult tissues, and testing the involvement of key genes in directing their differentiation to form de novo whole CNS.

Conclusions
As the closest living relative of vertebrates, the tunicates serve a critical role in understanding developmental processes that are comparable to those of vertebrates. Our study has revealed a detailed road map of complete CNS regeneration, in which the activation and inhibition of candidate neurogenesis factors are documented. CNS removal initiates the proliferation of progenitor cells that differentiate into neurons and trigger the expression of known stem-cell transcriptional regulators, including sox2, hes1, and c-myb. As recruiting endogenous and exogenous NSC for CNS regeneration has become one of the main goals of regenerative medicine [72,73], the identification of conserved NSC-activating genes in this simple model system offers a promising approach for developing therapeutic applications in vivo.

Supplementary Materials:
The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/cells11233727/s1, Figure S1: P. mytiligera de novo transcriptome assembly summary; Figure S2: Gene ontology categories enriched at the different stages of CNS regeneration; Figure S3: Expression profile for genes involved in nerve regenerative functions; Figure  S4: Expression profile for genes involved in regenerative functions; Table S1: Complete dataset of significantly differentially expressed genes (FDR < 0.05) and transcript counts at different time points of CNS regeneration (logFC and p-values); Table S2: Binary gene lists expressed; Table S3: Pathway enrichment expression along P. mytiligera CNS regeneration stages based on GeneAnalytics analysis.

Data Availability Statement:
The datasets generated for this study can be found in the NCBI Sequence Read Archive under accession number PRJNA773979.