scRNA-Seq atlas of human developing eyes and retinas
To better understand the emergence of RPCs, their heterogeneity and differentiation to retinal cell types, we performed scRNA-Seq of 11 embryonic and fetal eyes spanning 7.5-15 PCW and 13 retinal samples from 10-21 PCW (Table S1). Each sample was embedded using Uniform Manifold Approximation and Projection (UMAP) and clustered using Seurat graph-based clustering (Figure S1, Table S1). In our earlier study of human retinal development 2, we reported expression of RPC markers (e.g. VSX2) at the peripheral margin of developing human retina as early as 6.5 PCW, and at the outer neuroblastic zone of central retina at 7.8 PCW, marking this time period an important developmental window of human RPC divergence. To better understand this process at the single cell level, we first focused our scRNA-Seq analysis in the 7.5-8.5 PCW embryonic eyes, revealing the presence of retinal and non-retinal cell clusters (Figure S1A). At this stage of human development, the neuroepithelium-derived optic vesicle is surrounded by periocular mesenchyme (POM), which is derived from neural crest and mesoderm and contributes to anterior and non-neural ocular tissue development. In accordance, we identified neural crest cell clusters adjacent to POM and corneal endothelial and corneal stromal cells 17 (Figure S1A). The POM cell clusters of 7.5-8.5 PCW eyes displayed high expression of characteristic markers described in other species such as collagen chains (COL3A1, COL5A1, COL1A1), proteoglycan decorin (DCN) 18 and latent Transforming Growth Factor Beta Binding Protein 1 (LTBP1) 19. Other non-retinal cell clusters including extraocular muscle, ocular surface epithelium, red blood cells, monocytes and macrophages were also identified in most eye samples of this developmental window (Figure S1A, Table S1).
Amongst the retinal cell clusters, RPCs and retinal ganglion cells (RGCs) were identified in all 7.5 - 8.5 PCW eyes, consistent with our previous findings of RGC presence at the basal side of the inner neuroblastic zone at 8 PCW 2. A recent scRNA-Seq study has documented the presence of early and late RPCs and their transcriptional signature 5. In accordance, we found that RPCs clusters of 7.5 and 7.7 PCW eyes displayed high expression of characteristic early RPC markers such as SFRP2, RAX, PAX6, VSX2, ZIC2, ZIC1, SIX3, SIX6 etc. 20 (Table S1). While the expression of these markers was maintained in RPCs found in the 8 and 8.5 PCW eye specimens, expression of late RPC markers such as CCND1, ASCL1 and HES6 5 became prominent, suggesting the co-existence of early and late RPCs. A proliferating cone photoreceptor cluster marked by high expression of proliferation (TOP2A, PCNA), neurogenic progenitors (HES6) and cone markers (RXRG, GNB3) was detected in two eye specimens obtained at 7.7 and 8 PCW (Figure S1A) as well as one of the 10 PCW samples (Figure S1B). In the latter sample, the proliferating cone cluster was distinct to the mature cones. Notably a small cluster of horizontal cells was detected in two specimens of 8 and 8.5 PCW of development (Figure S1A), corroborating with the first reported emergence of immature horizontal cells at day 59 of human development 4. A retinal pigment epithelium (RPE) cell cluster with high expression of characteristics markers namely PMEL, TYRP1, MITF and GJA1 was identified in all eyes of this stage (Figure S1A, Table S1).
A recent scRNA-Seq study provided evidence of the existence of neurogenic RPCs in addition to the early and late RPCs, based on gene expression profiling 5. Sridhar and colleagues 6 further defined the neurogenic into intermediate transitional T1, T2 and T3 cell populations with the capacity to give rise to defined types of retinal neurons at specific developmental windows. Our analysis of 10-14 PCW samples demonstrated the presence of some of these neurogenic transitional progenitors; however, a better definition was obtained from the analysis of retinal samples (> 10 PCW) due to a higher number of analysed cells within the retina per se (Figure S1B, Table S1). From 14 PCW, we were able to detect all types of retinal neurons apart from bipolar and Muller glia cells, which were detected from 14 and 16 PCW respectively (Figure S1B,C, Table S1). Published topographical studies report the highest density of microglia in the retinal periphery and retinal margin at ~ 8 PCW. Their presence in the central retina is not observed until 12 PCW (Diaz-Araya et al., 1995). In accordance with these studies, presence of cell clusters with high expression of microglia markers were observed from 14 PCW retinal samples (Figure S1B, C, Table S1).
Pseudo-temporal analyses reveal transition of RPCs to T1-T3 neurogenic cell clusters and retinal neurons
To identify gene regulatory networks (GRNs) that control RPCs specification and differentiation we integrated the transcriptomes of all retinal cells from 7.5-16 PCW human embryonic/fetal eyes with those of 10-21 PCW retinas. 48,856 cells were embedded using UMAP and clustered with Seurat (Figure 1A). Forty-three clusters were identified: of these, thirty-eight were composed of RPCs and retinal neurons, one of microglia, one of fibroblasts and two expressed markers of more than one cell type (Table S2). These last four clusters and an amacrine cell cluster with high expression of mitochondrial markers (cluster 41) were removed prior to further analysis. For illustrative purposes, all cell clusters defined as the same cell type, are annotated with the same color (Figure 1A).
Six RPC clusters were identified with enriched expression of typical markers including SFRP2, DAPL1, HES1, HMGA1, PAX6, RAX etc. Adjacent to the RPCs, a Muller glia cell cluster was found with high expression of RLBP1, SLC1A3, APOE, VIM, CLU, SOX9 (Table S2). Eight clusters named proliferating RPCs were also found next to RPCs and were characterized by high expression of proliferation markers (e.g., MKi67, TOP2A) and RPCs markers (SFRPP2, SOX2, HES1, ID3 etc.) (Figure 1A). All differentiated retinal cell types were easily identified based on the expression of characteristic cell markers (Table S2). Between the RPCs and retinal neurons, we identified a progenitor cell cluster with high expression of ATOH7, HES6 and DLL3, which was defined as the transitional T1 cluster. Two other cell clusters emanating from T1 were identified leading to the horizontal and amacrine cells clusters, and to photoreceptors and bipolar cells. These were defined as T2 and T3 based on the high expression of PTF1A and PRDM13 and FABP7 and OTX2 respectively (Table S2).
To fully define the transitions from RPCs to T1, T2 and T3 neurogenic progenitors, we performed pseudo-temporal analysis of gene expression changes following cell cycle regression (Figure 1B). This analysis revealed bimodal densities of RPCs, mirroring the transition between early and late RPCs described by Lu and colleagues 5, followed by the transitional T1 and then T2 and T3 neurogenic progenitors (Figure 1C, D). Early and late RPCs showed overall similar expression patterns of markers, but also some differences in the expression level. For example, the inhibitor of Wnt pathway (SFRP2) was more highly expressed in the early RPCs, while Muller glia cell markers (VIM, FOS) were more prominent in the late RPCs (Figure 1E, Table S2). High expression of well described neurogenic markers such as ATOH7, HES6 and OLIG2 was characteristic of the T1 and to a lesser extent of the T3 progenitor cluster. In accordance with the position of T2 between horizontal and amacrine cells, high expression of key transcription factors (ONECUT2, TFAP2C) required for differentiation to these two lineages was observed (Table S2). The transitional cluster T3 shared the expression of highly expressed T1 markers, and additionally displayed ethe expression of typical photoreceptor precursors (RXRG, NRL) and bipolar cell markers (VSX1), consistent with its positioning between the T1, and photoreceptor and bipolar cell clusters (Figure 1A, C, Table S2). To tease out the transcriptional machinery that controls the specification of each differentiated retinal cell type and their transition, separate pseudo-temporal analyses of RGCs, horizontal and amacrine, and photoreceptor and bipolar cells were conducted (Figures S2, Table S2). These analyses show that RGCs go through a T1 transition (data not shown), horizontal and amacrine cells through a T1-T2 (Figure S2A) and photoreceptor and bipolar cells through a T1-T3 transition state (Figure S2B), corroborating recently published scRNA-Seq data on few stages of human fetal retinal development 6.
During the course of these analyses many transcription factors highly expressed on defined progenitor subtypes or retinal neurons were identified. Human horizontal cells express the well-known markers such as ONECUT1, ONECUT3, LHX1, PROX1 and show considerable overlap in gene expression with amacrine cells which are characterised by high expression of TFAP2A, LHX9 and MEIS2 (Figure S2A, Table S2). The cone and rod precursors were characterised by high expression of THRB, and NRL and NR2E3 respectively, but interestingly also shared high expression of RGC and horizontal and amacrine cell markers (Figure S2B, Table S2). For example, high expression of RGC marker SNCG was found in the cone precursor cluster, while rod precursors displayed high expression of PROX1, a marker of retinal progenitors, horizontal and AII amacrine cells 21, 22. These findings were further corroborated by immunofluorescence analysis showing co-expression of cone photoreceptor marker RXRG with RGC marker SNCG in the 8 and 11 PCW, which tailed off at 12 PCW (Figure S3). Notably, these co-expressing cells were in a proliferative state at 8 PCW, but not in the later stages of development.
ST analyses reveal the location of early RPCs in the CMZ of 8-13 PCW developing human eyes
Our pseudo-temporal analysis above and a recent scRNA-seq analysis of human retinal development identified clear transcriptional signatures of early and late RPCs 5, however their spatial location in the developing retina has not been defined to date. To investigate this in detail, we performed ST analysis on an 8 PCW eye sample (Figure 2A-D) revealing the presence of 12 spatially organised cell clusters (Table S3) including POM, lens, vitreous, corneal stroma, extraocular muscle and RPE (Figure 2C). Notably, the optic stalk (cluster 9, Figure 2C) was clearly defined by the ST analysis and characterised by high expression of PAX2, ZIC1, HES1, SOX2, LHX2, THY1, PAX5 and ID3 (Table S3). Histologically, a single cell layer was present at the peripheral retina, but outer and inner neuroblastic layers were present in the central retina (Figure 2A, C). The cells residing in the inner neuroblastic layer of central retina (cluster 5) were characterised by high expression of RGC markers such as GAP43, PRPH, SNCG, INA, NEFL and NEFM, corroborating our previous immunofluorescence findings 2.
High expression of genes typically marking the ciliary margin (e.g., WNT2B 23), eye field (RAX and PAX6), earlyRPCs (ID3, HMGA1, MDK, SFRP2) and pigmented cell (PMEL, TYR and TYRP1) markers (Table S3) were characteristic of the cluster 4, which was defined as CMZ. Subclustering of cluster 4 resulted in two further cell subclusters (Figure 2E, F) with subcluster 0 expressing at high level marker genes associated with early RPCs (FGF19, SFRP2, DAPL1, ZIC1, ID1, ID3, HMGA1, EEF1A1, TPT1, TMSB4X), and subcluster 1 expressing at high level ciliary body (KCNJ8, TPM2, ADGRA2A) and iris pigment epithelium signature markers (TYR, TRYP1, DCT, SILV, MLANA, PMEL). In contrast, the progenitor cells residing in the outer neuroblastic layer of the central retina (cluster 10, Figure 2C) were characterised by high expression of late RPC markers such as RORB, CKB, HES1, HES5, ASCL1, HES6 and NEUROD1 (Table S3). This led us to hypothesise that the early and late RPCs could be spatially located in different regions of the developing human retina. To investigate this further, we generated a gene expression signature for early and late RPCs from our pseudo-temporal analysis. Violin plots (Figure 2G, I) indicated the highest expression of early RPCs markers in the subcluster 0 of CMZ (cluster 4), while the highest expression of the late RPCs markers was observed in the outer neuroblastic layer of central retina (cluster 10). These findings were further corroborated by the spatial mapping, showing the predominant localisation of early RPCs in the CMZ, the late RPCs in the central retina (Figure 2H, J), and the transitional neurogenic clusters T1, T2 and T3 predominantly in the outer neuroblastic layer of the central retina (cluster 10, Figure S4). Interestingly, an early RPCs expression signature was also observed in and around the optic stalk (Figure 2H), albeit less intense than in the CMZ.
Similar cell clusters with enriched expression of RPCs and pigmented epithelial cell markers were consistently found in the ciliary margin zone of 10, 11 and 13 PCW eye specimens analysed with the ST approach (Figure 3A-C, Table S3); but the frequency of early RPCs was very much reduced compared to the 8 PCW sample. We were unable to fit larger size eyes to the current Visium ST slides, thus it was not possible to extend the ST analyses to fetal eyes older than 13 PCW. However, we utilised the scRNA-Seq data and plotted the ratio of early to late RPCs, showing a significant reduction in fraction of early RPCs from 8 PCW onwards (Figure 3D).
To complement the ST analyses and obtain more insights into early and late RPC localisation into the CMZ during retinal development, we performed RNA-Scope investigations using a marker of early RPCs (ZIC1), iris pigmented epithelium (TFPI2), late RPCs/neurogenic progenitors marker (HES6) and ciliary body (OPTC) in eye samples from 6.3-16 PCW. Consistently with data obtained above, we observed ZIC1 expression in the CMZ as well as rest of retinal neuroepithelium of 6.6 PCW retinas (Figure S5A, B). In contrast, the expression of HES6 was first seen in the central retina (Figure S5B) spreading to the periphery, up to CMZ, from 8 PCW until the last 16 PCW specimens examined (Figure S5D, G, Figure S6C, D, G). Although ZIC1 expression was still present in the CMZ of 8 PCW retinas (Figures S5C, D) it was reduced and by 10 PCW, it was localised in small patches at the posterior end of CMZ of less than 50% of eye sections assessed (Figure S5E, F, G). In 13 and 16 PCW, ZIC1 expression was completely absent from the CMZ (Figure S6A-F), however its expression in the rest of retinal neuroepithelium persisted in a co-expressing pattern with HES6 (Figure S6C, D, G). Together these data corroborate localisation of early RPCs in the CMZ of developing human retinas with decreasing frequency from 8 PCW of development and reveal the propagation of neural retina differentiation from the centre to the periphery.
Characterisation of chromatin accessibility during human retinal development
To investigate chromatin accessibility during human retinal development, we performed scATAC-Seq analyses of two developing human eyes (8.5 PCW), and ten retinal samples dissected from 10-21 PCW (Table S4). 118,883 cells were captured using the 10XGenomics Chromium Single Cell ATAC Library and Gel Bead Kit. The corresponding scRNA-Seq datasets were used as reference maps to identify ATAC-Seq clusters for each sample (Figure S7). Similarly, to scRNA-Seq analysis, defined clusters of corneal stroma, epithelium, endothelium, and periocular mesenchyme, CMZ, extraocular muscle, red blood cells, microglia and optic nerve were found in the 8 PCW human eye samples alongside retinal clusters comprised of RPCs, RGCs, horizontal and amacrine cells and transient neurogenic cluster T1 (Figure S7A). All the retinal cell types were present between 10-14 PCW, albeit some at the precursor stage (for example rod precursors at 10 PCW), in addition to RPCs and the three types of transient neurogenic clusters T1, T2, T3 (Figure S7A). From 16 PCW, small clusters of microglia and the last-born cell type, Muller glia cells were evident (Figure S7B), corroborating the scRNA-Seq data and previously published sequential order of retinal cell emergence 2, 5, 6.
54,045 retinal cells from 12 samples and 9 developmental timepoints encompassing 8-21 PCW with a 3694,28 average median fragments per cell were integrated. The cells were clustered using the chromatin accessibility peaks near previously known marker genes resulting in 22 clusters (Figure S7C, Table S4). These included the abundant cell types (for example rods, cones) as well as the rarer cell types (e.g., microglia) and several types of amacrine cells (gabaergic, glycinergic and starburst). The DNA accessibility peaks were classified using annotation from cellranger and associated with promoters (if found within -1000 to +100 bp of the transcription start sites), exons, introns, distal (if found within 200 kb of the closest transcription start site), or intergenic regions (if not mapped to any genes) (Figure 4A, Table S5). This analysis enabled us to identify cell type specific regions of accessibility in RPCs, transient neurogenic progenitors T1, T2, T3 and retinal neurons (Figure 4B). All clusters displayed scATAC-Seq marker peak enrichment of cell type specific marker genes (Figure 4C).
We then went on to predict transcription factor (TF) binding motifs within the scATAC peaks using Signac (Figure 5A, Table S6), followed by foot printing validation analysis (Figure 5B, C). In RPCs, we identified binding motifs for TFs expressed in the optic stalk (VAX1, VAX224), eye field (RAX25) and RPCs (LHX6, VSX2, SOX226, 27, 28) as well as TF binding motifs that were shared with other cell types (for example FOS, NFI family members and SOX6 binding motifs were shared with Muller glia cells, Figure 5A). Notably we identified binding motifs for TFs not previously associated with RPCs, for example EN1, GBX2, LBX2, SHOX2 and GSX1 (Figure 5A,B, Table S6).
In accordance with sequential emergence of T2 and T3 neurogenic clusters from T1, we observed a set of shared TF binding motifs for NEUROG2, NEUROD1, HAND2, TAL1:TCF3, PITF1A, ATOH7, BHLHA15, MYF5 and ATOH1 (Figures 5A). The earliest born retinal cell types, RGCs displayed the POU4F and EBF family members binding motifs, whilst horizontal and amacrine cells were characterised by ONECUT and CUX, and TFAP2 and MEIS family members binding motifs respectively (Figure 5A, C). Rod and cone photoreceptors are derived from the T3 neurogenic progenitors; hence shared binding motifs were identified for TFs such as OTX2 29, CRX 30 and DMBX1 31 (Figure 5A, C), which are well described in the literature for their role in photoreceptor specification. Importantly, shared binding motifs in T3 and photoreceptors were also discovered for PITX1 and GSC TFs not associated previously with a role in rods or cones (Figure 5A). Bipolar cells are also derived from the T3 progenitors; hence enrichment of T3 TF binding motifs (such as OTX2) was observed in addition to less well characterised TFs such ZBTB18 (Figure 5C). Enrichment of NFI binding motif family members was observed for glia cells (Muller glia and microglia) in accordance with their role in regulating specification of the late-born cell types in the retina 32.
scATAC-Seq enables prediction of gene regulatory networks and novel TFs that govern retinal development
To identify GNRs governing the RPCs differentiation to transient neurogenic progenitors (T1, T2, T3) and the retinal neurons, the upstream regulator tool in IPA was used, combined with overlay analysis of DA peaks. A large number of upstream regulators including TGFβ, IGF-1, FGF-2, Sonic hedgehog were predicted to be activated in RPCs (Table S7). Their main common characteristic was predicted activation of key genes encoding proteins important for RPCs proliferation (HES1 33, CCND1 34, ID3 35) and inhibition of transcription factors (ATOH7 36, NEUROD4, NEUROD1 37, PTF1A 38, and HES6 39) (Figures 6A, B) that define the RPC competence to RGCs, horizontal and amacrine cell fates, and photoreceptor differentiation. Amongst the predicted inhibited upstream regulators, we found key transcription factors such as PAX6 40and ASCL1 41, shown to control the timing and specificity of retinal neurogenesis (Figures 6C, D). Together these findings suggest the presence of a finely tuned balance between activation of proliferation regulators and inhibition of neurogenic cues to maintain RPC self-renewal. Conversely in the transient neurogenic progenitors, we predicted the activation of upstream neurogenic regulators (e.g., ASCL1 in T1, FOXA2 in T2, GTF2IRD1 in T3), which control the expression of TF necessary for retinal cell type specific differentiation (Figure S8A-C) and inhibition of regulators (e.g. BMP4, LIN28A, IGF-1) that govern RPCs fate (Figure S8D-F).
Performing the same analysis on the differentiated retinal cells revealed some interesting insights and predicted novel upstream regulators (Table S7). For example, Eomesodermin (EOMES), a target gene of Pou4f2, required for RGC and optic nerve development in mouse 42, was predicted to be activated in horizontal cells (Figure S9B), resulting in activation of LHX1 26, 43, a transcription factor that specifies horizontal cells, while suppressing LHX9, a transcription factor required for amacrine cell subtype specification 26. We identified a novel putative upstream regulator in RGCs, namely KLF2 (Figure S9A), which is predicted to inactivate Notch signalling, an important event required for RGC differentiation 44. Our scRNA-Seq data have shown that the basic-helix-loop-helix PTF1A is highly expressed in the T2 progenitors, which give rise to amacrine and horizontal cells. A similar role for PTF1A has been demonstrated during mouse retinal development 38, albeit transient expression in all types of amacrine cells has been demonstrated in the zebrafish 45. In accordance with these findings, the PTF1A upstream regulator was identified in amacrine cells, resulting in activation of TFs that promote amacrine cell specification and function (Figure S9C).
Our analysis was also able to predict more complex upstream regulators, which are ligand-activated and able to play a role in the control of gene expression in a cell type specific manner. For example, the peroxisome proliferator-activated receptor γ (PPARG), was identified as a putative novel upstream regulator in cone photoreceptors (Figure S9D) and predicted to activate amongst others PRDM1, which has been demonstrated to stabilise photoreceptor cell fate in OTX2+ progenitors by preventing bipolar cell induction 46. In contrast, the upstream regulator TCF7 in bipolar cells is predicted to suppress PRDM1 and activate expression of key genes important for bipolar cell function such as GNAO1 and TRPM1 47 (Figure S9F). G-protein coupled receptors (GPCRs) play a significant role in many tissues by transducing complex signalling networks that coordinate gene expression. In accordance, one of the putative activated upstream regulators in rods, was the GPCR Rhodopsin (RHO), the most abundant protein in rods which functions as the primary photoreceptor molecule of vision (Figure S9E). Our data suggests that RHO may regulate the transcription of rod specific phosphodiesterases (e.g., PDE6G, PDE6B) and transducins (GNAT1, GNGT1).
scATAC-seq trajectory analysis reveals a role for TEAD transcription factors in retinal lamination and RPC, photoreceptor and RGCs specification
The TF binding motif analysis in RPCs (Figure 5A) revealed an enrichment for TEA domain (TEAD) TFs in RPCs (Figure 7A, Table S6). TEADs play an essential role in mediating YAP-dependent gene expression resulting in transcription of target genes responsible for cell proliferation, inhibition of apoptosis or retinal neurogenesis 48. During retinal development, YAP’s expression is restricted to the outer neuroblastic layer where RPCs reside 49. The TEAD TFs display different expression patterns in the retina, with Tead2 being highly expressed in the proliferating cells located at the basal side of the outer neuroblastic layer of mouse retina and Tead3 being highly enriched in the inner neuroblastic layer (genepaint.org).
To assess if TEADs binding is essential for RPC development and differentiation, we treated PSCs-derived retinal organoids with TEAD palmitoylation inhibitor MGH-CP1, that blocks TEAD2/4-YAP interaction and suppresses the expression of their target genes in cancer cell lines 50. Treatment was started at day 25, coinciding with immunofluorescence detection of RPCs in the retinal organoids. Three different doses (2.5, 5 and 10 µM) of MGH-CP1 were added every two days to the culture media for 14 days. Under control conditions (treated with vehicle) retinal organoids develop a bright phase retinal neuroepithelium by day 28 (day 3 of treatment) which expanded over time. By the end of treatment the majority of control retinal organoids (> 75%) had full coverage with bright phase retinal neuroepithelium, with a minority of organoids displayed partial coverage (15%) or no coverage (7%) at all (Figure S10A, B). In contrast, the great majority of retinal organoids treated with MGH-CP1, displayed either partial or full loss of bright phase neuroepithelium in a dose dependent manner. These findings were fully corroborated by immunofluorescence analyses, which showed the presence of partial retinal neuroepithelium harbouring VSX2+ RPCs in organoids treated with 5 µM MGH-CP1 and much reduced and mislocalised RPCs in the organoids treated with 10 µM MGH-CP1 (Figure 7B, F). The latter were also often characterised by the presence of SNCG+ RGCs in the apical layer of the organoids and internal rosettes with VSX2+ RPCs or Ki67 proliferating cells (Figure 7C-E). Quantitative analysis demonstrated attenuated RPCs, photoreceptor, and RGCs specification in a dose-dependent manner in retinal organoids treated with MGH-CP1 (Figure 7f, H-J). In accordance with TEADs role in activating transcription of genes important for cell proliferation, we noted a significant reduction in the fraction of Ki67+ cells in the retinal organoids treated with the highest dose of MGH-CP1 (Figure 7G), which also showed a small but significant increase in the percentage of apoptotic cells (Figure S10C).