Single Cell Transcriptomic Analyses Reveal the Impact of bHLH Factors on Human Retinal Organoid Development

The developing retina expresses multiple bHLH transcription factors. Their precise functions and interactions in uncommitted retinal progenitors remain to be fully elucidated. Here, we investigate the roles of bHLH factors ATOH7 and Neurog2 in human ES cell-derived retinal organoids. Single cell transcriptome analyses identify three states of proliferating retinal progenitors: pre-neurogenic, neurogenic, and cell cycle-exiting progenitors. Each shows different expression profile of bHLH factors. The cell cycle-exiting progenitors feed into a postmitotic heterozygous neuroblast pool that gives rise to early born neuronal lineages. Elevating ATOH7 or Neurog2 expression accelerates the transition from the pre-neurogenic to the neurogenic state, and expands the exiting progenitor and neuroblast populations. In addition, ATOH7 and Neurog2 significantly, yet differentially, enhance retinal ganglion cell and cone photoreceptor production. Moreover, single cell transcriptome analyses reveal that ATOH7 and Neurog2 each assert positive autoregulation, and both suppress key bHLH factors associated with the pre-neurogenic and states and elevate bHLH factors expressed by exiting progenitors and differentiating neuroblasts. This study thus provides novel insight regarding how ATOH7 and Neurog2 impact human retinal progenitor behaviors and neuroblast fate choices.


INTRODUCTION
As an integral component of the central nervous system, the vertebrate neural retina retains a highly conserved laminar structure that senses, processes, and delivers visual information to the brain. Classic cell birth dating and lineage tracing studies have established that the seven major neuronal cell types constituting the retinal network are generated in a temporal order from a common ocular progenitor pool during development (Young, 1985;Livesey and Cepko, 2001;Wong and Rapaport, 2009). The ensuing research has ruled out a rigid deterministic cell fate specification mechanism, but instead supports a view that multipotent progenitors progressively evolve through different competence states to enable the sequential production of distinct cell types (Cepko et al., 1996;Livesey and Cepko, 2001). Cumulative molecular genetic studies have uncovered important roles of cell intrinsic factors involved in retinal development (Xiang, 2013). Among these, transcription factors containing the basic helix-loop-helix (bHLH) motif have emerged as important players regulating the production and differentiation of various retinal cell types (Akagi et al., 2004). Multiple bHLH factors are expressed during retinal development either in proliferating progenitors or in postmitotic neurons; however, their dynamic regulation and function in specific cellular contexts remain to be fully elucidated.
The bHLH factors Atoh7 and Neurog2 are both expressed early in the developing vertebrate retinal epithelium. Atoh7 plays a critical role in the development of an early born retinal neuronal type, the retinal ganglion cells (RGCs), which project axons through the optic nerve to multiple higher visual centers (Young, 1985;Hoon et al., 2014). Both Atoh7 and Neurog2 mRNAs are expressed by subsets of early retinal progenitors, with some coexpression at the protein level during certain time windows (Brown et al., 1998;Fu et al., 2009;Miesfeld et al., 2018a). Loss of Atoh7 function in mouse, zebrafish, and humans results in a severe reduction of RGCs leading to a diminished optic nerve and blindness (Brown et al., 2001;Kay et al., 2001;Wang et al., 2001;Ghiasvand et al., 2011;Miesfeld et al., 2020). Atoh7 deficiency also causes a minor abnormality in the production of cone photoreceptors, another early-born cell type in the retina (Brown et al., 2001). In contrast, genetic ablation of Neurog2 yields a transient stall of neurogenesis but without severe lasting deficits (Hufnagel et al., 2010). In the mouse retina, Atoh7 protein is not detected in fully differentiated RGCs (Fu et al., 2009;Miesfeld et al., 2018a), suggesting that its main biological activity is transiently required in uncommitted early progenitors. Ectopic expression of Atoh7 in different late stage retinal progenitors either redirects progenitors toward an RGC fate (Mao et al., 2013) or fails to specify the RGC fate . Therefore, Atoh7 is thought to confer a competent state of progenitors to adopt early cell fates (Brzezinski et al., 2012). In the absence of Atoh7, co-expression of two downstream transcription factors Islet1 and Pou4f2 is sufficient to rescue the RGC production deficit and ensure full execution of the RGC differentiation program in the mouse retina (Liu et al., 2001;Mu et al., 2008;Pan et al., 2008;Li et al., 2014;Wu et al., 2015). We have shown previously that viral mediated expression of human ATOH7 in the developing chicken retina results in precocious neurogenesis and a significant increase in RGC production (Zhang et al., 2018), supporting a hypothesis that a critically high threshold of ATOH7 expression triggers uncommitted early progenitors to exit the cell cycle and predominantly adopt the RGC fate .
In the human retina, RGC development occurs during the first trimester and remains a minor cell population in the mature retina (Hendrickson, 2016;Hoshino et al., 2017;Mellough et al., 2019). This scarcity of human RGCs has hindered research on human RGC development as well as blinding diseases caused by RGC loss. The advancements of pluripotent stem cell technologies in the preceding decade have led to robust stem cell based retinal organoid culture systems (Meyer et al., 2011;Nakano et al., 2012;Reichman et al., 2014;Ohlemacher et al., 2016), thus providing an excellent opportunity to produce and study human RGC development in vitro. Recent single cell transcriptomic analyses have revealed that primate retinas, including humans, show distinct molecular features and RGC subtype proportions compared to rodents (Liang et al., 2019;Lo Giudice et al., 2019;Lukowski et al., 2019;Menon et al., 2019;Peng et al., 2019;Cowan et al., 2020;Hoang et al., 2020;Lu et al., 2020;Sridhar et al., 2020;Yan et al., 2020). In this study, we have examined the function of bHLH factors ATOH7 and Neurog2 on human RGC development using embryonic stem cell (ESC)-derived 3D retinal organoids. Our results demonstrate that elevating these two factors in uncommitted human retinal progenitors asserts powerful neurogenic effects. By performing single cell transcriptomic analysis, we identify distinct statuses of human retinal progenitors, including a population poised to exit the cell cycle. Our results show that ATOH7 and Neurog2 participate and regulate an interactive gene network to accelerate progenitors through two transitional stages to adopt postmitotic neuronal identities.

