Longitudinal scRNA-seq of retinal organoids derived from Stargardt disease patient with ABCA4 mutation

ABCA4 mutations represents the most prevalent genetic cause of inherited retinal diseases in outbred populations, accounting for more than 17% of cases¹. The ABCA4 gene encodes a crucial transport protein that eliminates toxic bisretinoids from photoreceptor cells, once disrupted, leads to the toxic metabolite accumulation and progressive retinal degeneration in STGD1, the predominant form of STGD². The clinical landscape of STGD1 is particularly complex due to ABCA4’s extensive genetic heterogeneity, with over 2,200 variants identified across 50 exons, resulting in diverse disease onset and severity pattern based on the affected protein domains³.

While mouse models have been extensively developed to study ABCA4-associated pathology, they present significant limitations in replicating human disease characteristics. The primary constraint is the absence of a macula in mice, the distinct anatomical structure that serve as the primary site of pathology in STGD1 patients⁴. Furthermore, mouse models with disease-causing ABCA4 missense mutations exhibit distinct patterns of protein expression and localization that may inadequately represent human pathophysiology⁵. These complexities underscore the critical need for ideal model systems and innovative approaches to elucidate the mechanistic relationships between specific missense mutations and their pathogenic consequences in STGD1.

Induced pluripotent stem cells (iPSCs) have emerged as an invaluable tool for modeling inherited retinal diseases, particularly for investigating the diverse spectrum of disease-causing genetic variants. The successful establishment of multiple STGD1-specific iPSC lines has provided renewable and reproducible cellular resources for mechanistic investigations^6,7,8,9,10. These iPSC platforms offer two complementary modeling approaches: the generation of three-dimensional retinal organoids (ROs) for assessing variant pathogenicity¹¹, and the derivation of RPE cells¹². Together, these approaches enable comprehensive examination of disease mechanisms within an authentic human genetic context.

In this study, we established RO models using patient-derived iPSCs harboring the most prevalent ABCA4 variants identified in the Chinese population (Asn965Ser/Arg18Pro) (Fig. 1a). Single-cell RNA sequencing was employed to characterize the transcriptional landscape during retinal organoid development, enabling detailed analysis of cell type specification and developmental trajectories (Fig. 1b,c). The development time points were selected based on previously established correlation between healthy iPSC-derived retinal organoids and human developing retina¹³, and were validated through immunostaining of representative cell types and light microscopy (Fig. 2a). Following technical validation of dataset quality, unsupervised clustering analysis identified photoreceptor cell populations that closely correlated with STGD1 pathological changes.

The study workflow of ROs generation and scRNA-seq. (a) Generation of Patient-specific ROs. Patient-derived iPSCs were differentiated into RO following an established protocol. The timeline illustrates critical developmental stages and corresponding culture conditions throughout the differentiation process. (b) Single-Cell RNA sequencing. Single cell suspensions were processed using 10X Chromium Controller for microfluidic encapsulation in Gel Beads-in-Emulsions (GEMs). Following encapsulation, single-cell cDNA libraries were constructed and subsequently sequenced. (c) Bioinformatic pipeline. The computational workflow encompassed multiple sequential steps: initial data quality control, reference genome alignment, unique molecular identifier (UMI) counting, gene and cell filtering, data normalization, and cell clustering analysis.

Full size image

Morphological Characterization of ABCA4-Mutant ROs. (a) Brightfield and immunofluorescence images of retinal organoids at different develop-mental stage (d40-d200). Marker used: CRX (photoreceptor precursors), HuC/D (RGC), Ki67(proliferating cells), Chx10(retinal progenitors), ARR3 (cone photoreceptors), NRL (rod photoreceptors). DAPI stains for nuclei. (b) Immunofluorescence images showing L/M-opsin (green) for outer segment, HSP60 (red) for inner segment, and DAPI (blue) for nuclei staining in ABCA4 MT-RO (top row) and control HiPS-RO (bottom row) at day 200. Scale bar: 50 μm. (c–f) Immunofluorescence staining for ARL13B (green) as connecting cilia, HSP60 (red) as inner segment, and DAPI (blue) in ABCA4 MT-RO (c,d) and HiPS-RO (e,f). Panels d and f show higher magnification views of outlined in c and e, respectively. Scale bars: 50 μm (c,e), 10 μm (d,f). (g) Schematic diagram of a cone photore-ceptor, illustrating the localization of key proteins: L/M opsin in the outer segment (OS), ARL13B in the connecting cilium (CC), and HSP60 in the inner segment (IS).

