Myelomeningocele (MMC) is a severe form of spina bifida associated with substantial neurologic morbidity. In vitro modeling systems of human spinal cord development may help to elucidate the underlying pathophysiology of the MMC spinal cord. We developed spinal cord organoids (SCO), defined as self-organized, three-dimensional clusters of spinal tissue, that were derived from human amniotic fluid induced pluripotent stem cells. Here, we used a variety of analyses, including immunofluorescent and single-cell transcriptomic approaches, to compare SCOs from healthy and MMC fetuses. Organoids contained a diverse range of neural and mesodermal phenotypes when cultured for up to 130 days in vitro. Multielectrode arrays revealed functional activity with evidence of emerging neuronal networks. Fetal spina bifida modeling was successfully established by culturing MMC SCOs in second- and third-trimester amniotic fluid for 3 weeks. Taken together, we show that functional SCOs can recapitulate the cellular identity of the fetal spinal cord and represent a novel research platform to study the interplay between cellular, biochemical, and mechanical cues during human MMC neural tube morphogenesis.

Single-cell RNA sequencing library preparation

scRNA-seq libraries were generated using the Chromium Next GEM Single Cell 3ʹ Reagent Kits v3.1 (Dual Index). Per manufacturer’s protocol, a maximum volume of 43.3 µL/sample was used for processing to target up to 10,000 cells. Cells were combined with RT reagents and loaded onto 10X Next GEM Chip G along with 3’ v3.1 gel beads. Approximately 100µL of emulsion was retrieved from the chip and incubated (45 min at 53°C, 5 min at 85°C, cool to 4°C), generating barcoded cDNA from each cell. Samples were cleaned using 0.6X SPRIselect beads. Then 10uL of amplified cDNA was carried into library preparation. Fragmentation, end repair, and A-tailing were completed and samples were cleaned up with SPRIselect beads. Adaptor ligation was followed by a 0.8X cleanup, and then amplification was performed with PCR. Libraries were sequenced on the Illumina NovaSeq 6000 using v1.5 kits, targeting 50K reads/cell. Demultiplexing and FASTQ generation was completed using Illumina’s BaseSpace software.

Single-cell RNA sequencing analysis - Preprocessing & Quality Control

UMI count matrices were generated by aligning fastq files to GRCh38 human genome reference assembly with the 10x Genomics Cell Ranger 6.1.1 count function [48]. Preprocessing and QC filtering were then performed with the R package Seurat 4.2.1 [49]. The Read10X and CreateSeuratObject functions were used to convert count matrices into a Seurat object. Samples 29-5 (MMC1) and 47-4 (MMC2.1) were filtered to remove cells with gene counts below 200 or over 8000, and over 15% mitochondrial expression. Sample 42-2 (Control) was filtered to remove cells with gene counts below 200 or over 7500, and over 20% mitochondrial expression. Sample 47-6.2 (MMC2.2) was filtered to remove cells with gene counts below 200 or over 7500, and over 15% mitochondrial expression. Quality control cutoffs were based on violin plots showing feature count and mitochondrial percentage distribution for each sample, cutoffs tailored to each sample to properly remove dead cells and doublets while maximizing input. After filtering the cell and gene counts for each sample were as follows: 29-5; 4752 cells expressing 26193 genes, 47-4; 5193 cells expressing 26311 genes, 42-2; 9543 cells expressing 25581 genes, 47-6.2; 3621 cells expressing 25570 genes. Next, the data was cell-cycle scored and log normalized and the top 2000 variable genes were identified with the FindVariableGenes function. This function uses the variance stabilizing transformation method of identifying top variable genes with a regularized negative binomial regression model. We then scaled the data, regressed out cell cycle score, and then used PCA for dimensional reduction. After using an elbow plot to visualize the percentage variance explained by each principal component, the top 15 principal components were utilized for clustering and UMAP reduction.

Single-cell RNA sequencing analysis – Cell Identification

We identified cell types through a combination of manual and computational techniques. First, we loaded the matrix and metadata on the developing mouse spinal cord and converted the gene IDs to human using biomaRt [36, 50]. We then created a Seurat object from the Delile data and transferred cell labels using the FindTransferAnchors and TransferData functions. We looked at percentage of each assigned label per cluster to get broad cell identification of neuron, skin, mesoderm, and progenitor. Next, to determine more specific cell subtypes we performed differentially expressed gene (DEG) analysis and gene set enrichment analysis (GSEA). Additionally, we examined collections of marker genes to identify neuronal and mesoderm cell subtypes. Between these different methods of analysis, we were able to identify both genes and pathways to distinguish glia, mesoderm, progenitor, glial progenitor, neural progenitor, neuron, and actively dividing cell population as clusters within our organoids and renamed them accordingly. We were also able to define the mixed cell populations within these clusters, as genes and pathways indicate the presence of both astrocytes and oligodendrocytes within our glial populations and the presence of both immature and mature neurons in the neuron population. We determined rostro-caudal domains by

performing a specific TransferData of the Neuron clusters to identify dorsal and ventral neurons. To examine cervical, brachial, thoracic, and lumbar regionalization, we analyzed HOX genes and their overlap through Seurat feature plots.