Time Course Transcriptome Analysis of Spina Bifida Progression in Fetal Rats

A better understanding of the transcriptomic modifications that occur in spina bifida may lead to identify mechanisms involved in the progression of spina bifida in utero and the development of new therapeutic strategies that aid in spinal cord regeneration after surgical interventions. In this study, RNA-sequencing was used to identify differentially expressed genes in fetal spinal cords from rats with retinoic acid-induced spina bifida at E15, E17, and E20. Gene ontology, KEGG, and protein–protein interaction analysis were conducted to predict pathways involved in the evolution of the disease. Approximately 3000, 1000 and 300 genes were differentially expressed compared to the control groups at E15, E17 and E20, respectively. Overall, the results suggest common alterations in certain pathways between gestational time points, such as upregulation in p53 and sonic hedgehog signaling at E15 and E17 and downregulation in the myelin sheath at E17 and E20. However, there were other modifications specific to gestational time points, including skeletal muscle development at E15, downregulated glucose metabolism at E17, and upregulated inflammation at E20. In conclusion, this work provides evidence that gestational age during spina bifida repair may be a significant variable to consider during the development of new regenerative therapeutics approaches.


Introduction
Myelomeningocele (MMC), the most significant form of spina bifida, is a devastating congenital malformation of the spinal cord associated with severe morbidity and mortality [1]. This "two-hit" process occurs during the folding of the neural plate into the neural tube during early development, 3-4 weeks of gestation, resulting in the lack of sclerotomal coverage leaving the neural tissue directly exposed to the amniotic fluid [2,3], called the "first hit" followed by an in utero acquired neurodegeneration by the chemical action of the amniotic fluid to the neural tissue or "second hit". This can result in severe consequences, including decreased mobility and limb paralysis, bladder and intestinal incontinence, and stunted neurological function [4]. Evidence indicates that consequences progress in severity if not corrected in utero; however, fetal repair only stops progression in most cases, as it does not reverse the existing damage [5][6][7].
Despite the advances in surgical techniques for spina bifida repair in utero, effective regenerative treatments for the devastating neurotological alteration have not yet been developed. To develop better therapeutic approaches, it is extremely important to understand the molecular changes present in the neural tissue once the defect has occurred and during this progressive degeneration in utero. The analysis of transcriptome studies at different time points throughout gestation would provide in-depth knowledge of the regulatory changes present in the neural tissue of the neurodegenerative progression in utero after ated using the RNA Nano6000 Assay Kit for the Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, CA, USA).
Three micrograms of total RNA per sample were used to generate sequencing libraries (n = 2 animals per group per time point). Libraries were generated using the NEBNext ® UltraTM RNA Library Prep Kit for Illumina ® (NEB, Ipswich, MA, USA) following the manufacturer's recommendations. cDNA products were purified using the AMPure XP system and library quality was determined using the Agilent Bioanalyzer 2100 system. The cBot Cluster Generation System using the HiSeq PE Cluster Kit cBot-HS (Illumina, San Diego, CA, USA) was used to cluster the samples according to the manufacturer's instructions. After this, the library products were sequenced on an Illumina HiSeq (Illumina, San Diego, CA, USA) and 125 bp/150 bp paired-end reads were generated per sample.

RNA-Seq Data Processing and Analysis
Raw FastQ files were processed through in-house Perl scripts, where clean reads were obtained by removing reads containing adapter or poly(N) as well as low-quality reads. TopHat v2.0.12 was used to align clean reads to the reference genome, which was built using Bowtie v2.2.3. HTSeq v 0.6.1 was used to calculate the fragments per kilobase of exon per million base pairs mapped (FPKM) based on the gene's length and reads mapped to that gene. Differential expression analysis was performed using the DESeqR software v1. 18.0. The resulting p values less than 0.05 were considered differentially expressed. This analysis identified differentially expressed genes (DEGs) between control, vehicle, and MMC groups at E15, E17, and E20. Principal component analysis was conducted using the AltAnalyzer package v2.0 (http://www.altanalyze.org, accessed on 18 February 2021) [13][14][15]. The GOseqR package was used to perform gene ontology (GO) enrichment analysis and GO terms with a p value less than 0.05 were considered statistically significant. RNAseq data were deposited in the Sequence Read Archive (SRA) and can be found via BioProject (PRJNA683230) and SRA(PRJNA683793). The protein-protein interaction network (PPI) was constructed and illustrated using the search tool for the retrieval of interacting genes/proteins (STRING) (https://string-db.org/, accessed on 27 May 2021) database to reveal the relationships of the top 25 DEGs based on a minimum required interaction.

