Understanding the complexity of retina and pluripotent stem cell derived retinal organoids with single cell RNA sequencing: current progress, remaining challenges and future prospective.

Single-cell sequencing technologies have emerged as a revolutionary tool with transformative new methods to profile genetic, epigenetic, spatial, and lineage information in individual cells. Single-cell RNA sequencing (scRNA-Seq) allows researchers to collect large datasets detailing the transcriptomes of individual cells in space and time and is increasingly being applied to reveal cellular heterogeneity in retinal development, normal physiology, and disease, and provide new insights into cell-type specific markers and signaling pathways. In recent years, scRNA-Seq datasets have been generated from retinal tissue and pluripotent stem cell-derived retinal organoids. Their cross-comparison enables staging of retinal organoids, identification of specific cells in developing and adult human neural retina and provides deeper insights into cell-type sub-specification and geographical differences. In this article, we review the recent rapid progress in scRNA-Seq analyses of retina and retinal organoids, the questions that remain unanswered and the technical challenges that need to be overcome to achieve consistent results that reflect the complexity, functionality, and interactions of all retinal cell types.

Introduction

Single-cell RNA sequencing (scRNA-Seq) is a powerful tool for exploring the transcriptional heterogeneity of large cell populations by measuring gene transcription in each individual cell. Elucidating the cellular heterogeneity within a biological system is an essential requisite for understanding development and evolution, how a biological system is homeostatically regulated, and responds to external stimuli. Today, with the rapid development of the single-cell-based technologies it is possible to determine the mRNA levels of several thousand transcripts in thousands of cells as well as map their spatial distribution through spatial transcriptomics techniques.[1] The first paper which described the feasibility of single cell-RNA sequencing was published 10 years ago where the authors profiled only eight cells.[2] Only 7 years later, 10× Genomics released a data set of more than 1.3 million cells (https://support.10xgenomics.com/single-cell-gene-expression/datasets). The number of high impact studies has increased rapidly during the last few years, and in addition to transcriptional data, chromatin accessibility profiles can now be obtained at the single-cell level, enabling gene regulation studies and identification of cell type transcription factors and their targets.[3] ScRNA-Seq is able to detect and distinguish the transcriptome of individual cells; thus, the similarities and differences between cells can be studied. Previously unknown heterogeneity has been observed with distinct and rare cell populations being revealed[4,5] in neurons[6], embryos,[7,8] and immune cells.[9,10] Relationships between cell types can also be inferred allowing for insights into developmental progression and lineage tracing.[8,11-15] Furthermore analysis of gene expression may reveal gene regulatory networks that control cell specificity, fate, and function.[16,17]

These rapidly emerging technologies also enable the detection and characterization of specific cell populations related to human development, health, and disease.[18] In parallel to the technological developments, multiple methods of bioinformatics analysis have been developed, including read quantification, quality control, normalization, clustering, cell annotation using reference maps, and identification of differentially expressed genes, cell surface markers, and ligand-receptor pairs. There has been rapid development of methods to reduce noise, improve sensitivity and, notably the throughput for single-cell transcriptomics.[19-23] To date, scRNA-Seq has been applied to mouse, primate, human fetal and adult retina,[24,25] and to retinal organoids generated from pluripotent stem cells.[26-28] In this review, we take the opportunity to review this fast progress and summarize the current challenges and unanswered questions for the field of retinal development and disease.

scRNA-Seq technologies

There are a number of methods that have been developed for scRNA-Seq with differing approaches to cell capture, cDNA amplification, and library generation in order to isolate single cells and reduce biases.[5,23,29-38] The method used has an effect on the level of information obtained, such as the number of genes detected, number of transcripts per gene, expression of genes of interest, and resolution of differential splicing.[39] Despite this diversity in scRNA-Seq techniques, they all rely on the basic principles of viable single-cell preparation, isolation, lysis, capture of polyadenylated mRNA, reverse transcription, cDNA amplification, cDNA library preparation, and next-generation sequencing (Figure 1).

Figure 1.

Schematic presentation of single-cell capture, lysis, RNA/cDNA preparation and amplification, and library generation.

Cells can be isolated in a number of ways. Tens to hundreds of cells can be captured in wells of a multi-well plate by fluorescence-activated cell sorting or smaller “capture sites” using microfluidics: this allows nanoliter volumes to be used, which can improve sensitivity.[40] Alternatively, droplet-based platforms such as Drop-Seq[37] and inDrop[36] are able to capture thousands of cells in microdroplets, each containing a bead with a unique barcode and molecular identifier allowing generation of many scRNA-Seq libraries.[39] The number of cells and the methods employed for reverse transcription (RT), cDNA amplification, and sequencing library generation affects the ability to detect relative differences of transcript expression and reveal the biological variation between cells. However, as the amount of isolated RNA is relatively small from a single cell, scRNA-Seq is also prone to technical variability.[16,21] This variation can occur during the conversion of a particular RNA into cDNA and its subsequent presence in the sequencing library, as a fraction of mRNAs is lost resulting in dropouts, or be affected by the fidelity of the amplification step introducing noise, or how accurately the mRNA input correlates to the number of reads.[41-43] The Smart-Seq method was developed to optimize the RT and amplification steps, and full-length transcript coverage. Further improvement led to Smart-Seq.2[32,33] These methods are able to generate full-length cDNA libraries facilitating greater transcript coverage and detecting more genes per cell than alternative counting methods that employ unique identifying sequences.[43] A large portion of scRNA-Seq methods utilize these additional sequences in the RT primers. This reduces the transcript coverage but allows the inclusion of such sequences as barcodes, unique molecular identifiers, and adaptor sequences to identify a single RNA molecule and the cell of origin, as well as multiplexing and NGS library generation. The Unique Molecular Identifier (UMI) tags the mRNA at the RT stage prior to amplification allowing certain protocols such as CEL-Seq, inDrop, MARS-Seq, Drop-Seq, and SCRB-Seq to increase the throughput of library generation.[5,23,36-38] This can also distinguish cDNA and library amplification duplicates[44] and improve the quantification of mRNA transcripts,[35,42] which has been utilized in STRT-Seq,[35] CEL-Seq,[42] CEL-Seq2,[34] Drop-Seq,[37] inDrop,[36] MARS-Seq,[5] and SCRB-Seq.[38] Amplification-associated noise can also be reduced by in vitro transcription (IVT) rather than PCR amplification, which is utilized in CEL-Seq2 and MARS-Seq.[42,43] Therefore, full-length Smart-Seq methods have fewer dropouts but greater amplification noise due to the use of PCR amplification. Methods utilizing IVT amplification (CEL-Seq2 and MARS-Seq) or UMIs (SCRB-Seq, CEL-Seq2, Drop-Seq, and MARS-Seq) have less amplification-associated noise.[42,43] STRT-Seq enriches for the 5ʹ end of mRNA. CEL-Seq, CEL-Seq2, MARS-Seq, SCRB-Seq enrich for the 3ʹ end. All incorporate cell-specific barcodes and UMIs, facilitating pooling of cDNA for library generation, shortening the procedure. MARS-Seq increases the CEL-seq2 method throughput through the use of a liquid-handling platform.[5] If the aim is the quantification of transcriptomes from a large number of cells with a low sequencing depth then droplet-based approaches, e.g., Drop-Seq, are recommended. Whereas other methods such as SCRB-Seq and Smart-Seq2 are preferable for the quantification of fewer cells and greater sensitivity.[43]

Miniaturization of the CEL-seq2 and Smart-Seq reactions to nanoliter volumes, as demonstrated by chip-based microfluidic systems, such as the Fluidigm platform, can improve sensitivity over standard scRNA-Seq.[45] The commercialization of these methods with proprietary hardware such as the Fluidigm C1 platform, as well as a number of droplet-based platforms, such as Chromium from 10x Genomics, ddSEQ from Bio-Rad Laboratories, InDrop from 1CellBio, and μEncapsulator from Dolomite Bio/Blacktrace Holdings is facilitating robust scRNA-Seq methodology for the masses. An alternative approach to scRNA-Seq is the isolation of single nuclei (sn) for snRNA-Seq. Studies have shown that despite the reduced number of transcripts from nuclei there is sufficient number to type them into broad classes of cells. Isolation of single nuclei may have some advantages over single cells as they are potentially less prone to any dissociation induced transcriptional changes and can be more easily isolated from complex and frozen tissues.[46-48]

Computational methods and challenges

Single-cell RNA-Seq measures gene expression at the cellular level, meaning that distinct gene expression profiles of rare cell types are not masked by average expression. This gives the potential to answer questions that cannot be addressed using bulk RNA-Seq analysis. The analysis of such datasets can be used to identify cell populations using statistical clustering methods, to study changes from one developmental time point to another and pinpoint key regulatory genes.

Alignment and quantification

The analysis begins with the quantification of RNA by alignment of reads to a reference genome to create a gene by cell expression matrix. This process is very similar to bulk RNA-Seq analysis and many of the same tools are applicable to single-cell experiments. However, some specialized tools such as STARsolo which is an extension of the popular aligner STAR,[49] and Alevin, which is part of the Salmon toolkit are available for quantification of the reads detected. Additionally, a number of pipelines are available such as CellRanger,[50] which is distributed by 10x genomics for analysis of 10x datasets and DropEst[51] which can be used for the analysis of data from other platforms.

After the expression matrix has been created, the analysis methods start to deviate from bulk RNA-Seq analysis. Single-cell data are fundamentally different from bulk data and many of the assumptions made by statistical methods designed for bulk analysis do not hold true.[52] Single-cell data are sparse, with many genes either not detected or detected at very low levels; there are no replicates as each cell can only be measured once and the data is inherently noisy and prone to variation caused by technical artifacts. These qualities mean that a different analysis approach is required. Since 2015, the number of tools and analysis approaches has grown rapidly and there are now a rich array of methods, which can be applied to this data, the majority of which are written in R or Python. scRNA-tools[53] is a comprehensive database cataloging the available methods. The omictools.org repository and ‘awesome-single-cell’[54] site also catalog a number of single-cell tools. The Bioconductor repository[55] hosts numerous packages for single-cell analysis as well as number of workflows explaining how these packages can be used in combination. In June 2019 Luecken and Theis[56] published a detailed article outlining the best practices for single-cell analysis. A general workflow for processing single-cell data is explained in Figure 2. Briefly, the workflow consists of a pre-processing step involving quality control checks, normalization, and dimension reduction followed by statistical analysis to identify relationships between the single cells and key regulatory genes. These processes are described in more detail below.

Figure 2.

Bioinformatics pipeline used for scRNA-Seq analysis.

Pre-processing – quality control

The analysis of the expression matrix begins with a quality control step. Dead cells or debris can obscure the interpretation of downstream analysis, so the gene expression must be filtered to remove problematic data. Poor quality data can be detected using metrics such as the number of genes detected per cell or the total counts per cell compared to the other cells in the experiment and filtering these cells from the analysis. The Scater[57] and Seurat[58] R packages and the ScanPy[59] Python package provide tools and a framework for performing these initial quality control steps. Another source of technical noise comes in the form of doublets, where two or more cells are captured in the same bead or droplet. A number of tools now exist to predict doublets so that they can be filtered from the data.[60-63]

Normalization and dimension reduction

Technical factors such as the number of genes detected per cell or the total number of reads per cell can contribute a significant proportion of cell-to-cell variability. Such factors must therefore be removed from the data by an appropriate normalization method before it is possible to measure the relationships between cells.[64,65] scRNA-Seq data can be normalized by estimating a size factor[66] that is used to scale the data or using a probabilistic method to model the counts.[64,67,68]

Where multiple samples are present in an analysis, a batch correction step may then be applied to the data. Batch correction methods designed for bulk samples can confound biological differences between cells with technical differences. A number of approaches have been designed for single-cell data. These include Mutual Nearest Neighbors[69] (MNN), Canonical Correlation Analysis[70] (CCA), and Harmony.[71]

After normalization and removal of batch effects, it is then necessary to reduce the dimensions of the dataset. This is because many of the genes measured show similar expression profiles across all cells so that despite the large number of data points captured for each cell, the distance between cells becomes very small, a phenomenon known as the “Curse of dimensionality”. To overcome this problem, the number of data points must be reduced so that only relevant sources of variation contribute to the interpretation of the data. This is achieved by a combination of feature selections where only the most biologically informative genes are used in the analysis and dimension reduction techniques. One method of feature selection that is commonly used is selecting highly variable genes. These are genes which show more variability than would be expected from the average expression of that gene.

Dimension reduction involves projecting the expression matrix into a low dimension space using techniques such as principal component analysis (PCA). These steps not only reduce the noise in the datasets but also help to reduce the computation time.[72]

Cluster analysis

One of the key goals of scRNA-Seq is to identify cell types and subpopulations of cells within a dataset. This can be achieved using clustering analysis. Multiple tools exist for identifying clusters of cells in single-cell data, in fact there are currently 136 analysis tools listed in the scRNA-tools database. Some of the more popular packages include Seurat[58,70] a graph-based method written in R, and Backspin – a biclustering[73] method written in Python. Kiselev et al[72] published a comprehensive review with advice on selecting the most appropriate methods for single-cell clustering.

As many clustering methods involve estimating the number of clusters a priori or by using tuning parameters to set the resolution of the clusters, multiple clustering solutions can be found. Assessment and interpretation of the clusters require an understanding of the biology of the cells being studied. It is worth remembering that clusters can consist of true cell types (for example, Rods or Cones) or cell states, which reflect transitory functional differences between the cells. Clusters can be validated statistically[74] by measuring stability of clusters using methods such as silhouette analysis[75] or the elbow method. One useful tool for this task is Clustree[76] which is designed to compare and evaluate multiple clustering solutions. Due to the complexity of scRNA-Seq datasets, it is useful to visualize in lower dimensional space – and the results of cluster analysis can be visualized using t-distributed stochastic neighbor embedding (t-SNE)[77] or Uniform Manifold Approximation and Projection for Dimension Reduction (UMAP).[78]

Trajectory analysis

Clear groups of cells may not always be present within a sample. During development, cells may be transitioning from one state to another and it is therefore preferable to analyze this journey rather than perform clustering analysis, which forces partitions onto the data. Trajectory analysis can be used to study these continuous processes in cells. The analysis takes advantage of the asynchronous nature of cell development, and each cell is treated as a snapshot of the developmental pathway. The process begins by defining a path through the data and the cells are then ordered along this pathway. This creates a so-called pseudotime ordering to the cells in the experiment, allowing inferences to be made about the developmental trajectory of a process. More than 70 tools[79] have already been developed for trajectory analysis capable of modeling linear, bifurcating or tree-shaped, or more complex topologies. Choosing the correct tool for a particular process is not an easy task. A comprehensive review of the available software was made by Saelens et al. in April. 2019[79] It concluded that the performance of a particular method was linked to the trajectory type being studied and that this should be used to guide the selection of the most appropriate method. For example, a linear trajectory may be sufficient when modeling the development of a single-cell type whereas a more complex tree-shaped trajectory may be necessary to analyze developmental processes of multiple cells.

Identification of differentially expressed genes

Differential expression analysis may be used to identify gene signatures for clusters. General statistical tests such as Wilcoxon rank-sum test as well as specialized methods, such as MAST,[68] ZINBWAVE,[80] and BaSICs,[81] which take into account zero inflation and high levels of technical noise specific to single-cell data can be applied, whereas for trajectory analysis a regression analysis or graph-autocorrelation analysis can be used to identify genes which change across a trajectory.

Cell type annotation and cell atlases

Whilst the above statistical methods provide insight into the level of heterogeneity and number of sub-populations of cells within a sample, to fully understand the data, the complex gene expression profiles must be annotated with known morphology or locations of cell types. This is a challenging process, which requires matching of marker genes detected by differential gene expression with known gene signatures for specific cell types, typically found through intensive literature searches. The annotation must then be validated by immunocytochemistry or other spatial methods.

In order to speed up this lengthy process, the Human Cell atlas project identified the need to create a detailed molecular map of every cell type in the human body, which could be used as a reference for future research. Cell atlases have the potential, when combined with cluster transfer tools such as scMap[82] and Garnet,[83] to dramatically speed up the cell annotation process.

Application of single-cell sequencing to studies of fetal and adult retina: cell type and marker identification

The retina is the innermost layer that lines the back of the eye and is vital for light sensing and image processing. The retina is derived from a germinal zone in the optic vesicle in which neuroepithelial cells proliferate to give rise to the six principle types of retinal neuronal cells (cones, rods, bipolar, amacrine, horizontal and ganglion cells) and one glial cell type (Müller glia), organized within three nuclear layers. All these cell types derive from retinal progenitor cells in an orderly spatio-temporal manner that has been well studied in vertebrates[84-86] and more recently in humans. Most of the retinal cells are sub-classified into subtypes; however, the molecular signature of each cell subtype remained unknown until the recent application of single-cell sequencing methods. The first study by Macosko et al. in 2015[87] used Dropseq to capture and analyze 44,808 cells from the retina of 14 old day mice, which identified 39 different cell clusters. This analysis was able to identify the known retinal cell types based on the specific expression of well-established retinal markers. Importantly, other cell types corresponding to astrocytes, microglia, endothelial cells, and fibroblasts were also identified and the relative abundance of each cell type was similar to previous estimates of microscopy–based techniques. Amacrine cells are the most diverse morphologically neuronal class being identified as inhibitory (utilizing GABA or glycine as neurotransmitter), excitatory (releasing glutamate) or amacrine cells with an unidentified neurotransmitter (nGnG amacrine cells). Macosko et al. were able to identify 21 clusters of amacrine cells in the adult mouse retina: of those, 12 were identified as GABAergic, 5 as glycinergic, 1 as excitatory, and 3 had low levels of GABAergic, glycinergic and glutamatergic markers, and hence these could represent the nGnG amacrine cells. Selective markers were then identified for each of the 21 amacrine cell subtypes, with the majority being further validated by immunohistochemistry.

Previous morphological and functional criteria have suggested the presence of around 100 cell types in the retina: this stands at odds with the above study, which identified only 39 cell types. This could be due to the number of cells captured and analyzed and the sequencing depth, which was perhaps insufficient to reveal diversity in the rare cell types (for example, retinal ganglion cells, which comprise <1% of retina). To overcome this, Rheaume et al. captured[88] and analyzed 6225 retinal ganglion cells (RGCs) from the left- and right – eye mouse retina, using immuno-panning with the pan-RGC marker, Thy1. This analysis identified 40 RGC subtypes, corroborating a recent electrophysiological study[89] and showed that all RGCs expressed the RBPMS marker, albeit at varying level. Most of the RGC subtypes were distributed proportionally between the left and the right eye; however, three clusters showed preferential distribution, indicating that some RGC related functions may be predominant in one eye. Using differential gene expression analysis, the study was able to identify subtype-specific markers, some of which were validated by fluorescence in situ hybridization (FISH). Interestingly the authors were also able to identify a new marker (Zic1) for the RGCs enriched in the right eye as well as transcription factor combinations for all the 40 RGC subtypes, which are documented in the RGC subtypes gene-browser: https://health.uconn.edu/neuroregenerationlab/rgc-subtypes-gene-browser.

Similarly, Shekhar et al.[90] used Dropseq to study mouse retinal bipolar cells (BCs) at single-cell resolution. BCs comprise around 7% of the retinal cells and are classified into rod and cone bipolar cells depending on the photoreceptor type from which they receive their synaptic input and further categorized into ON and OFF types based on the response to increases or decreases in light levels. Using a transgenic line that expresses GFP into All BCs and Müller Glia cells, the authors were able to capture and analyze around 25,000 cells from which 15 BC clusters could be identified, including all the know cell types and two novel populations, one of which displayed a non-typical BC morphology. The authors were also able to relate candidate marker expression to BC types using a combination of lentiviral labelling and FISH, which facilitated the assignment of new BC clusters, for example, BC1B with amacrine like morphology, which could have been incorrectly identified as an amacrine cell cluster, if only morphological-based classifications would have been applied.

Recently, Lukowski et al. were able to capture and analyze 20,009 cells from three adult human donor retinas and to identify 18 transcriptionally different cell populations.[91] The authors were able to show that there were no obvious variations in all major retinal cell types obtained from retina 6–14 hours post death with the exception of rod photoreceptors, which appeared to be more sensitive to degeneration. The single-cell sequencing data enabled the authors to identify a population with low MALAT1 expression, which correlated to degenerating rods. In an independent study, Voigt et al.[92] were able to capture 8,217 cells from three human adult neural retina obtained within 5.5 hours post death. Planned and reliable procurement of human adult retina is impossible; hence, authors cryopreserved the disassociated cells, which enabled stringent inclusion criteria and processing of samples obtained at different time points within one single cell-sequencing run. The authors were able to identify all the known retinal cell types, endothelial cells, and microglia within 17 clusters and to provide molecular markers, which should be of immense value for optimizing retinal organoid differentiation protocols as well as identification of novel candidate disease – genes. Both of these studies identified less cell clusters in human retina compared to the mouse. This is unlikely to be due to different cell type composition between human and mouse retina, but more likely to the deeper data mining within each cell cluster in the mouse retina, revealing cell type sub-specificity resulting in higher cluster number.

Voigt et al. were able to perform a regional difference expression analysis by characterizing 3578 cells from the fovea and 4639 cells from the peripheral human retina. The fovea is responsible for high acuity color vision and is characterized by a pit in which the inner retinal layers are displaced peripherally. Over the pit center, cone photoreceptor density is the highest in the retina, reaching over 200,000 cones/mm2. Rods are absent from the central 300um, which is called the rod free zone. Blood vessels are also absent forming the foveal avascular zone. A previous study has used bulk RNA-Seq analysis to characterize the development of the human fovea compared to the peripheral region;[85] nonetheless, Voight et al. provide the first regional single-cell analysis, identifying a total of 148 differentially expressed genes, 23 of which are enriched in foveal cells. Interestingly when the top 10 most enriched foveal and peripheral genes were compared to the single-cell study performed in M. fascularis,[93] 85% similarity was observed, indicating the very high similarity in regional differences in retina transcriptome between humans and macaques. However, when the single-cell RNA-Seq data were compared between macaque and mouse retinal cell types, a tight correlation was found for photoreceptors, BCs, and amacrine cells; nonetheless, the RGCs were divergent between the two species, despite displaying conserved transcription factor codes. These findings are of great interest, and while suggesting a reliance on the mouse models for studies related to photoreceptor BC and amacrine cells function, they do question the validity of this model for investigating the function of genes related to RGC and function.

While the majority of the studies described above have focussed on adult human retina, a recent study published by Hu et[25] al has performed single-cell sequencing of 2421 cells from human fetal neural retina and the retinal pigment epithelium, encompassing the 5th until the 24th gestational week. Twenty-one clusters were identified comprising all the known retinal cell types as well as microglia, fibroblasts and blood cells. By analyzing peak marker expression, this study documented the emergence of RGCs and retinal progenitor cells (RPCs) at 5 weeks, followed by horizontal cells at 9 weeks and the later emergence of amacrine, photoreceptors, bipolar and Müller glia cells, corroborating bulk RNA-Seq studies of developing human retina.[85,86] Importantly, the study identified cell-type specific transcription factors and their targets, which should facilitate immunohistochemical detection of these cell types both during in vivo and in vitro retinogenesis. Through mapping of the expression of known retinal disease genes to each cluster, the authors of this study were able to identify each gene to retinal-specific cells, thus shedding more light on the pathology of inherited retinal dystrophies.

Lineage tracing studies have shown that RPCs are multipotential, exhibit great variations in size, composition and division mode and can generate several if not all retinal cell types.[94-99] RPC fate is predominantly determined by the intrinsic cellular program with fine-tuning from extrinsic/environmental cues.[100-102] To better understand the transcriptional changes associated with changes in RPC competency, Clark et al.[103] performed a single-cell RNA-Seq of the developing mouse retina. Whilst this study identified a significant differential expression between RPCs at different stages, there were no sub-clusters within each individual age, indicating that RPCs undergo profound transcriptional changes during retinal development, which are consistent with the changes in their competency. The study identified embryonic day 16–18 as primary subdivision for early and late RPCs, with selective genes expressed in either of these two RPC subsets. The pseudotime analysis identified Nfi factors as being highly expressed within late-stage RPCs, bipolar and Müller glia cells. Accordingly, overexpression of NF1A/B/X resulted in an increase in the late-emerging retinal cells (BCs and Müller glia cells), while conditional loss of function alleles of Nfi genes in RPCs leads to a disruption of Müller glia marker staining and the outer limiting membrane. Together these new findings demonstrate that single-cell sequencing has the potential to identify new regulators of retinal histogenesis, in addition to cell type identification and generation of cell type markers (Figure 3).

Figure 3.

scRNA-Seq analysis of retina and pluripotent stem cell-derived retinal organoids: current insights.

scRNA-Seq of pluripotent stem cell-derived retinal organoids: new insights into their cellular composition and cell fate determination

Three-dimensional (3D) structures combining multiple cell types, which interact with each other, named organoids have been used in classical developmental biology experiments since the early 1960s. These 3D organoids can be generated from pluripotent, adult and cancer stem cells and are widely used for studies of disease modeling, gene editing, drug testing/repurposing, and cell-based replacement therapies due to their ability to mimic human development and organogenesis.[104] To date, several protocols have been developed to optimize generation of various organoids including those of gut, kidney,[105,106] liver,[107,108] brain,[109] and retina. Traditionally the study of these organoids has relied on bulk transcriptomics and/or immunofluorescence characterization of various cell types based on known markers with the latter comprising a laborious and time-consuming technique, which is often difficult to implement on multiple time points and experimental groups. The power of scRNA-Seq to provide transcriptomic information of thousands of cells simultaneously has revolutionized the field of organoids and has started to generate new insights into their composition and cell maturation[110] as well as novel gene expression signatures as shown recently by Harder et al. and Czerniecki al., for human pluripotent stem cell-derived kidney organoids. Equally, scRNA-Seq has been used to compare 2D and 3D culture conditions for generation of liver bud-like structures from human pluripotent stem cells, showing that those generated under 3D culture conditions have a striking similarity to fetal liver cells.[107] Multiple studies have also applied scRNA-Seq to brain organoids and have provided evidence for the existence of progenitors and differentiated cells of both neuronal and mesenchymal lineages, resembling the fetal neocortex.[111] These studies have revealed that cerebral organoids have the ability to self-organize alongside the dorsal-ventral front and to give rise to mature astrocytes and oligodendrocytes in additional to neuronal lineages.[112] In some reported instances these brain organoids contain retina. Quadrato et al. performed scRNA-Seq on these photosensitive brain organoids and identified Müller glia, photoreceptors, retinal ganglion, bipolar, and amacrine cells of the neural retina as well as retinal pigment epithelium (RPE).[113]

In a remarkable tour de force, Sasai and his group generated 3D aggregates from mouse embryonic stem cells, which under minimal culture condition went on to form optic vesicles and later optic cups displaying at later stages of differentiation process retinal stratification with typical apical basal polarity and correct spatial positioning of neural retina and RPE.[114] This process was soon replicated in human pluripotent stem cells[115] with multiple groups using a combination of 2- and 3D protocols to generate human retinal organoids comprising all the retinal cell types within a laminated fashion that resembles the adult retina and responds to light.[116-119] Identification of each retinal cell type within these complex organoids is not easy as cell-type surface markers are far and few in between. However, gene-editing techniques have been used by Collin et al.,[120] and in 2017 Phillips et al.[121] to generate reporter cell lines harboring the green fluorescent protein (GFP) or td-Tomato at the Cone-Rod Homeobox (CRX) gene locus, a key transcriptional factor in retinal development. In both cases, scRNA-seq was performed allowing assessment of retinal organoid culture. Collin et al. enriched the CRX+ve precursors and used an unsupervised clustering approach, SC3,[122] which combines different clustering outcomes into a consensus matrix and displays each cell and differentially expressed genes within clusters. This resulted in the identification of 2 clusters with cluster 1 comprising 72% of the cells and showing gene expression akin to early cone photoreceptors and cluster 2 showing a transcriptional profile associated with an “early photoreceptor precursor”-like state. Importantly, when the authors transplanted the CRX+ve in an animal model of retinal degeneration, the transplanted cells showed a cone phenotype in vivo, corroborating the scRNA-Seq study and indicating the power of these technologies to predict cell fate determination prior to lengthy and costly in vivo studies. In the study performed by Phillips et al., the traditional unbiased principal component analysis (PCA) and hierarchical clustering methods were used but they did not reveal clearly and defined cell populations. Instead a new method, which generates graphical clusters based on the relative expression levels of bait genes, which ranks the single cells based on the level of expression that most closely matches that of the bait was used, bypassing the requirement for cell sorting and enrichment of photoreceptors prior to the application of scRNA-Seq.

The possibility to sequence higher cell numbers (1976 cells) allowed Collin et al.[27] to resolve the complexity of retinal organoids in a dataset combining day 60, 90 and 200 of differentiation, identifying 9 clusters and revealing the presence of precursor and mature cell types including cone and rod photoreceptors, Muller glia and RGCs as well as RPE cells. The rarer cell types such as horizontal, amacrine and bipolar cells were not identified despite being detectable by immunohistochemistry, suggesting that a higher cell number need to be sequenced in order to detect the rarer cell types within these complex organoids. A similar study was also performed by Kim et al.[123] in 8-month old cone-rich retinal organoids. Single-cell sequencing analysis indicated that 60% of the cells were photoreceptors with a ratio of cones to rods of 1.4:1. Moreover, this analysis demonstrated similarity between cones in the organoids and human macula, indicating the power of scRNA-Seq technologies to decipher the cell type, transcriptome, and maturity between in vitro and in vivo retinogenesis.

In addition to photoreceptor and full organoids, single-cell sequencing technologies have also been applied to the analysis of RGCs, the projection neurons of the retina that transmit the visual information to midbrain through the optic nerve. Daniszewki et al.,[124] analysed 1,714 embryonic stem cell-derived RGC (enriched on the basis of reporter expression (Brn3b-mCherry)) and identified 3 distinct subpopulations of cells that included progenitors as well as mature RGCs, which contained genes related to the axon guidance together with the semaphoring interactions, extracellular matrix proteins and a downregulation of cell cycle genes. In a similar study, Langer et al.[125] were able to distinguish several RGC subtypes as well as molecular markers, which could be invaluable for improving and assessing the efficiency of differentiation protocols. It is known that each RGC subtype plays a different role in the visual pathway; in fact, molecular markers remain unknown for primates, including humans.[126] Langer et al. identified DCX as a new candidate gene for the direction-selective RGCs (DS-RGCs) and provided new insights on the abundance of each RGC subtype based on expression of RGC specific subtype markers. Taken together all these data point towards the utility of scRNA-Seq in deciphering the cell composition of pluripotent stem cell retinal organoids and retinal cells, discovering novel subtype-specific cell markers and assessing their similarity and maturation state in comparison to the equivalent cell types in the developing and adult retina. As the scRNA-Seq techniques improve and costs get lower, their application to drug repurposing/discovery and disease-modeling studies becomes more feasible and tangible with potential to deliver important insights for patients with retinal disease in a short time frame.

Conclusions and final remarks

Since the early work of Ramon y Cajal, depicting and drawing the various types of retinal cells and their interactions,[127] a lot of thought and effort has gone towards the retinal cellular classification by a combination of approaches including light and transmission electron microscopy, molecular biology, FISH, optical coherence tomography, etc. While distinctions at the higher levels of the cellular hierarchy, based on morphology, physiology, and gene expression, are relatively straightforward to make, it is more difficult to achieve the finer divisions of cellular subtypes with precision. This is where single-cell sequencing methods are making a huge and rapidly growing contribution to cell type specification of many human organs from development to adulthood as well generation of cell-type specific markers as summarized in this current review. Despite this remarkable progress, there remain unresolved questions. For example, how similar/different do the transcriptional profiles of the cells have to be in order to be associated/segregated into the same/different subtype? Are there profound differences in cell subtypes from the left and the right retina? How do we match the cell number and subtypes identified by different methods? How many methods are needed to validate the results of single-cell sequencing? From how many donors does one need to perform single-cell capture and analysis to achieve consistent and reliable results? What is the ideal cell number and sequencing depth for each tissue to achieve the best cellular and subcellular classification? How does single cell disassociation and time from retrieval affect the viability of different retinal cell types and the single-cell data? To date, the handful of studies performed in mouse, primate and human retina, as discussed in the main body of the review, have not managed to dissect all 100-cell types shown by microscopy and electrophysiological findings. This indicates that the cell number and sequencing depth reported in these initial studies are probably insufficient to reveal the complexity of adult retina. The success in classifying BCs[90] with shallow sequencing suggests that capturing more cells and sequencing at a lower read depth can be used for comprehensive cellular classification. This can be followed by resequencing of particular subsets, which are more dynamic and complex in their transcriptional profile, at a higher depth. In addition, advances made in the field of single-cell sequencing and continuously decreasing sequencing costs will permit single-cell studies in many labs with increased cell numbers. Combined analysis of datasets from multiple donors performed in different laboratories with appropriate computational methods that correct for batch sequencing effects, coupled with extensive validations, could provide a solid foundation for generating the cellular hierarchy and retinal cell-specific markers both during development and in adulthood.

Exciting developments in the field of spatial transcriptomics[128] and light sheet microscopy[129] have started to provide the necessary high spatiotemporal resolution. The complexity of these methods, the skills and equipment and the specificity of each tissue response to application of such techniques make it difficult to envisage how they could be commonplace practice in every laboratory. However, as seen with scRNA-Seq, rapid advancement, optimization, and commercialization (e.g., 10x Genomics) will lead to wider adoption of such techniques. Furthermore, strategic initiatives such as the Human Cell Atlas (HCA), which bring together excellent groups of scientists and bioinformaticians, will provide the technological and data-driven online accessible resource for all those interested in further studies. These could be directed at confirmation of gene expression at protein level and how this changes during development; others could compare the developmental and adult single-cell expression data to retinal organoids grown in their own laboratories to stage and assess the developmental maturity and cell type needed for transplantation or toxicology-based studies. Importantly, the wealth of data provided by single-cell sequencing together with regional and spatial maps provides an excellent resource for mapping the expression of known disease genes to certain cell types, which will facilitate better disease understanding as well as therapeutic interventions. Finally, yet importantly, identification of novel regulators during retinal development is going to provide the necessary baseline for identifying new disease genes, which has not been possible to date due to mutations in non-coding sequences or lack of understanding about the expression of the disease gene in particular retinal cell types. scRNA-Seq coupled with single-cell chromatin accessibility (for example ATAC-Seq) and bioinformatics analyses (pseudotime, etc.) are starting to identify the gene regulatory networks and their dynamics across species, cell surface markers and ligand-receptor pairs at a scale that was probably unimaginable during the previous decade. The question now is not when and at what scale, but how biologists, clinicians, and lawmakers, will best use this information to help treat patients with inherited and age-related retinal disease. Perhaps the experience set by the UK in exploiting the data generated from the 100,000 genomes could provide an enlightening example of how technological revolutions can bring about improvement in healthcare and life quality of many patients worldwide and bring the bench research closer to the bedside.