Viral Mediated bHLH Factor Expression Affects Progenitor Proliferation in 3D Human Retinal Organoids
To investigate the roles of bHLH neurogenic factors during development of the human retina, we established H9 embryonic stem cell (ESC)-derived 3D retinal organoid cultures (Kuwahara et al., 2015;Ohlemacher et al., 2016). These human retinal organoids showed typical morphology of the retinal neural epithelium and co-expressed PAX6 and VSX2 transcription factors (Supplementary Figure 1), a characteristic feature specific to retinal progenitors. To regulate gene expression during retinogenesis, we constructed lentiviral vectors that encode the doxycycline (Dox) inducible TetO promoter upstream of the human ATOH7 cDNA fused with the Flag epitope tag (LV-ATOH7f) or co-expressing EGFP and puromycin resistant genes (LV-AEP) ( Figure 1A). We also produced a previously described lentiviral TetO vector co-expressing the mouse Neurog2, EGFP, and puromycin resistant genes (LV-NEP) (Zhang et al., 2013; Figure 1A). Immunohistochemistry confirmed Dox-induced bHLH protein expression following in vitro co-infections of H9 ESCs with LV-rtTA and either LV-ATOH7f, LV-AEP, or LV-NEP ( Figure 1B). We next carried out co-infection of retinal organoids with LV-rtTA and either LV-AEP or LV-NEP at the onset of retinogenesis followed by Dox inductions (Figure 1C). Live imaging of the EGFP reporter showed that at 24 hours post Dox induction, significant numbers of LV-AEP or LV-NEP infected cells had already migrated toward the inner retina compared to cells from the control LV-GFP infected organoids. This trend continued and became more obvious as the Dox induction times lengthened (Supplementary Figure 2). By 6-days after Dox induction, the majority of LV-AEP and LV-NEP infected cells were located in the inner layer of the organoids (Figures 1D,E), indicating that viral mediated ATOH7 and Neurog2 expression impacted retinal organoid development.
To investigate whether virally expressed bHLH factors affected retinal progenitor proliferation, we performed BrdU pulse-labeling and examined progenitor marker expression by immunohistochemistry. In the control LV-GFP infected retinal organoids, BrdU-labeled and the proliferating cell nuclear antigen (PCNA)-labeled progenitors occupied the ventricular zone, whereas phospho-histone 3 (PH3) labeled M phase cells were detected at the ventricular surface of the organoids (Figures 2A-C). Furthermore, PCNA showed co-labeling with PAX6-expressing progenitors in the ventricular zone, but was absent from the PAX6 + cells located in the inner layer, where postmitotic POU4F + RGCs resided (Figures 2B,C). In contrast to LV-GFP infected retinal organoids, in which GFP + cells were distributed throughout the neural epithelium, the majority of GFP + cells in LV-AEP and LV-NEP infected organoids were located in the inner retinal layer and devoid of co-labeling with BrdU or PCNA (Figures 2D,E). Quantification of dissociated FIGURE 2 | Influences of viral mediated neurogenic factor expression on cell proliferation. (A-C) Immunofluorescence labeling of cross sections from retinal organoids infected with LV-GFP. At day 34, PH3 + dividing cells are located at the ventricular surface (A), while PCNA + progenitors occupy the ventricular zone and are co-labeled with PAX6-expressing progenitors (B). Note that PAX6 + PCNApostmitotic neurons are located to the inner layer of the retinal organoid. At day 40, proliferating progenitors labeled with BrdU are distributed in the ventricular zone, whereas POU4F + postmitotic retinal ganglion cells (RGCs) reside in the inner layer of the retinal organoid (C). Scale bars, 20 µm. (D,E) Confocal images of cross sections from retinal organoids infected by LV-GFP, LV-AEP, or LV-NEP after 6-day Dox induction. Co-staining of GFP with BrdU at day 36 (D, 3-h labeling) or PCNA at day 40 (E) shows that most LV-AEP and LV-NEP infected cells are postmitotic. Scale bars, 20 µm. (F) Quantification of double labeled PCNA + GFP + progenitor cells among total viral infected GFP + cells at day 47. Mean ± S.E.M. of n = 3 independent samples. One-way ANOVA, **p < 0.01. retinal organoids confirmed that percentages of GFP + PCNA + double labeled cells were reduced significantly from 52.7 ± 8.0% for LV-GFP infection to 9.0 ± 3.8% and 2.0 ± 1.9% for LV-AEP and LV-NEP infections, respectively ( Figure 2F). These results demonstrate that viral driven ATOH7 or Neurog2 expression promoted cell cycle exit.

Elevated Neurogenic Factor Expression Promotes RGC Production in 3D Retinal Organoids
The retinal projection neurons (RGCs) are among the earliest neuronal cells produced during retinogenesis (Young, 1985;Hoshino et al., 2017). In our human retinal organoid cultures, RGC genesis was detected as early as day 25 and continued through day 60 ( Figure 2C). By day 40, retinal organoid-derived neurons exhibited voltage-gated Na + , K + , and Ca 2+ channels as well as spontaneous and provoked electrophysiological excitability in vitro during whole cell patch clamp recording, characteristic of native RGCs (Supplementary Figure 3). TTXsensitive Na + channels, TEA-sensitive K + channels, and Cd 2+sensitive Ca 2+ channel currents were recorded. Average Na + current amplitude was 219 ± 76 pA (n = 12) but in cells with large Na + currents (∼1 nA), multiple action potentials were observed, while cells expressing smaller Na + currents (200-400 pA) typically produced a single spike. Average K + current amplitude at + 40 mV was 438 ± 49 pA (n = 40).
To determine whether viral mediated ATOH7 or Neurog2 expression promoted RGC production, we performed immunohistochemical analyses of retinal organoid sections using known RGC markers (Figure 3). In LV-ATOH7f infected retinal organoids, signals of the viral reporter Flag closely correlated with POU4F-expressing RGCs in the inner retina as detected by a pan-POU4F/BRN3 antibody ( Figure 3A). Compared to control LV-GFP infected organoids, both LV-AEP and LV-NEP infected cells showed increased co-labeling with the RGC markers NF145, NeuN, and DCX in the inner retina ( Figure 3B). Similarly, in attached organoid cultures that displayed extensive neurite outgrowth, the ATOH7f-expressing cells showed extensive co-labeling for the RGC markers POU4F, NF145, and RBPMS ( Figure 4A).

