Single-cell sequencing of individual retinal organoids reveals determinants of cell fate heterogeneity

Abstract

With a critical need for more complete in vitro models of human development and disease, organoids hold immense potential. Their complex cellular composition makes single-cell sequencing of great utility; however, the limitation of current technologies to a handful of treatment conditions restricts their use in screens or studies of organoid heterogeneity. Here, we apply sci-Plex, a single-cell combinatorial indexing (sci)-based RNA-seq multiplexing method to retinal organoids. We demonstrate that sci-Plex and 10x methods produce highly concordant cell class compositions and then expand sci-Plex to analyze the cell class composition of 410 organoids upon modulation of critical developmental pathways. Leveraging individual organoid data, we develop a method to measure organoid heterogeneity, and we identify that activation of Wnt signaling early in retinal organoid cultures increases retinal cell classes up to six weeks later. Our data show sci-Plex's potential to dramatically scale-up the analysis of treatment conditions on relevant human models.

Figure 1.. Comparison between sci-Plex and 10X.

A) Retinal organoids were cultured for 29, 78, or 185 days before being dissociated into single-cell suspensions. The cells were split and either processed by a standard 10x RNA-seq pipeline or prepared by sci-Plex. In sci-Plex, the cells are lysed in age-specific wells, leaving the nuceli intact. Single-stranded polyadenylated DNA oligos, called “hashes”, are added to each well. The hashes contain well-specific barcode sequences and have an innate affinity for the nuclei which helps recruit the oligos to the surface of the nuclei where we permanently crosslink the hashes to the nuclei. After crosslinking, the nuclei from all wells are pooled and subjected to sci-RNA-seq. B) Integrated UMAP of sci-Plex and 10x datasets derived by Seurat’s CCA-MNN methodology. Left: Cells are colored by cell type. RPC: retinal progenitor cell, MG: Muller glia, Npre: retinal neuronal precursor, RGC: retinal ganglion cell, RGCpre: RGC precursor, HApre: horizontal/amacrine precursor, PBpre: photoreceptor/bipolar precursor. Right: Cells are faceted by organoid age and technology. Cells are colored by age. C) Dot plot of the genes used to define cell types in B. The size of the dot indicates the percent of cells that express the gene of interest. The color indicates the log10 mean UMIs per cell. D) A representation of the established view of retinal neuron developmental trajectories. E) The timeline of retinal cell differentiation in human organoids. F) Immunostaining of D28, D63, and D189 organoids for progenitor (VSX2-D28), RGC (POU4F2), photoreceptor/bipolar (OTX2), neuronal (MAP2), amacrine/horizontal/RGC (HuC/D), bipolar cells (VSX2-D189) and photoreceptor (RCVRN) markers. Scale bar represents 100 μm. G) Using the cell type assignments from the integrated dataset, cell types were split into Brain/Non-neuronal (Non-neuronal), RGC+Horizontal+Amacrine (RGC+H+A), and Photoreceptor+Bipolar (PR+BP) categories. Cells were counted and the % Cell Type was determined for each technology and organoid age.

Figure 2.. sci-Plex applied to individual organoids after signaling pathway modulation.

A) An overview of the retinal organoid differentiation protocol and treatment timing. On Day 6, BMP4 was added to the cultures. BMP4 was diluted in half on Day 8 and again on Day 10 before being removed completely on Day 12. SAG, SU5402, and/or CHIR were added on Day 14–18. The organoids were cultured for an additional 10 days before collection. B) A summary of the seven treatments performed in this experiment. Note that the All, BMP, CHIR, SU5402, and SU5402:CHIR treatments include BMP4. C) A box plot of the number of cells recovered per organoid for each treatment. The number above the boxplot indicates the number of individual organoids in which > 50 cells were recovered for that treatment. D) Dot plot of the genes used to define cell types for D28 organoids. The size of the dot indicates the percent of cells that express the gene of interest. The color indicates the log10 mean UMIs per cell. E) Left: UMAP generated by Monocle3 of the recovered cells with cell type annotations from D28 organoids treated with signaling modulators. Right: For each of the UMAPs only the cells from the specified treatment are colored by treatment. RPC_OV: retinal progenitor cell/optic vesicle, RPE: retinal pigmented epithelium, OS: optic stalk, RP: roof plate, FP: floor plate. F) A stacked bar plot in which the mean cell type composition of the D28 organoids are displayed for each treatment. G) Violin plots of size-factor normalized cell counts for each cell type from D28 organoids. Plots are colored by treatment.

Figure 3.. sci-Plex allows the detection of significant changes in cell type abundance.