RT-qPCR Analysis
Utilizing the RT 2 First Strand Kit (Qiagen, Germantown, MD, USA), 1 µg RNA /sample was reverse transcribed into cDNA. Four to six samples from the MMC and vehicle groups were analyzed. A 1 µg cDNA sample was then used as a template for RT-qPCR employing TaqMan ® gene expression assays (Applied Biosystems, Foster City, CA, USA) ( Table 1) in the 7500 Fast Real-Time PCR system. Samples were run in duplicate for target genes and were normalized using HPRT1 as an endogenous control.

Statistical Analysis
All graphs were performed in GraphPad Prism9 software (GraphPad Software Inc., La Jolla, CA, USA). A fold change > 1.5 and p value < 0.05 was considered statistically significant.
Relative quantification of transcript expression from RT-qPCR was performed using the 2 −∆∆Ct method comparing MMC and vehicle, where Ct represents the threshold cycle. Error bars indicate the standard error of the mean.

Differentially Expressed Genes in the Fetal Spinal Cord after MMC Occurs
To identify spinal cord changes in gene expression after MMC occurs in rat fetuses, we compared RNA-sequencing data between MMC, control, and vehicle groups at E15, E17, and E20. Principal component analysis indicates segregation of transcriptomes from the MMC group compared to the vehicle and control groups at each time point (Figure 1). The number of DEGs, identified as >1.5-fold change and p < 0.05 for each comparison at each time point is found in Table 2. Interestingly, the number of DEGs between MMC and either control or vehicle were greater than the control, compared to the vehicle, and the number of DEGs decreased between each comparison as gestation progressed. These trends are illustrated as volcano plots found in Figure 2

Pdgfa
Platelet-Derived Growth Factor Polypeptide Rn00709363_m1 *Probe codes from Thermo Fisher Scientific.

Statistical Analysis
All graphs were performed in GraphPad Prism9 software (GraphPad Software Inc., La Jolla, CA, U.S.A.). A fold change > 1.5 and p value < 0.05 was considered statistically significant. Relative quantification of transcript expression from RT-qPCR was performed using the 2−ΔΔ Ct method comparing MMC and vehicle, where Ct represents the threshold cycle. Error bars indicate the standard error of the mean.