Single Cell RNA-Sequencing Analysis Reveals Accelerated Neurogenesis
To elucidate the influences asserted by viral driven ATOH7 and Neurog2 expression on transcriptome, we performed single cell RNA-sequencing (sc RNA-seq) analysis. We first used fluorescent activated cell sorting to enrich for LV-GFP, LV-AEP, and LV-NEP infected retinal organoid cells between days 45 and 48 (Supplementary Figure 4), followed by 10X Genomics automated single-cell capture, mRNA barcoding, and cDNA library preparation. The high throughput DNA sequencing resulted in 192,932, 153,154, and 132,651 mean readings per cell for LV-GFP, LV-AEP, and LV-NEP samples, respectively. After aligning to the reference human genome and eliminating poor quality cells, the final single cell datasets had a mean gene range between 2,935 and 3,079 per cell, and consisted of 3,004 cells for LV-GFP, 2063 cells for LV-AEP, and 3909 cells for LV-NEP infected samples.
The sc RNA-seq datasets were subjected to Seurat cell clustering analysis (Butler et al., 2018), which resolved into 12-15 clusters, and visualized as 2-dimensional UMAPs ( Figure 5A). We applied dot plot analysis using known genes to assign various cell clusters into seven categories or states; each was represented with one or two highly or uniquely expressed genes ( Figure 5B and Supplementary Figure 5). The pre-neurogenic progenitor (PNP), neurogenic progenitor (NP), and cell cycle-exiting progenitor (EP) categories all expressed CCDN1, encoding cyclin D1, indicating that they all belonged to proliferative cells. However, the pre-neurogenic progenitor cluster cells showed distinctively high levels of PLOD2 transcripts (procollagen-lysine, 2-oxoglutarate 5-dioxygenase 2), whereas the exiting progenitor cells exhibited upregulation of GADD45A (growth arrest and DNA damage inducible 45 alpha) and downregulation of the retinal progenitor marker VSX2. The exiting progenitor cells and postmitotic neuroblasts (NB) were the main populations expressing ATOH7 mRNA (see Figure 9). Markers identifying the early born retinal neuronal types included ISL1 and POU4F2 for RGCs, CRX, and RCVRN for photoreceptors, and PRDM13 and TFAP2A for horizontal and amacrine cells (HC/AC) ( Figure 5B and Supplementary Figure 5). These sc RNA-seq analyses thus not only allowed identification of the early postmitotic neuronal types in the human retinal organoids at the time of analysis, but also revealed sub-classes of neural progenitors with distinct molecular signatures, and transitional states including cell cycle exiting progenitors and postmitotic neuroblasts.
Next, we combined different cell clusters assigned to each of the seven cell categories ( Figure 6A) and performed a pseudotime trajectory analysis ( Figure 6B). Without defining the starting and ending points, the trajectory of LV-GFP infected cells revealed that the pre-neurogenic progenitors were closely related to the neurogenic progenitors, which in turn produced exiting progenitors that developed into the postmitotic neuroblasts and a single trajectory including three types of retinal neurons. In both LV-AEP and LV-NEP infected samples, the pseudotime trajectory For (C-E) Mean ± S.E.M. of n = 3 independent samples. One-way ANOVA, ****p < 0.0001, ***p < 0.001, **p < 0.01. "*" is used to represent levels or ranges of probability, i.e. statistical significance.
displayed a more clearly defined progression from the neurogenic progenitors toward the existing progenitors, which in turn gave rise to neuroblasts. The neuroblast cells showed a single node for bifurcated trajectories separating the RGC and HC/AC branch from the photoreceptor branch ( Figure 6B). Quantification of the different cell categories/states further demonstrated that LV-AEP and LV-NEP infection significantly reduced the percentage cells in the pre-neurogenic progenitor category from 25.2% to less than 4% ( Figure 6C). Neurog2 expression also caused a significant reduction of neurogenic progenitors from 45 to 26.9%. Concomitantly, ATOH7f and Neurog2 induction increased the percentage of exiting progenitors from 2.2 to 4.5% as well as the percentage of neuroblasts from 5.0 to 14.9% and 20.2%, respectively. Moreover, consistent with marker analysis, elevated ATOH7f and Neurog2 expression resulted in increased proportion of RGCs from 7.0 to 22.1% and 27.7% of total cells, respectively ( Figure 6C). Interestingly, Neurog2, but not ATOH7f, enhanced the photoreceptor population from 3.4 to 9.8%; whereas ATOH7f, but not Neurog2, promoted HC/AC proportions from 3.1 to 4.6% ( Figure 6C). These sc RNA-seq data demonstrated that elevating ATOH7f and Neurog2 expression promoted transitions from the pre-neurogenic to the neurogenic state and enhanced neurogenesis of early retinal cell types in retinal organoids.