A) Heatmap of fold change in cell type abundance compared to “None” for D28 organoids. A beta-binomial model was fit to the data and used to determine whether there were significant changes in cell type abundances. A in the center of the box indicates a q-value of less than 0.05 after Benjamini Hochberg correction. RPC_OV: retinal progenitor cell/optic vesicle, RP: roof plate, OS: optic stalk, RPE: retinal pigmented epithelium. B) Heatmap of fold change in cell type abundance compared to “BMP” for D28 organoids. A beta-binomial model was fit to the data and used to determine whether there were significant changes in cell type abundances. A* in the center of the box indicates a q-value of less than 0.05 after Benjamini Hochberg correction. C) Immunostaining of D28 organoids with VSX2 (RPCs/Npre, green), PAX2 (OS, red), and DAPI (blue). Scale bar represents 50 μm. Arrow = RPCs/Npre, empty arrow head = non-optic-stalk (OS), filled arrow head = OS. D) Quantification of D28 immunostained organoids. Between 3–7 organoids were stained per condition (All = 5, BMP = 3, CHIR = 5, None = 5, SAG = 7, SU5402 = 5, SU5402:CHIR = 5). Significance was determined by performing t-tests between the treatment of interest and None. Benjamini Hochberg correction was performed to adjust the p-value. * represents a p < 0.05.

Figure 4.. Detection of retinal organoid heterogeneity and organoid subtypes.

A) Bright phase images of organoids the day before they were collected for sci-flex (D27). Scale bars are 200 μM. B) Cell type legend for C-I. C) Stacked bar plot in which each bar represents an individual D28 organoid exposed to the None treatment. The bars are colored by the percent a given cell type makes up the composition of the individual organoid. Same as C, but for the D) All, E) BMP, F) SU5402, G) CHIR, H) SU5402:CHIR, and I) SAG treatments. J) A UMAP generated using the size-factor normalized cell type by organoid matrix generated from the D28 individual organoid dataset. Organoids are colored by cluster. K) The UMAP from J with the organoids colored by treatment. L) The mean and coefficient of variation (CV) were modeled for each cluster and each treatment using a gamma distribution. Displayed are the modeled relationships colored by treatment. M) Output of likelihood ratio tests (Irtest) with or without inclusion of the treatment term for D28 organoids. N) Model from L faceted by treatment. The shaded area marks the range within 2 standard deviations from the mean. Points represent the individual cluster values.

Figure 5.. D63 organoids demonstrate continued perturbation by signaling modulation.

A) UMAP from Monocle3 with cell type annotations of all cells recovered from individually hashed D63 organoids exposed to BMP, CHIR, and SU5402:CHIR treatments. RPC: retinal progenitor cell, RPCprolif: proliferating RPC, RPCproneural: proneural RPC, Npre: retinal neuron precursor, HApre: horizontaVamacrine precursor, RGC: retinal ganglion cell, RGCpre: RGC precursor, PBpre: photoreceptor/bipolar precursor, ONHprog: optic nerve head progenitor, RPE: retinal pigmented epithelium, CN: cortinal neuron, RG: radial glia, CP: choroid plexus. B) Heatmap of fold change in cell type abundance compared to BMP. A * in the center of the box indicates a q-value of less than 0.05 after Benjamini Hochberg correction. C) Violin plots of size-factor normalized cell counts for D63 organoids. Dark gray = the treatment that was used for fold-change calculations in B, light gray = no significant change in cell type abundance, teal = decrease in cell type abundance, orange = increase in cell type abundance. D) Stacked bar plots in which each bar represents an individual organoid exposed to the BMP, CHIR, or SU5402:CHIR treatment. The bars are colored by the percent a given cell type makes up the composition of the individual organoid. Right: cell type color legend. E) Immunostaining of cryosectioned D63 organoids with OTX2 (photoreceptors/RPE/CP, blue), MAP2 (retinal and cortical neurons, green), POU4F2 (RGC, red), and DAPI (gray). Arrow = Non-retinal neurons, filled arrow head = OW+ retinal tissue, empty arrow head = Pou4f2+ retinal tissue. Scale bar represents 100 μm. F) Violin plots of % Area measurements for RPE as calculated from bright phase microscopy images of D62 organoids. Plots are colored by significance as in C. G) Violin plots of size-factor normalized RPE cell counts for D63 organoids. Plots are colored by significance as in C. H) Representative brightfield images of D62 organoids from each group quantified in G-H. Arrow = non-retinal tissue, filled arrow head = retinal tissue, empty arrow head = RPE. Scale bar represents 200 μm.