Differentially Expressed Genes in the Fetal Spinal Cord after MMC Occurs
To identify spinal cord changes in gene expression after MMC occurs in rat fetuses, we compared RNA-sequencing data between MMC, control, and vehicle groups at E15, E17, and E20. Principal component analysis indicates segregation of transcriptomes from the MMC group compared to the vehicle and control groups at each time point ( Figure  1). The number of DEGs, identified as >1.5-fold change and p < 0.05 for each comparison at each time point is found in Table 2. Interestingly, the number of DEGs between MMC and either control or vehicle were greater than the control, compared to the vehicle, and the number of DEGs decreased between each comparison as gestation progressed. These trends are illustrated as volcano plots found in Figure 2         Volcano plot illustrating RNA-seq data of fetal spinal cords isolated from MMC and vehicle at (A) E15, (B) E17, and (C) E20. Red identifies genes that were greater than 1.5fold change and p < 0.05.   Hierarchical clustering analysis indicates that vehicle and control groups cluster together compared to the MMC group at each time point studied ( Figure 5). These data indicate that at E15 and E17, genes involved in neurological systems function, such as synaptic transmission, synapse, neuron projection, and neurotransmitter transport, are downregulated in MMC compared to either control or vehicle ( Figure 5A,B). Additionally, at E20, inflammatory genes, such as those involved in MHC class II antigen presentation, are upregulated, and genes involved in the development of the myelin sheath are downregulated in MMC compared to the control groups ( Figure 5C). Because of this clustering pattern, we chose to focus our further analyses on the comparison between MMC and vehicle at each time point; however, the comparisons between MMC and control and control compared to the vehicle can be found in the supplemental information. Hierarchical clustering analysis indicates that vehicle and control groups cluster together compared to the MMC group at each time point studied ( Figure 5). These data indicate that at E15 and E17, genes involved in neurological systems function, such as synaptic transmission, synapse, neuron projection, and neurotransmitter transport, are downregulated in MMC compared to either control or vehicle ( Figure 5A,B). Additionally, at E20, inflammatory genes, such as those involved in MHC class II antigen presentation, are upregulated, and genes involved in the development of the myelin sheath are downregulated in MMC compared to the control groups ( Figure 5C). Because of this clustering pattern, we chose to focus our further analyses on the comparison between MMC and vehicle at each time point; however, the comparisons between MMC and control and control compared to the vehicle can be found in the Supplemental Information.
The 25 genes that are most significantly upregulated and downregulated between MMC and vehicle groups at E15, E17 and E20 are listed in Tables 3-5 respectively. Additionally, the top 25 upregulated and downregulated DEGs between MMC and control at E15, E17, and E20 are listed in Supplementary Tables S1-S3 respectively. Furthermore, the top 25 upregulated and downregulated DEGs between control and vehicle at each time point are listed in Supplementary Tables S4-S6. Collectively, these results provide initial evidence that downregulated genes at E15 and E17 may lead to a decline in neural function and upregulated genes at E20 result in inflammation and disruption in myelin sheath development in rat fetuses with RAinduced MMC. The 25 genes that are most significantly upregulated and downregulated between MMC and vehicle groups at E15, E17 and E20 are listed in Tables 3, Table 4 and Table 5 respectively. Additionally, the top 25 upregulated and downregulated DEGs between MMC and control at E15, E17, and E20 are listed in Supplementary Tables S1-S3 respectively. Furthermore, the top 25 upregulated and downregulated DEGs between control and vehicle at each time point are listed in Supplementary Tables S4-S6.

GO Analysis of Differentially Expressed Genes
The top 10 enriched GO biological processes, cellular components, and molecular functions from upregulated and downregulated DEGs between MMC and vehicle at E15 are listed in Tables 6 and 7, E17 in Tables 8 and 9, and E20 in Tables 10 and 11. Additionally, the enrichment score of each GO listed between MMC and vehicle at E15, E17, and E20 is depicted in Figures 6-8, respectively. At E15, upregulated DEGs were enriched in biological processes, such as glial cell migration and regulation of mesoderm development, cellular components, such as contractile fiber and T-tubule, and molecular functions, such as extracellular matrix binding and caspase regulator activity (Table 6, Figure 6A). At E15, downregulated DEGs were enriched in biological processes, such as protein localization to synapse and synaptic vesicle maturation, cellular components, such as cell junction and ionotropic glutamate receptor complex, and molecular functions, such as calciumdependent protein binding and GABA receptor binding (Table 7, Figure 6B).       At E17, upregulated DEGs were enriched in biological processes, including proximal/distal pattern formation and negative regulation of cell fate specification, cellular components, including myosin complex and collagen, and molecular functions, such as extracellular matrix structural constituent, actin filament binding, and growth factor binding (Table 8, Figure 7A). At E17, downregulated DEGs were enriched in biological processes, such as cellular glucose homeostasis and negative regulation of tissue remodeling, cellular components, such as the ionotropic glutamate receptor complex and myelin sheath, and molecular functions such as GABA receptor binding and neuropeptide receptor binding (Table 9, Figure 7B).
At E20, upregulated DEGs were enriched in biological processes, such as cellular response to interferon-gamma, and regulation of cytokine production cellular components, such as keratin filament, and molecular functions, such as GTPase activity and cytokine receptor binding (Table 10, Figure 8A). At E20, downregulated DEGs were enriched in biological processes, such as regulation of the cellular ketone metabolic process and negative regulation of neurogenesis, cellular components, such as the extracellular matrix and neuronal cell body, and molecular functions, such as anion transmembrane transporter activity and voltage-gated ion channel activity (Table 11, Figure 8B).
The top 10 GO biological processes, cellular components, and molecular function enriched from DEGs that are upregulated or downregulated in MMC spinal cords, compared to the control at E15, E17, and E20, can be found in Figure S1. Several common GOs were identified when comparing MMC to either control group, but not when comparing control and vehicle groups ( Figure S2). For example, at E15, both comparisons indicated that the biological process innervation and molecular function neurexin binding were downregulated. Additionally, at E17, both comparisons indicated a downregulation in neuronal action potential propagation and upregulation in the extracellular matrix. Finally, at E20, antigen processing and presentation of peptide or polysaccharide antigen via MHC class II and MHC class II protein complex were upregulated and myelination was downregulated.