Neurogenic Factors Promote Transitions of Distinct Developmental States
Since viral mediated expression of ATOH7 and Neurog2 affected transitions between developmental states, we explored the key characteristics of each cell state. We compiled differentially expressed genes (DEGs, adjusted p < 0.05) of each cell category for LV-GFP, LV-AEP, and LV-NEP samples (Supplement Tables 1-3), and constructed heatmaps for the top 10 DEGs displaying more than Log 1.5-fold change of expression levels ( Figure 7A and Supplementary Tables 4-6). In the control LV-GFP infected sample, the top 10 DEGs in the preneurogenic progenitor category included SLC2A1, encoding glucose transporter protein type 1 (GLUT1), and GPI, glucose-6 phosphate isomerase, and were quite distinct from those of the neurogenic progenitors. In LV-AEP and LV-NEP infected retinal organoids, the pre-neurogenic progenitors were not only reduced, but also shared genes, such as SFRP2, IFITM3, and VIM with the neurogenic progenitors. In all three virus infected samples, the exiting progenitors shared top DEGs, including HES6, a suppressor of HES1 (Bae et al., 2000), and HMGB2, a member of the chromosomal high mobility group of proteins (Bianchi and Agresti, 2005), indicating common processes involved in cell cycle withdrawal. The heatmaps also revealed that neuroblasts in each virus infected sample expressed some genes associated with the RGC, photoreceptor, and HC/AC cell lineages, suggesting that the newly postmitotic neuroblasts were poised at an intermediate developmental state and in the process of committing to specific cell fates. To further decipher the changes of gene expression and the key biological processes that occur during normal developmental transitions, we constructed volcano plots ( Figure 7B) and performed gene ontology (GO) analyses ( Figure 7C) using significant DEGs from the LV-GFP infected cell categories. The pre-neurogenic progenitor state was associated with higher levels of transcripts for several glycolytic pathway enzymes such as SLC2A1, GPI, ENO1, and PGK1, as well as low levels of mitochondrial respiratory chain components and proteins required for rapid cell proliferation such as PCNA and TOP2A ( Figure 7B). This was consistent with the dominant GO pathways for the pre-neurogenic state ( Figure 7C). In contrast, cells in the neurogenic progenitor state expressed high levels of SFRP2, FGF19, and SOX2 (Figure 7B), correlating with the biological responses to growth stimulation and the mitotic cell cycle ( Figure 7C). The dominant GO terms associated with the exiting progenitors included cell cycle arrest and Notch signaling ( Figure 7C). The exiting progenitor to neuroblast transition was marked by elevated expression of the Notch signaling inhibitor HES6 and the antiproliferation factors BTG1 (Liu et al., 2015), BTG2, and GADD45A (Barreto et al., 2007; Figure 7B), consistent with the GO pathway analysis ( Figure 7C). The predominant GO pathways for the three neuronal types reflected their phenotypic differentiation with high axonal growth and light transduction processes associated with RGCs and photoreceptors, respectively ( Figure 7C).

ATOH7 and Neurog2 Assert Differential Effects on Cell Cycle and Neuronal Differentiation
To examine the potential impact of neurogenic factor expression on the cell cycle, we combined the LV-GFP, LV-AEP, and LV-NEP sc-RNA-seq datasets and carried out cell clustering analysis ( Figure 8A). Based on feature plots of known genes (Supplementary Figure 6), cell clusters in the combined UMAP were identified as distinct cell categories. Furthermore, the progenitor cell clusters were assigned to different phases of the cell cycle ( Figure 8A) based on known functions of cell cycle genes, including CDK2 associated with G1/S phase (Horiuchi et al., 2012), MCM4 as a DNA replication licensing factor (Yabuta et al., 2003), CCNB2 encoding M phase cyclin B2 (Liu et al., 1999), and PLK1 involved in spindle assembly and   profiles since cluster 6 was poised to exit the cell cycle ( Figure 8A). As expected, individual UMAPs of ATOH7 and Neurog2 virus infected samples displayed significantly reduced pre-neurogenic progenitor populations (Figures 6C, 8B). In addition, quantitative analyses showed that ATOH7 and Neurog2 virus infection altered cell cycle distributions among progenitors compared to controls (Figure 8C). For example, LV-AEP increased G1 distribution by 5% while decreasing S phase cells by 3.7%, whereas LV-NEP caused a 7% S phase cell reduction and a 4.5% increase of M phase cells, which included the exiting progenitors ( Figure 8C).
The Slingshot analysis demonstrated that the postmitotic neuroblast pool, cluster 7, served as the root source giving rise to three neuronal cell lineages ( Figure 8A). Interestingly, UMAPs of individual virus infected retinal organoids revealed differential effects of ATOH7 and Neurog2 expression on neuronal fate specification. ATOH7 elevation resulted in clear enhancement of the RGC cluster 8, without affecting the photoreceptor cluster 10 ( Figure 8B, also see Figure 6C). In contrast, Neurog2 overexpression not only significantly enhanced the production of photoreceptor cluster 10 and RGC cluster 8, but also generated a distinct RGC cluster 9, which was largely absent in LV-GFP and LV-AEP infected retinal organoids (Figure 8B, also see Figure 6C). Comparison analyses revealed that the two RGC clusters consisted of cells with differential gene expression levels. For example, GAP43 and NSG1 were expressed by cells in both clusters, however, they were frequently detected at higher levels in cluster 8 (Figures 8D,E). Although both clusters expressed ISL1 and POU4F2 (Supplementary Figure 6), we often observed lower levels of POU4F2 and SNCG in many cells of cluster 9 ( Figure 8D). These results suggest that ATOH7 and Neurog2 might differentially influence neuroblast fate specification and/or neuronal differentiation.