Full size image

This dataset represents the first single-cell transcriptomic analysis of STGD1, revealing significant transcriptional heterogeneity within photoreceptor cells caused by ABCA4 mutations. Our advances in long-term ROs culture have established improved disease models with enhanced translational relevance compared to conventional cell lines and rodent models. This comprehensive dataset not only enables detailed investigation of STGD1 pathological mechanisms but also provides a valuable platform for therapeutic drug discovery.

Human subject

PBMCs were isolated from blood samples collected from a 22-year-old male patient diagnosed with Stargardt disease (STGD1). The patient provided written informed consent to participate in the study, including specific consent for sharing and publishing the genetic data in public databases. The patient was fully informed about the potential risks associated with sharing genetic information. The study was conducted in compliance with the Declaration of Helsinki and approved by the ethics committee of the Eye and ENT Hospital of Fudan University (Approval ID: 2019025). The experimental protocols were also reviewed and approved by the local ethics committee, ensuring that all procedures adhered to the appropriate ethical standards and guidelines.

Isolation of PBMCs

Peripheral blood samples were processed using Ficoll (cytiva, Cat No. 17-1440-03) and centrifuged to isolate PBMCs. The isolated mononulcear cells were subsequently seeded on to 24-well plates (500 μL of complete medium per well) and cultured in PBMC Medium, which consisted of StemSpan^TM SFEM II (STEMCELL Technologies, Cat No.0-9655) supplemented with a cytokine cocktail containing 100 ng/ml stem cell factor (SCF), 20 ng/ml interleukin-3 (IL-3), 20 ng/ml interleukin-6 (IL-6), 100 ng/ml FMS-like tyrosine kinase 3 ligand (Flt-3L). To maintain the cells in an undifferentiated state and promote proliferation, half of the expansion medium was replaced daily for a total of 4 days.

Reprogramming and characterization of induced pluripotent stem cells

Following the expansion period, the cells were subjected to the reprogramming process to generate induced pluripotent stem cells (iPSCs). iPSCs were generated from PBMCs by CytoTuneTM Sendai Reprogramming Kit (Thermo Fisher Scientific, Cat No. A16517) according to the manufacturer’s instructions. Briefly, 5 × 10^5 PBMC were plated in a well of 6-well plate in PBMC Medium and transduced with appropriate multiplicity of infection (MOI) of Sendai virus overnight. On day 6 post-transduction, the cells were plated onto Matrigel-precoated plates and cultured with mTeSR^TM Plus medium (STEMCELL Technologies, Cat no.100-0276). The colonies emerged between day 15–21 post-transduction and were manually picked for expansion. The commercially available iPSC line (DYR0100) was used as a healthy control. The established iPSCs were maintained in mTeSR Plus medium at 37 °C and 5% CO2 with daily medium changes. Passaging was performed every 4–5 days using ReLeSR (STEMCELL Technologies, Cat No.0-5872). The generated iPSCs were thoroughly evaluated through a comprehensive set of characterization assays: (1) alkaline phosphatase staining to confirm pluripotency status (Fig. S1c), (2) karyotype analysis to ensure chromosomal stability (Fig. S1e), (3) immunostaining of key pluripotency markers including NANOG,SSEA4,TRA-1-81 (Fig. S1d), (4) assessment of in vitro differentiation capacity through immunostaining of markers for all three embryonic germ layers (Fig. S1f), and (5) teratoma formation assay to verify in vivo differentiation potential (Fig. S1g). The iPSC line successfully passed these quality control assessments, confirming its pluripotent status and suitability for disease modeling.

Generation of Retinal Organoids from induced Pluripotent Stem Cells