KEGG Pathway Enrichment Analysis of Differentially Expressed Genes
The top 10 upregulated and downregulated KEGG pathways were determined and reported as percent total changed genes in each pathway. All pathways identified were greater than 1.5-fold change and p < 0.05. Interestingly, there were several common pathways changed at E15 and E17 when comparing MMC and vehicle groups ( Figures 9A and 10A). For example, at both time points, p53 signaling and hedgehog signaling were upregulated; however, fewer genes were changed in these pathways at E17 ( Figure 10A). Additionally, at E17, metabolic pathways including glutathione and pyruvate metabolism were also downregulated and ECM receptor interaction was upregulated ( Figure 10A). At E20, immunological pathways, such as staphylococcus aureus infection, were upregulated, while metabolic pathways including PPAR signaling, ether lipid metabolism and glycerophospholipid metabolism, were downregulated ( Figure 11A). Additionally, at this time point, gap junction, and regulation of actin cytoskeleton were also downregulated ( Figure 11A).
The top 10 upregulated and downregulated KEGG pathways identified between the MMC and control groups and control and vehicle groups are found in Figures S3 and S4, respectively. Several common pathways were identified at each time point between MMC and either of the control groups. For example, the cell cycle and ECM receptor interactions were upregulated in both comparisons at E15 and E17. Furthermore, both comparisons demonstrated upregulation in autoimmune diseases, such as rheumatoid arthritis, type 1 diabetes mellitus, and autoimmune thyroid disease. Brain Sci. 2021, 11

Protein-Protein Interaction Network Analysis of Differentially Expressed Genes
To investigate the interactions between DEGS, protein-protein interactions (PPI) were determined for each time point by submitting the top 25 up and downregulated

Protein-Protein Interaction Network Analysis of Differentially Expressed Genes
To investigate the interactions between DEGS, protein-protein interactions (PPI) were determined for each time point by submitting the top 25 up and downregulated DEGs between vehicle and MMC groups to the string. At E15, we identified a cluster of interactions between upregulated DEGs involved in structure and skeletal muscle development including Acta1, Actn2, and Myog ( Figure 9B). Additionally, interactions between various downregulated transporters, such as Slc22a6 and Slc6a13, were determined as well as proteins involved in CNS cell fate determination, including Olig1 and S100b ( Figure 9C); however, these interactions were also identified when comparing MMC groups to the control group ( Figure S5). Furthermore, interactions between the same downregulated transporters were also identified at E17 in addition to proteins involved in neuronal cell identity, such as those in the Hox family, and neuron and axon maturation, including Calca and Nefh ( Figure 10C). Interestingly, similar interactions were also identified when comparing MMC and control groups ( Figure S6). Specific to the comparison between MMC and vehicle groups were the interactions between upregulated Dbx1 and neurog, both involved in neurogenesis, as well as Cebpa and Wt1, which are involved in the differentiation and survival of non-CNS cell types, including adipocytes and granulocytes ( Figure 10B). Furthermore, interactions between downregulated Uts2b, Gng13, CCL19, and Gng4 were demonstrated ( Figure 10C). At E20, PPI analysis indicated interactions between a host of downregulated DEGs involved myelination and oligodendrogenesis ( Figure 11C). Similar to interactions observed at E17, the majority of these interactions were also found when comparing MMC and control groups ( Figure S7); however, specific to the comparison of MMC and vehicle were the interactions between GPR17 and Bcas1 with Lims2 ( Figure 11C). Interactions involving upregulated DEGs associated with MHC class II immune reactions, CD74 and RT1-Da were identified at E20 when comparing MMC and the control group ( Figure S8B), with the addition of RT1-Bb and Eef1a1 in the comparison between MMC and vehicle groups ( Figure 11B).