Neurogenic Factors Modulate Early Retinal Gene Network
Next, we examined the effects of viral mediated ATOH7f and Neurog2 elevation on endogenous bHLH gene expression. Since the lentiviral vectors did not contain poly-A sequences associated with ATOH7 and Neurog2 cDNAs, and the single cell cDNA libraries were constructed using oligo dT priming, we were able to analyze expression of the endogenous genes using the sc RNA-seq datasets, while excluding transgenes encoded by the viruses. We first examined expression of ASCL1, an early onset bHLH neurogenic factor. Consistent with previous studies (Lu et al., 2020;Sridhar et al., 2020), feature plots showed that ASCL1 was predominantly expressed in pre-neurogenic and neurogenic progenitors (Figure 9A). Viral mediated ATOH7 and Neurog2 expression resulted in significant reduction of ASCL1 + cell from 29.7 to 17.1% and 13.2%, respectively ( Figure 9C); but neither affected the median levels of ASCL1 expression ( Figure 9B). In retinal organoids, ATOH7f and NEUROG2 transcripts were detected in exiting progenitors, neuroblasts, and some postmitotic neurons ( Figure 9A). Both virally expressed ATOH7f and Neurog2 increased the number of endogenous ATOH7-expressing cells from 9.1 to 22.0% among total cells (Figure 9C), and LV-NEP infection also increased the median ATOH7 expression level ( Figure 9B). Interestingly, LV-NEP but not LV-AEP infection caused a twofold increase of endogenous NEUROG2-expressing cells as well as elevated the median expression level (Figures 9B,C), suggesting that Neurog2 positively regulates endogenous NEUROG2 expression. In addition, elevated levels of ATOH7f and Neurog2 both correlated with the increased percentages of cells expressing OLIG2, another bHLH gene in retinal organoids (Hafler et al., 2012; Figure 9C).
In addition to ASCL1, ATOH7, NEUROG2, and OLIG2, sc RNA-seq analysis identified additional Notch signaling pathway genes and bHLH genes affected by ATOH7f and Neurog2 expression (Supplementary Table 1-3 and Supplementary  Figure 7). Violin plots and quantification showed that the Notch ligands DLL1 and DLL4 were upregulated, especially by Neurog2 (Figures 9D,E). In addition, the percentage of cells expressing the Notch signaling effector HES1 was significantly reduced compared to the control LV-GFP infected cells (Figure 9E), while changes in median expression levels were mild ( Figure 9D). In contrast, not only the percentages of HES6 and HEY1 expressing cells were increased (Figure 9E), the median HES6 expression levels were significantly elevated by bot LV-AEP and LV-NEP ( Figure 9D). Since HES1 was predominantly expressed by progenitors, whereas HES6 and HEY1 were upregulated in exiting progenitors and neuroblasts (Supplementary Figure 7), these changes reflected the trend of enhanced neurogenesis. Data from sc RNA-seq also showed that viral expression of ATOH7f or Neurog2 caused substantial increases of several bHLH genes involved in neuronal fate specification or differentiation, including NEUROD1, NEUROD4, NHLH1 and NHLH2 (Figures 9D,E and  Supplementary Figure 7). In all cases, elevated Neurog2 strongly impacted both the number of cells expressing these genes as well as their expression levels.
Together, these results suggested that viral mediated ATOH7f and Neurog2 expression in the developing retinal organoids affected an interactive gene network that plays important roles in cell cycle exit and cell fate specification. We therefore performed STRING network analysis by including bHLH genes and a few selected genes known to be involved in RGC development. The resulting network model included novel regulatory relationships revealed in this study as well as previously reported molecular interactions ( Figure 9F). We also summarize the observed temporal progression of gene expression as retinal organoid cells advanced through the developmental states as defined by our sc RNA-seq analysis ( Figure 9G).