ROs were generated following a modified protocol described previously¹⁴. iPSCs were expanded in Matrigel-coated dishes using mTeSR medium. On day 0, iPSCs colonies were dissociated into small pieces to form embryoid bodies (EBs) and transferred to Neural Induction Medium (NIM) containing DMEM/F12 (Thermo Fisher,11330032), 1%N2 supplement (Thermo Fisher, 17508-048), 1x MEM nonessential amino acids (NEAA) (Thermo Fisher, 11140050), 1x GlutaMAX (Thermo Fisher, 35050061) and 2 μg/ml heparin (STEMCELL Technologies, 7980), 10 μM Y27632(abmole, M1817). Y27632 was removed on day 1, and the medium was changed every other day until day 6. On day 6, 1.5 nM BMP4 (R&DSystems,314-BP-010) was added to fresh NIM. EBs were plated on Matrigel at a density of 200 Ebs/well in a 6-well plate on day 7. Half the NIM was replaced on days 9, 12 and 15. On day 16, the medium was switched to Retinal Differentiation Medium (RDM) containing DMEM: F12 3:1, 2% B27 supplement (Thermo Fisher, 17504-044), MEM NEAA, 1% penicillin/streptomycin). 3D ROs emerged on days 25–30 and were dissected with a surgical knife. ROs were maintained in low attachment 3.5-cm dishes and fed twice weekly feeding with 3D-RDM consisting of DMEM: F12 3:1, 2% B27 supplement, 1 × MEM NEAA, 1 × antibiotic-antimyocotic (Thermo Fisher, 15240-062) and 5% FBS, 100 μM taurine, and 1:1000 chemically defined lipid supplement (Thermo Fisher, 11905031). 1 μM all-trans retinoic acid (Sigma, R2625) was supplemented until day 100.

Immunostaining and imaging

Immunofluorescence was performed either on fixed cells or on retinal organoid cryosections. Slides were then incubated for 1 h with a blocking buffer (PBS, 5% bovine serum albumin, 0.5% Triton X-100) at RT, and then incubated with primary antibody with appropriate dilution at 4 °C overnight. On the next day, accordingly secondary antibodies conjugated with either AlexaFluor-488, 594 or 647 diluted appropriately were adopted. To counterstain the nuclei, slides washed with PBS were then incubated with DAPI diluted at 1:1000. Finally, slides were washed two times in PBS before mounting with Fluoromount-G.

Fluorescent staining signals were captured with Leica confocal microscope. Control of the confocal microscope and image acquisitions were conducted by using the Olympus Fluoview software version 4.2 at a resolution of 1024 × 1024 pixels, with a scan rate of 8–10 μs pixel−1. Images were acquired sequentially, line by line, in order to reduce excitation and emission crosstalk, step size was defined according to the Nyquist-Shannon sampling theorem. Exposure settings that minimized oversaturated pixels in the final images were used. Twelve-bit images were processed with ImageJ or FIJI, Z-sections were projected on a single plane using maximum intensity under the Z-project function.

Single-cell RNA sequencing of retinal organoids

Preparation of single-cell suspension

ROs were harvested at multiple timepoints during the reprogramming process (day 40, 90, 150, 200, and 260) and subjected to PBS washes twice at 37 °C. To dissociate the organoids into single cells, an activated papain solution was prepared by combining papain (23.2 μl/ml) with a mixture of 1.1 mM EDTA, 5.5 mM cysteine-HCl, 0.067 mM 2-mercaptoethanol, and 0.005% DNase I and incubating at 37 °C for 30 minutes prior to use. The ROs were then incubated in the activated papain solution for 10 minutes at 37 °C to facilitate digestion. The enzymatic digestion was terminated by the addition of 4 volumes of ‘stop solution’ (DMEM/F12 supplemented with 5% BSA). The digested samples were centrifuged at 800 rpm for 10 minutes at 4 °C, and the resulting cell pellet was resuspended in 2 ml of 1640 culture medium. The centrifugation was repeated once to ensure the removal of any residual debris. The supernatant was carefully aspirated, and the cell pellet was resuspended in 100 μl of 1640 medium, yielding a single cell suspension suitable for scRNA-seq.

Single-cell Sequencing