Validation by qRT-PCR
To validate the sequencing results, the gene expression of 7-8 genes was randomly measured using qRT-PCR. The Log2fold change was determined between spinal cords isolated from MMC and vehicle groups at each time point and analyzed compared to results from the RNA sequencing ( Figure 12). While the absolute values were not identical between each measurement type, the trend changes between each gestational time point were congruent with the results obtained from RNA sequencing. These data support that the RNA sequencing method provided reliable quantifications.

Discussion
In this study, we employed RNA sequencing to characterize the chronology of most significant gene expression alterations in the spinal cord after defect developm and exposure to the amniotic fluid in the congenital retinoic acid rat model of spina b da. Models created from single or double mutants in mice do not always match the p notype of human neural tube defect with the same mutation [12]. For example, a mu tion in Vang-like protein 1 (VANGL1) was associated with neural tube defects in mans; however, there was not a neural tube defect phenotype found in homozygo mouse VANGL1 mutants [12]. It is suggested that human neural tube defects are a res of a combination of interacting mutations. The retinoic acid model is similar to hum MMC both developmentally and anatomically, so this model may be also associa with interacting mutations [16]. Therefore, the use of retinoic acid is a translationally evant model to study human MMC. By comparing the gene expression of spinal co isolated from fetuses with retinoic acid-induced spina bifida to those of its control s lings or fetuses from vehicle-treated dames, we identified over 3000 differentially pressed genes (DEGs) at E15, over 1000 DEGs at E17, and over 300 DEGs at E20. Us this well-established and cost-efficient teratogenic animal model and the matched in nal controls and vehicle controls allowed us to identify the changes during the seco hit MMC pathophysiology (i.e., the acquired degenerative effect after amniotic fluid posure) and not due to specific genetic mutations. We used this information to then p dict the potential role of the DEGs using GO analysis, KEGG pathway analysis, and p tein-protein interaction analysis at each time point. To our knowledge, this is the f study of its kind that identifies key DEGs and potential pathways that could be involv in the neural alterations of spina bifida in this model in utero. We propose that this formation lays a foundation for the further study of novel pathways that could pot tially be involved in the advancement of spina bifida and specific clinical outcomes. A

