DOI: 10.5061/dryad.4j0zpc8mj
title: Myelomeningocele spinal cord organoids scRNAseq

description: 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.
description: <strong>Single-cell RNA sequencing library
 preparation</strong> 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.
 <strong>Single-cell RNA sequencing analysis - Preprocessing
 & Quality Control</strong> 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. <strong>Single-cell RNA sequencing
 analysis – Cell Identification</strong> 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.
description: # **README: SCO scRNAseq** Data for *Spinal Cord Organoids from Human
 Amniotic Fluid iPSC Recapitulate the Diversity of Cell Phenotypes During
 Fetal Neural Tube Morphogenesis* Single-cell RNA sequencing data of spinal
 cord organoids generated from amniotic fluid cells collected from one
 control and two myelomeningocele fetuses. We have uploaded 3 stages of the
 data for this paper: 1. Raw FASTQ files generated & demultiplexed with
 Illumina's BaseSpace software. Organoid libraries prepped with
 Chromium Next GEM Single Cell 3' Reagent KIits v3.1 (Dual Index) and
 sequenced with Illumina NovaSeq 6000, targeting 50K reads/cell. 2. UMI
 Count Matrices from aligning FASTQ files to GRCh38 human genome reference
 assembly with 10X Genomics Cell Ranger 6.1.1. 3. Processed & clustered
 SeuratObject with identified cell types & metadata. ## **Description
 of the data and file structure** The four cell lines are labeled as
 follows: H29_5, H42_2, H47_4, and H47_62. Sample H42_2 is the Control, all
 other are MMC. For each cell line we have uploaded the following files:
 1\. Raw fastq files: `SCO-SAMPLEID_R1_001.fastq.gz` &
 `SCO-SAMPLEID_R2_001.fastq.gz`. For
 example: SCO-42-2_S21_L001_R1_001.fastq.gz,
 and SCO-42-2_S21_L002_R1_001.fastq.gz. These fastq files can be input into
 cell ranger for alignment and initial processing. 2\. The UMI count
 matrices: `SAMPLEID_matrix.mtx.gz`, `SAMPLEID_features.tsv.gz`, and
 `SAMPLEID_barcodes.tsv.gz`. For
 example: H42.2_matrix.mtx.gz, H42.2_features.tsv.gz,
 and H42.2_barcodes.tsv.gz. These files are the output of cell ranger and
 can be input into Seurat. 3\. The SeuratObjects: `SAMPLEID_named.Rds`. For
 example: h42_2_named.Rds. These are the fully processed and annotated
 Seurat files, and can be loaded in directly to generate graphs and DEGs.
 ## **Sharing/Access Information** This data is not uploaded elsewhere. ##
 **Code/Software** All analysis was performed using R code with the
 packages and functions described in the Methods section of the manuscript.