Single-cell RNA sequencing libraries were prepared using the Chromium Controller Instrument and Single Cell 3′ V3.1 Reagent Kits (10X Genomics, Pleasanton, CA). Single-cell suspensions were adjusted to a concentration of approximately 1000 cells/µL and encapsulated into Gel Bead-In-Emulsions (GEMs) for barcoded reverse transcription. Following GEM dissociation, the barcoded cDNA was purified and amplified. Library construction proceeded with fragmentation, A-tailing, adapter ligation, and index PCR amplification. Library quality was assessed using the Qubit High Sensitivity DNA assay (Thermo Fisher Scientific) for quantification and the Bioanalyzer 2200 High Sensitivity DNA chip (Agilent) for size distribution analysis. Sequencing was performed on an Illumina NovaSeq 6000 platform (Illumina, San Diego, CA) using 150 bp paired-end configuration.

Quantification and statistical analysis

Generating single-cell gene expression matrices

Single-cell RNA sequencing data analysis was performed using the NovelBrain Cloud Analysis Platform (NovelBio Co., Ltd., www.novelbrain.com). Raw sequencing data underwent quality control using fastq (v0.21.0) with default parameters to remove adapter sequences and low-quality reads. Quality filtering criteria excluded reads containing >5% ambiguous bases (N) and reads with >40% low-quality bases (quality score <15). The filtered clean data were aligned to the human reference genome (GRCh38 Ensemble: version 104) using CellRanger (v7.1.0)¹⁵ with include-introns mode to generate feature-barcode matrices.

Quality control and filtering

The matrices were converted to Seurat objects individually using the Seurat R package (version 5.1.0)¹⁶. The data were filtered to remove cells with fewer than 200 unique genes per cell or greater than 10000 genes per cell, or more than 15% reads mapped to mitochondrial genes. Filtered gene matrices were normalized with the ‘NormalizedData’ function using the ‘LogNormalize’ method. The top 2000 variable genes were identified using the ‘vst’ method in the ‘FindVariableFeature’ function¹⁷.

Integration, dimensional-reduction and clustering

Batch effects between samples were corrected using the standard integration algorithm with “FindIntegrationAnchors” function in Seurat to minimize the influence of technical variability. Dimensionality was reduced using principal component analysis (PCA). The optimal number of principal components (15) was determined through elbow plot analysis, jackstraw procedure (p < 0.05), and clustree-based cluster stability assessment. In this way, a neighborhood graph was calculated using these 15 principal components. Following a Uniform Manifold Approximation and Projection (UMAP) embedding, cell clusters were detected based on the Louvain clustering implementation using a resolution parameter of 0.25.

Cell-type annotation

To facilitate the cell type annotation, DGEs was calculated by using ‘FindAllMarkers’ function on normalized count data to identify markers of each cluster. Cell identities were assigned based on canonical markers assembled from relevant literature¹³.

All the original data have been deposited to the NCBI Sequence Read Archive (SRA), with the number of SRP554832¹⁸. Additionally, the immunofluorescence data mentioned in the methods are included in the manuscript as part of the results section, specifically within Fig. 2 and Fig. S1.

Morphological characterization of STGD1 RO models

We firstly validate the successful establishment of RO through morphological observation during development. The ROs exhibited three distinct developmental stages consistent with previous reports: Stage 1 (Day 40) showed continuous, phase-bright outer neuroepithelial rim; Stage 2 (Day 90) developed a phase-dark core with a reduced phase-bright outer rim; and Stage 3 (Day 150 and Day 200) displayed characteristic hair-like surface appendages (Fig. 2a). Cell-type specific immunostaining confirmed the presence and distribution of key retinal cell populations, including RGCs (HuC/D), photoreceptor precursor (CRX), proliferating progenitor (Ki67, Chx10), cones (ARR3), and rods (NRL) (Fig. 2a).

Comparative analysis between ABCA4-mutant and control ROs revealed distinct morphological alterations. Notably, ABCA4-mutant ROs demonstrated significantly reduced cone populations compared to control ROs (Fig. 2b), while maintaining comparable inner nuclear layer composition as indicated by HSP60 staining. This observation aligns with the previous findings showing opsin inhibition by Stargardt disease byproduct¹⁹. Intriguingly, synaptic structures in ABCA4-mutant ROs exhibited punctate appearance (Fig. 2c,d) rather than the classic “horseshoe” morphology observed in controls (Fig. 2e,f), suggesting that ABCA4 mutations substantially impact cone morphology, particularly affecting outer segment and ciliary structural integrity.

Single-Cell RNA sequencing quality control