Discussion
In this study, we employed RNA sequencing to characterize the chronology of the most significant gene expression alterations in the spinal cord after defect development and exposure to the amniotic fluid in the congenital retinoic acid rat model of spina bifida. Models created from single or double mutants in mice do not always match the phenotype of human neural tube defect with the same mutation [12]. For example, a mutation in Vang-like protein 1 (VANGL1) was associated with neural tube defects in humans; however, there was not a neural tube defect phenotype found in homozygous mouse VANGL1 mutants [12]. It is suggested that human neural tube defects are a result of a combination of interacting mutations. The retinoic acid model is similar to human MMC both developmentally and anatomically, so this model may be also associated with interacting mutations [16]. Therefore, the use of retinoic acid is a translationally relevant model to study human MMC. By comparing the gene expression of spinal cords isolated from fetuses with retinoic acid-induced spina bifida to those of its control siblings or fetuses from vehicle-treated dames, we identified over 3000 differentially expressed genes (DEGs) at E15, over 1000 DEGs at E17, and over 300 DEGs at E20. Using this well-established and cost-efficient teratogenic animal model and the matched internal controls and vehicle controls allowed us to identify the changes during the second hit MMC pathophysiology (i.e., the acquired degenerative effect after amniotic fluid exposure) and not due to specific genetic mutations. We used this information to then predict the potential role of the DEGs using GO analysis, KEGG pathway analysis, and protein-protein interaction analysis at each time point. To our knowledge, this is the first study of its kind that identifies key DEGs and potential pathways that could be involved in the neural alterations of spina bifida in this model in utero. We propose that this information lays a foundation for the further study of novel pathways that could potentially be involved in the advancement of spina bifida and specific clinical outcomes. Additionally, these results may also lead to the future tailoring time-specific treatments in combination with the current standard of care that would potentially enhance current surgical spina bifida repair strategies as well as improve overall clinical outcomes.
Our results suggest that a spina bifida defect is associated with alterations in the spinal cord gene expression that regulates aspects of cell survival and positioning, neuron function, and skeletal muscle development during E15 and E17 gestational points. Compared to fetuses with normal development, genes involved in p53 signaling were upregulated in fetuses with MMC at E15, a time point similar to gestational weeks 5-6 in humans [17]. This potential role of p53 signaling is supported by previous evidence that p53 mRNA is upregulated at E15 in the spinal cord of rat fetuses with retinoic acid-induced spina bifida [18]. This pathway is likely involved in the extreme apoptosis found in the spinal cord during spina bifida, which consequently leads to neuronal cell loss, contributing to the impairment of neurological functions, such as motor skills [8,18,19]. Similarly, KEGG pathway analysis also predicted hedgehog signaling to be upregulated, which was previously implicated in neural tube defects [4]. Abnormal hedgehog signaling decreases the survival of neuronal precursors and alters the position of motor neurons resulting in abnormalities in the structure of the motor column [20,21]. Motor neuron function may also potentially be impacted through hedgehog signaling effects on cellular retinoic acid-binding protein 1 (CRABP1) in which variants were observed in patients with neural tube defects [22,23]. Additionally, the glutamate receptor complex was downregulated at this time point, which also diminishes neural cell migration [24]. Furthermore, our results support modifications in neuron function potentially due to downregulation in synaptic vesicle maturation and GABA receptor binding, which plays a role in synapse formation [25]. Finally, DEGs involved in skeletal muscle development and cytoskeleton were upregulated at this time point. Previous studies support the accumulation of actomyosin machinery, which increases tissue stiffness and mechanically inhibits the normal closure of the neural plate, leading to a neural tube defect [26][27][28].
Interestingly, p53 and hedgehog signaling were also upregulated in fetuses with retinoic acid-induced MMC at E17, a time point similar to gestational weeks 20-26 in humans [17]. This was also associated with the downregulation of the glutamate receptor complex, supporting that neuronal survival and positioning continued to be impacted during gestation of this model. However, during this time point, further alterations were observed that diminish neuron function and may lead to clinical outcomes such as paralysis and neurological dysfunction. For example, downregulations in metabolic processes, such as those described by GOs, such as glutathione metabolism, pyruvate metabolism, and cellular glucose homeostasis, were also evident. Deficiencies in molecules involved in glutathione, a major antioxidant, metabolism have been linked to neuronal cell loss in the brain and cognitive impairment, a significant clinical outcome in spina bifida patients [29,30]. In addition, abnormalities in glucose metabolism may lead to impairments in neuron function, as neurons not only break down glucose to meet energetic demands, but they also consume pyruvate that is released by astrocytes after glycolysis [31,32]. Preclinical and clinical studies indicate that diabetic mothers or mothers with alterations in genes that regulate glucose metabolism have an increased risk of neural tube defects, supporting a potential role for glucose metabolism in the progression of spina bifida [33,34]. Furthermore, at E17, GO analysis also indicated downregulation of processes involved in myelin sheath production, which is heavily controlled by oligodendrocytes [35]. This effect is likely due to dysregulation in hedgehog signaling at E15 and E17, as this pathway is a major regulator of oligodendrocyte differentiation and function [36].
The dysregulation of processes involved in metabolic pathways and myelination was also apparent in fetuses with retinoic acid-induced MMC at E20, a time point similar to gestational week 34 in humans [17]. In contrast to glutathione and glucose metabolism downregulated at E17, KEGG analysis indicated the downregulation of pathways involved in lipid metabolism, such as lipase activity, ether lipid metabolism, and glycerophospholipid metabolism later in gestation. While few studies have investigated a role for lipid metabolism in spina bifida progression, rare variants in lipid metabolism have been observed in patients with spina bifida phenotypes [37]. Additionally, modifications in this pathway have been connected to diseases associated with motor neuron loss, such as amyotrophic lateral sclerosis in which patients develop muscle paralysis [38]. Furthermore, while several DEGs involved in oligodendrogenesis and myelination were downregulated as a result of retinoic acid, but not associated with spina bifida, other such DEGs including Gpr17, Lims2, and Bcas1 were downregulated in the spinal cord of fetuses with MMC defects. These specific genes were also significantly reduced in prenatal neurosphere cultures after BMP2 treatment, indicating the potential role of BMP2 in diminishing myelination during the progression of spina bifida [39]. Finally, specific only to E20, our results provide evidence that inflammatory processes are upregulated in the spinal cord of fetuses with retinoic acid-induced MMC. Our previous data further support this claim, as an increase in activated microglia, characterized as Iba1+ and MHC class II+ cells, IL-1B, IL-6, and IFN-g, was observed in the spinal cord of fetuses with MMC at E20, but not at earlier time points [8].
We hypothesize that the changes in gene expression we observed after the spina bifida defect is at least in part due to exposure of the spinal cord to the amniotic fluid. Studies on amniotic fluid supernatant collected from pregnant women at the time of open defect identified alterations in pathways similar to our results including those associated with neuronal development, axonal development, and synapse formation [40]. Interestingly, the inflammation had some of the most prominent alterations [40]. Furthermore, this inflammatory response may be due to the toxic components of amniotic fluid, as the amniotic fluid of rats with retinoic acid-induced MMC contained higher amylase levels and activity compared to healthy controls [41]. While the specific role of amylase in spinal cord inflammation and spina bifida outcomes has not been investigated, serum amylase is elevated during pancreatitis and is associated with the elevation of pro-inflammatory cytokines, such as Il-1B and IL-6, which supports the potential pro-inflammatory role of amylase in amniotic fluid [42].