DISCUSSION
In this study we have used ESC-derived organoid cultures as a model system to investigate human embryonic retinal development. The 3D retinal organoids retain the unique molecular signature and correct tissue polarity of the retinal epithelium in vivo, and generate the expected early retinal cell lineages and functional neurons. By transducing human retinal organoids with inducible viral vectors, we were able to determine whether elevating ATOH7 or Neurog2 affected retinal progenitors and neuronal production, and examine the biological processes and the gene network involving these neurogenic factors.
We have identified three classes of retinal progenitors by performing single cell transcriptome analysis. The preneurogenic progenitors express PAX6, VSX2, LHX2, and SOX2, and thus clearly possess the identity of the retinal primordium. Their main GO pathways show cellular metabolic characteristics common to stem cells (Gu et al., 2016), with prominently featured DEGs of the glycolytic pathway, including SLC2A1, GPI, PGK1, and PLOD2, and low mitochondrial respiratory chain genes. This state likely represents a naïve pre-neurogenic state found in the retinal primordium before the onset of neurogenesis and later at the ciliary margin of the mature retina. These pre-neurogenic progenitors express cyclin D1, but have relatively low levels of PCNA and TOP2A, suggesting that they might be slow cycling cells. The pre-neurogenic progenitor category was significantly reduced in LV-AEP and LV-NEP transduced retinal organoids, indicating that expression of these neurogenic factors promotes the transition from the pre-neurogenic to the neurogenic state. This finding is consistent with our previous observation in the embryonic chicken retina, where ectopic human ATOH7 expression in the pre-neurogenic peripheral retina induces precocious neurogenesis ahead of the neurogenic wave front (Zhang et al., 2018). These results indicate that the onset of neurogenic factor expression among naïve retinal progenitors can trigger or accelerate the transition into a neurogenic state in the context of the retinal primordium. Noticeably, the transition from the pre-neurogenic state to the neurogenic state is also accompanied by upregulation of the chromosomal high mobility group genes, such as HMGB2 and HMGN2, reflecting the underlying epigenetic changes accompanying this transition.
Our sc RNA-seq analysis also identified a novel and distinct progenitor state among the proliferating neurogenic progenitors, the cell cycle exiting progenitors. This group of cells still express the characteristic M phase genes such as CDC20 that interacts with the mitosis complex, and MZT1 and PLK1 that organize the mitotic spindles. However, the exiting progenitors show significant upregulation of the Notch ligands DLL1 and DLL4, as well as the bHLH genes ATOH7 and NEUROD1. Both viral mediated ATOH7 and Neurog2 elevation caused expansion of the exiting progenitor state. The developmental trajectory analysis points to a progression from the exiting progenitors to the postmitotic neuroblasts, which as a group exhibits transcriptome profile partially overlapping with all three early neuronal lineages, including VSX1, NHLH1, ATOH7, ONECUT2, NEUROD1, and the growth arrest gene GADD45A. However, transcription profiles of individual cells among neuroblasts are heterozygous, likely reflecting the dynamic gene expression of neuroblasts that serve as an important transitional cell pool poised for terminal cell fate choices and differentiation. Interestingly, we have observed that relative to ATOH7, elevating Neurog2 asserted more potent effects in promoting cell cycle withdrawal, neuroblast expansion, and neuronal differentiation in retinal organoids.
Among the early bHLH factors detected by sc RNA-seq in retinal organoids, a higher percentage of cells expressed ATOH7 (9%) than NEUROG2 (3.2%) endogenously. Although elevating either ATOH7 or Neurog2 caused enhanced neurogenesis, we also detected differential effects of these two factors. For example, at day 39, ATOH7 promoted POU4F + cells, whereas Neurog2 increased ISL1 + cells; but both factors enhanced ISL1 + cells by day 47. Furthermore, the RGC cluster enhanced by ATOH7 elevation shared transcriptomic signatures with RGCs from the control organoids between days 45 and 48. In contrast, Neurog2 induction promoted two groups of RGCs that both expressed ISL1 and POU4F2, but exhibited different levels of known RGC markers including SNCG, POU4F2, and GAP43. Recent transcriptome analysis of mature human RGCs has indeed detected differential expression levels of RGC genes such as POU4F2 among RGC subtypes (Yan et al., 2020). Therefore, the promotion of the two RGC subtypes may reflect the more potent neurogenic effect of Neurog2, which might have accelerated the RGC differentiation program to acquire RGC subtype features. In addition, the effect of elevated Neurog2 to induce ISL1 + neurons preferentially earlier on may also influence outcomes of neuronal differentiation. However, we cannot rule out that virally expressed exogenous Neurog2 has promoted a hybrid neuronal cell type, which does not naturally occur during human retinogenesis. Future research using long term cultures may provide data to address these possibilities. In addition to the RGC phenotypes, ATOH7 increased the HA/AC lineage, consistent with the developmental trajectory. The lack of ATOH7 enhancement on cone photoreceptor production is somewhat unexpected, as ATOH7 is associated with cone cell lineage in the mouse retina (Brzezinski et al., 2012), and we have previously detected a mild enhancement of cone cell genesis by ATOH7 on in the chicken retina (Zhang et al., 2018). In contrast, Neurog2 significantly enhanced photoreceptor production without affecting the HC/AC lineage, which likely reflects the effect of Neurog2-induced NEUROD1 upregulation.
Expression of cell intrinsic factors can be profoundly influenced by extrinsic signaling events during neural development (Gouti et al., 2015). Notch mediated cell-tocell signaling plays important roles during retinogenesis, and disruption of Notch signaling can influence RGC and cone photoreceptor development (Akagi et al., 2004;Jadhav et al., 2006a,b;Yaron et al., 2006;Riesenberg et al., 2009b). In addition, differentiated RGCs have been shown to produce secreted signals that modulate progenitor behaviors. Among the signaling molecules released by differentiated RGCs are Shh, GDF11, and VEGF, all of which promote retinal progenitor proliferation while simultaneously suppressing RGC production (Zhang and Yang, 2001;Kim et al., 2005;Wang et al., 2005;Sakagami et al., 2009). Accumulating evidence indicate that the Notch signaling downstream effector Hes1 may serve as an important integration node for distinct signaling pathways that converge upon retinal progenitors to influence cell proliferation and control neurogenesis. Hes1 is known to directly control cell proliferation by repressing the CDK inhibitor p27 (Kip1) (Bae et al., 2000). The negative feedback of extrinsic signals on RGC genesis is in part mediated by Hes1 suppression of ATOH7 (Hashimoto et al., 2006;Maurer et al., 2014;Miesfeld et al., 2018b). At the present time, the temporal expression sequence and the complex regulatory relationships among the various bHLH factors in the retina are not fully understood, but their elucidation is likely crucial for a better understanding of cell fate selection and the progressive changes in progenitor competent states.
As revealed by our sc RNA-seq data, HES1 is predominantly expressed by both pre-neurogenic and neurogenic progenitors, and is down regulated among exiting progenitors. Consistent with observations in embryonic human retinas and retinal organoids (Lu et al., 2020;Sridhar et al., 2020), our data indicate that ASCL1 level is low in pre-neurogenic progenitors, but is upregulated and expressed by neurogenic progenitors, and then down-regulated in exiting progenitors. This observation is consistent with other studies using sc RNA-seq to analyze embryonic human retina and early stage human retinal organoids (Sridhar et al., 2020). In the developing mouse retina, Ascl1 can block cell cycle exit but not specify RGC fate (Hufnagel et al., 2013), while in the mature retina ASCL1 promotes cell cycle re-entry and regeneration of new neurons from Muller glia (Pollak et al., 2013;Ramachandran et al., 2015;Jorstad et al., 2017). Furthermore, forced expression of ASCL1 in pluripotent stem cells can promote neurogenesis . These data together support that ASCL1 plays a role in establishing competence for the neurogenic state. Intriguingly, our results show that viral mediated ATOH7 and Neurog2 expression decrease both HES1 and ASCL1 expression. Concomitantly, ATOH7 and Neurog2 significantly upregulate the expression of the HES1 inhibitor HES6 (Bae et al., 2000) among exiting progenitors. These data suggest that ATOH7 and NEUROG2 comprise a selective group of neurogenic factors that can dampen the effect of HES1 and ASCL1 on maintaining cell proliferation and relieve the inhibitory effect of HES1 on neurogenesis. Moreover, viral mediated ATOH7 and Neurog2 expression significantly increase transcripts of other bHLH factors including NEUROD1, which initiates its expression among exiting progenitors, as well as NEUROD4, NHLH1, and NHLH4, which are expressed by postmitotic neurons. Future investigations are necessary to determine whether the regulatory effects asserted by ATOH7 and Neurog2 are direct or indirect.
Results from this study support the hypothesis that high levels of AOTH7 or NEUROG2 trigger a withdrawal from the cell cycle which leads to the birth of a neuroblast. Critical questions remain regarding how the level of ATOH7 is regulated in neurogenic progenitors to result in higher levels among cell cycle exiting progenitors. It is known that the ATOH7 gene contains enhancers that mediate direct positive regulation by Pax6 (Riesenberg et al., 2009a). In Hes1 mutants, Atoh7 is precociously expressed along with the formation of RGC and HC/AC (Lee et al., 2005), indicating that Hes1 negatively regulates Atoh7 expression. Our sc RNA-seq analysis show that viral mediated ATOH7 expression elevates endogenous ATOH7 without affecting endogenous NEUROG2 expression, whereas Neurog2 elevation leads to increases in both endogenous ATOH7 and NEUROG2 expression. These results reveal the novel finding that ATOH7 and NEUROG2 are both under positive autoregulation, as well as confirm cross-regulation of ATOH7 by NEUROG2 (Matter-Sadzinski et al., 2001;Hufnagel et al., 2010). Since only a fraction of ATOH7 protein expressing cells appear to co-express Neurog2 protein during a given window of time (Miesfeld et al., 2018a), progenitors that have coincidental expression of both factors are more likely to reach higher level of Atoh7. It is known that ATOH7 dosage and expression levels can affect RGC production Chiodini et al., 2013;Zhang et al., 2018). We therefore propose that integrated cell extrinsic signals and interacting cellintrinsic factors could converge, resulting in stochastic expression of ATOH7, and thus enabling a limited subset of progenitors that have reached a threshold level of ATOH7 to exit the cell cycle and initiate downstream RGC and/or cone photoreceptor differentiation programs. Further investigations will be necessary to determine the ATOH7 protein threshold, the time, and cellular events involved. The sc RNA-seq analysis enables us to survey multiple genes and to construct a gene regulatory network model (Figure 9F) that integrates our new findings and previously known regulations. This model focuses on bHLH factors expressed during RGC development, but does not exclude other genes involved in the developmental process. In fact, RGC and other neuronal fate determinations are known to be regulated by multiple transcription factors (Rehemtulla et al., 1996;Koike et al., 2007;Oh et al., 2007;Jiang et al., 2013;Wu et al., 2013;Chang et al., 2017;Yamamoto et al., 2020). Interestingly, a recent report has shown that in the absence of Atoh7, mouse retinal progenitors can exit the cell cycle toward an RGC fate, but cannot fully differentiate into mature RGCs (Wu et al., 2021).
Recent single-cell analyses have shown evolutionarily conserved gene expression patterns during retinogenesis, but also revealed species-specific patterns between human and mouse retinas (Liang et al., 2019;Peng et al., 2019;Lu et al., 2020), including a role for ATOH7 in late stage photoreceptor specification. Our study shows that elevating ATOH7 and Neurog2 expression in human retinal organoids significantly enhances early retinogenesis, which can serve as a useful approach to produce authentic human RGCs for studying development and degenerative diseases. Furthermore, our single cell transcriptome analysis provides novel insights into the interactive network of bHLH factors and their functions during human retinogenesis.