We performed scRNA-seq analysis on ROs from five developmental time points (Fig. 1b). Sequencing metrics demonstrated high data quality (Table 1), with each sample achieving 300 million reads, barcodes efficiency of 97.1-98.5%, UMI Q30 percentages of 93.4–96.7%, and genome mapping rates of 92.8–94.1%. Quality control metrics (Table 2) confirmed minimal technical variation. From 43011 cells analyzed, median gene counts ranged from 3349 to 4000 per cell, with median UMI counts of 7050 to 9622. Following doublet removal and quality filtering, the final dataset demonstrated robust cell numbers and feature counts (Fig. 3a–c). We acknowledge that our bioinformatic analysis pipeline follows established Seurat workflows employing linear PCA and the FindIntegrationAnchors method for batch correction. While these methods provide reliable analysis for our current dataset description, we recognize that alternative approaches (such as non-linear dimensionality reduction techniques, newer batch correction algorithms like Harmony or LIGER, or alternative differential expression methods) may offer different perspectives on the data. To facilitate re-analysis using these or other emerging methodologies, we have made the raw data readily accessible and provided detailed documentation of our analytical pipeline in the code repository.

Quality Control for scRNA-seq of ABCA4-mutant ROs. (a) Cell numbers distribution across developmental time points (D40-D260) after quality filtering, showing the proportion of filtered cells, doublets, and low-quality cells removed. (b) Violin plots showing the distribution of detected genes per cell across different developmental stages (D40-D260). (c) Violin plots depicting the distribution of UMI (Unique Molecular Identifier) counts per cell across developmental stages. (d) UMAP dimensional reduction plot showing 14 distinct cell clusters (0–13) identified through unsupervised clustering analysis. (e) UMAP visualization of cells colored by developmental time points (D40-D260), demonstrating the temporal progression of cellular populations during retinal organoid development. Each panel represents a specific time point, revealing distinct cellular distribution patterns.

Full size image

Cell type identification and developmental trajectory

UMAP dimensionality reduction revealed cell distributions consistent with developmental progression (Fig. 3d). The cellular composition evolved from predominantly retinal progenitor cells at day 40, through emergence of ganglion cells, amacrine cells and photoreceptors at days 90–150, to rod and cone predominance at days 200–260 (Fig. 3e). Cell type identification using established markers: RPC (SFRP2, FOXP1), photoreceptor precursor (OTX2, PRDM1), RGC (GAP43, SNCG), AC (TFAP2A, PTF1A), HC (ONECUT2, ONECUT3), NRPC (ATOH7, GADD45A), bipolar (PRKCA, VSX2), rod (NRL, PDE6A), cone (PDE6H, THRB), RPE (BEST1, TYR), muller (RLBP1, GLUL) (Fig. 4a). Cell type-specific gene expression patterns were clearly delineated (Fig. 4b,c).

Cell type composition and marker gene expression of ABCA4-mutant ROs. (a) Stacked bar plot showing the proportional distribution of different cell types across developmental time points (D40-D260). The plot reveals the dynamic changes in cellular composition, including early retinal progenitor cells (RPC), photoreceptor precursors (PR_pre), cones, Müller glia, neurogenic retinal progenitor cells (NRPC), amacrine cells (AC), retinal ganglion cells (RGC), rods, and retinal pigment epithelium (RPE). (b) UMAP feature plots showing the expression patterns of key retinal cell type markers. Expression intensity is indicated by color gradient (blue to red). RCVRN and THRB (photoreceptor markers), NRL (rod photoreceptor marker), ISL1 (ganglion cell marker), TFAP2A (amacrine cell marker), ONECUT3 (horizontal cell marker), GAP43 (ganglion cell marker), and RPE65 (RPE marker) demonstrate distinct distribution patterns. (c) Dot plot showing the expression levels of selected marker genes across major retinal cell types. Dot size represents the percentage of cells expressing the gene, while color intensity indicates the average expression level.

Full size image

The scRNA-seq data processing pipeline were run on the Linux operating system. All R source code used for downstream data analyses and visualization are provided online (https://figshare.com/articles/journal_contribution/Code/27999230, where https://doi.org/10.6084/m9.figshare.27999230 is the DOI)²⁰.