Conclusions
In conclusion, through a comprehensive time-course transcriptomic analysis, this is the first study that characterized the progressive changes that occur in the neural tissue after exposure to the amniotic fluid in utero, using the congenital retinoic acid-induced spina bifida rodent model. Current standard surgical strategies only structurally repair the defect and inhibit further neurological damage due to amniotic fluid exposure. Likely due to limited understanding of the mechanisms that drive spina bifida degeneration in utero, there are no treatments available that regenerate healthy neural tissue and function after spina bifida diagnosis. Our results provide evidence that different mechanisms may play more important roles during specific periods throughout fetal progression. Therefore, it may be beneficial to tailor new therapeutic strategies to the gestational age at the time of treatment as well as to entertain an approach where a combination of pathways is targeted. Furthermore, future studies further elucidating the specific targets at different gestational ages during disease progression should be a focus in the spina bifida research community.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/brainsci11121593/s1, Figure S1: Top 10 GO biological processes, cellular components, and molecular functions enriched from DEGs that are upregulated (A-C) or downregulated (D-F) in MMC spinal cords compared to control at E15 (A,D), E17 (B,E), and E20 (C,F). Figure S2: Top 10 GO biological processes, cellular components, and molecular functions enriched from DEGS that are upregulated (A-C) or downregulated (D-F) in vehicle spinal cords compared to control at E15 (A,D), E17 (B,E), and E20 (C,F). Figure S3: Top 10 KEGG pathways upregulated and downregulated between MMC and control at E15 (A), E17 (B), and E20 (C) as determined by percent total changed genes in each pathway. Figure S4: Top 10 KEGG pathways upregulated and downregulated between vehicle and control at E15 (A), E17 (B), and E20 (C) as determined by percent total changed genes in each pathway. Figure S5: Protein-protein interaction network analysis based on top 25 downregulated (A) and upregulated (B) differentially expressed genes between MMC and control at E15. Figure S6: Protein-protein interaction network analysis based on top 25 downregulated (A) and upregulated (B) differentially expressed genes between MMC and control at E17. Figure S7: Proteinprotein interaction network analysis based on top 25 downregulated (A) and upregulated (B) differentially expressed genes between MMC and control at E20. Figure S8: Protein-protein interaction network analysis based on top 25 downregulated (A) and upregulated (B) differentially expressed genes between vehicle and control at E15. Figure S9: Protein-protein interaction network analysis based on top 25 downregulated (A) and upregulated (B) differentially expressed genes between vehicle and control at E17. Figure S10: Protein-protein interaction network analysis based on top 25 downregulated (A) and upregulated (B) differentially expressed genes between vehicle and control at E20. Table S1: Top 25 upregulated and downregulated genes between MMC and control groups at E15. Table S2: Top 25 upregulated and downregulated genes between MMC and control groups at E17. Table S3: Top 25 upregulated and downregulated genes between MMC and control groups at E20. Table S4: Top 25 upregulated and downregulated genes between control and vehicle groups at  E15. Table S5: Top 25 upregulated and downregulated genes between control and vehicle groups at  E17. Table S6: Top 25 upregulated and downregulated genes between control and vehicle groups at E20. Graphical Abstract (created with BioRender.com).