Lentiviral Construction and Production
The inducible lentiviral vector plasmids LV-GFP and LV-NEP (YS-TetO-FUW-Ng2-P2A-EGFP-T2A-Puro) were generous gifts from Dr. Thomas Sudhof (Zhang et al., 2013). The human ATOH7 cDNA was obtained as described previously (Zhang et al., 2018). The LV-ATOH7f vector was constructed by replacing the EGFP gene in the LV-GFP vector with ATOH7 cDNA fused to the Flag epitope tag at the c-terminus. The LV-AEP vector was constructed by custom synthesizing the continuous open reading frame of ATOH7f-P2A-EGFP-T2A-Puro by GenScript, and replacing the EGFP gene in the LV-GFP vector. All lentiviral vector plasmids were verified by DNA sequencing. The LV-rtTA lentiviral vector (FUW-M2rtTA) with the Ubi promoter was obtained from AddGene (Plasmid #20342). Lentiviral stocks were produced by co-transfection of HEK 293T cells with a given viral vector DNA and the third-generation lentiviral helper plasmids with VSVG pseudotyping as described Zufferey et al., 1998;Hashimoto et al., 2007). The 293T cell medium DMEM containing 10% fetal bovine serum (Sigma-Aldrich, 12103C) was changed to serum free CD293 (Thermo Fisher, 11913) 1 day post transfection, and viral supernatants were harvested every 24 h. Combined viral stocks were concentrated by ultracentrifugation as previously described (Hashimoto et al., 2007).

Lentiviral Infection and Transgene Induction
Retinal organoids were infected by different lentiviruses in conjunction with LV-rtTA three times between days 23 and 40. TetO promoter induction was carried out by adding doxycycline (Dox) to a final concentration of 2 µg/ml (Sigma-Aldrich, D3072) according to experimental designs as indicated in the results.

Immunohistochemistry and Imaging
For live imaging of whole mounts retinal organoids, EGFP signals were first captured using a Leica MZ10F fluorescent dissecting microscope, followed by image acquisition using an Olympus Flowview FV1000-IX81 (inverted) scanning laser confocal microscope. For whole mount immunolabeling, retinal organoids were fixed with 4% paraformaldehyde (PFA) in PBS for 30 min, followed by primary and secondary antibody incubations overnight at 4 • C, with extensive washes in between (Zhang et al., 2018). Cryosections (14 µm) or attached cultures were processed for immunolabeling as previously described (Zhang et al., 2018).

Cell Marker Quantification and Statistics
Lentivirus infected retinal organoids at day 39 were pooled (10-15 organoids in each sample n), dissociated into single cell suspensions using trypsin (Sigma-Aldrich, T-9935), and immunolabeled as previously described (Sakagami et al., 2009). Flow cytometry was performed using LSRII Analytic Flow Cytometer for cell marker analyses. Quantification of FACS data was performed using FlowJo software (Tree Star, Inc.). In addition, pooled retinal organoids at day 47 were dissociated with trypsin and plated as a monolayer for 3 h followed by immunolabeling for the various cell markers listed above. Monolayer cell quantification was performed by counting marker-positive cells in multiple fields of independent samples (n = 3) using captured confocal images. Bar graphs were constructed using Prism (Graphic Pad). Ordinary one-way ANOVA with Tukey's multiple comparison test was used for statistical analysis of cell marker quantifications, with p < 0.05 considered significant.

Single Cell cDNA Library Preparation and Sequencing
Distinct pools of H9 ES cell-derived retinal organoids (12-20 retinal organoids/pool) co-infected by LV-rtTA and LV-EGFP, LV-AEP, or LV-NEP were induced by Dox and dissociated between days 45 and 48 using trypsin and manual trituration. Dissociated cell suspensions were subjected to fluorescence activated cell sorting using FACSAriaII (BD Biosciences). Noninfected retinal organoid cells were used to set thresholds for selecting EGFP-positive cells. Sorted EGFP-positive cells were collected in HBSS without Ca 2+ and Mg 2+ (Thermo Fisher, 14170-112) containing 1% FBS and 0.4% BSA. The cells were washed with PBS containing 0.04% BSA, then counted with Countess II Cell Counter (Thermo Fisher).
Automated single-cell capture, barcoding, and cDNA library preparation were carried out using 10X Genomics Chromium Controller with Chromium Single Cell 3' Library & Gel Bead Kit v2 reagents, with 12 cycles of cDNA amplification and 12 cycles of library amplification, following the manufacturer's instructions. Qubit dsDNA Assay kit (Life Technologies) and TapeStation 4200 (Agilent) were used to assess the quality and concentration of the libraries. Illumina NovaSeq6000 S2 paired-end 2 × 50 bp mode was used to sequence the libraries.
Single Cell RNA-Sequencing Data Processing and Quality Control 10X Genomics Cell Ranger version 2.1.1 was used to demultiplex the raw base calls into FASTQ files (cellranger mkfastq). Spliced Transcripts Alignment to a Reference (STAR) version 2.5.1b (cellranger count) was used to perform sequence alignments to the reference human genome (GRCh38), barcode counts, and UMI counts to yield summary reports and t-Stochastic Neighboring Embedding (t-SNE) dimensionality reduction. For downstream analyses, cells with a number of unique molecular identifiers (UMI) > 2,500 per cell and < 0.1% mitochondrial gene expression were used. For LV-GFP, LV-AEP, LV-NEP samples, the mean reads per cell ranged from 139,000 to 195,000, with mean gene per cell ranging from 2935 to 3079. The resulting total single cell counts used for analysis were 3004 for LV-EGFP, 2063 for LV-AEP, and 3909 for LV-NEP infected samples.

Single Cell RNA-Sequencing Data Analysis and Visualization
The analysis of sc RNA-seq data was performed using Seurat R package 1 (Butler et al., 2018;Stuart et al., 2019). Clustering of cells was performed by using Seurat FindCluster function (top 20 principal components, resolution 0.8) that implements the shared nearest neighbor modularity optimization algorithm. Non-linear dimensionality reduction using UMAP (Uniform Manifold Approximation and Projection) was applied for the visualization of cells in two-dimensional space. Feature plots of known genes were used to designate clusters observed in the UMAP space into six major cell categories/states. Cell counts of each category were obtained using custom R code.
Pseudotime developmental progression of cell states was obtained by using the Monocle R package (version 2) to process the datasets with cell labels corresponding to the six cell categories and visualized as UMAPs. Pseudotime cell cycle progression and cell fate adoption analysis was performed using Slingshot R (version 1.6.1) by combining the LV-GFP, LV-AEP, and LV-NEP sc RNA-seq datasets and assigning the start point as neural stem cells and end points as differentiated neuronal cell types.
Differentially expressed genes (DEGs) were identified using edgeR with significantly enriched genes in each cell category defined as those with adjusted p < 0.05 for LV-GFP, LV-AEP, and LV-NEP data sets (Supplementary Tables 1-3). The top 10 enriched DEGs in each cell category were defined as those with adjusted p< 0.05, and log fold change > 1.5 ( Supplementary  Tables 4-6). Heatmaps of the top 10 DEGs were generated using the Seurat package. Volcano plots were generated using EnhancedVolcano R package to show p-values and fold changes of DEGs between two datasets. Gene Ontology (GO) Enrichment analysis was performed using ShinyGO v0.61 2 (Ge et al., 2020) and the Homo Sapiens background using p-value (FDR) cutoff at 0.05. The top 25 DEGs from each cell category within the LV-GFP dataset were used as inputs, and the redundancy of the output biological processes was manually reduced to the most predominant GO terms.
Feature plots of individual gene expression patterns in different cell clusters were presented as UMAPs. Violin plots for individual genes in all cell clusters were constructed to show expression levels and cell distributions. Kruskal-Wallis one-way ANOVA rank sum test and Tukey-Kramer-Nemenyi all-pairs test were used for statistical analysis, taken into consideration of both gene expression levels and cell numbers between different samples, with p < 0.05 considered significant. The statistical tests were performed on R Studio using "PMCMRPlus" (Kessler et al., 2020) and "FSA" (Derek H. Ogle et al., 2021) packages. Data were plotted using R Studio "ggplot2" (Hadley, 2016) and "ggsignif " (Ahlmann-Eltze, 2019) packages. Cells with gene expression level < 0.2 were exclude from the violin plots and statistical analyses.
STRING analysis exploring protein-protein association network was performed using human protein database 3 (version 11.0) by inputting relevant genes involved in retinal development. The schematic network model shows known molecular interactions reported previously and new regulatory relationships described in this study.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The authors have deposited the sc RNA-seq datasets into the public archive Dryad (doi: 10.5068/D1ZD5S).

AUTHOR CONTRIBUTIONS
X-JY conceived the project, designed lentiviral vectors, and wrote the manuscript. XZ carried out human retinal organoid and neuronal cultures, performed immunohistochemistry and imaging, and FACS analysis. XZ and KN prepared viral stocks and sc RNA-seq libraries. TN, IM, and MP performed bioinformatic analysis. TN performed statistical analysis. JG and SB performed electrophysiological recordings. XZ, TN, and X-JY prepared the figures and analyzed the data. All authors reviewed and commented on the manuscript. After 24-hour of Dox induction, both LV-AEP and LV-NEP infected cells showed a tendency toward localizing to the inner layer, compared to the control LV-GFP virus infected retinal organoids. This trend became more pronounced after 48-hour Dox induction. By 72-hour after the onset of Dox treatments, the majority of GFP + cells were concentrated in the inner layer of the retinal organoids. In contrast, most LV-GFP infected cells remained a ventricular zone distribution pattern after 48-and 72-hour induction. Scale bars, 200 µm for the lower magnification; 100 µm for the insets.