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

Zhang, Xiangmei; Mandric, Igor; Nguyen, Kevin H.; Nguyen, Thao T. T.; Pellegrini, Matteo; Grove, James C. R.; Barnes, Steven; Yang, Xian-Jie

doi:10.3389/fcell.2021.653305

ORIGINAL RESEARCH article

Front. Cell Dev. Biol., 13 May 2021

Sec. Stem Cell Research

Volume 9 - 2021 | https://doi.org/10.3389/fcell.2021.653305

This article is part of the Research Topic Mechanisms of Cellular Differentiation, Organ Development, and Novel Model Systems View all 25 articles

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

$\r\nXiangmei Zhang$ Xiangmei Zhang¹

Igor Mandric²

Kevin H. Nguyen¹

Thao T. T. Nguyen¹

Matteo Pellegrini²

James C. R. Grove¹

Steven Barnes^1,3

Xian-Jie Yang^1,4*

¹Department of Ophthalmology, Stein Eye Institute, University of California, Los Angeles, Los Angeles, CA, United States
²Department of Molecular, Cell, and Developmental Biology, University of California, Los Angeles, Los Angeles, CA, United States
³Doheny Eye Institute, University of California, Los Angeles, Los Angeles, CA, United States
⁴Molecular Biology Institute, University of California, Los Angeles, Los Angeles, CA, United States

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 co-expression 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 (Prasov and Glaser, 2012). 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 (Prasov et al., 2012).

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.

Results

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.

FIGURE 1

Figure 1. Inducible lentivirus-mediated gene expression in human ES cell-derived retinal organoids. (A) Schematics of lentiviral vectors used in the study. A vector expressing rtTA under the control of the constitutive Ubiq promoter was used to co-infect human ES cells with lentiviral vectors encoding the inducible TetO promoter. Puro, puromycin-resistant gene; P2A and T2A, intervening 2T sequences in frame with the coding sequences of the other genes in LV-NEP and LV-AEP. The Flag epitope tag was fused to the c-terminus of Neurog2 and ATOH7 coding regions. (B) Expression of Neurog2 and ATOH7 proteins as detected by immunofluorescent microscopy in H9 ES cells co-infected with LV-rtTA and either LV-NEP, LV-AEP, or LV-ATOH7f after 36 h of Dox induction. Scale Bar for all, 20 μm. (C) Experimental time course of human retinal organoid derivation and lentiviral infection. (D,E) Confocal images (D) and merged confocal and bright field images (E) show distribution of viral infected GFP⁺ cells in live retinal organoids at day 36 after 6 days of Dox induction. Names of lentiviruses are abbreviated as: GFP for LV-GFP, NEP for LV-NEP, AEP for LV-AEP. Scale Bars, 100 μm.

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 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.

FIGURE 2

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⁺ PCNA^– postmitotic 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.

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²⁺ channels as well as spontaneous and provoked electrophysiological excitability in vitro during whole cell patch clamp recording, characteristic of native RGCs (Supplementary Figure 3). TTX-sensitive Na⁺ channels, TEA-sensitive K⁺ channels, and Cd²⁺-sensitive Ca²⁺ 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).

FIGURE 3

Figure 3. Effects of viral mediated neurogenic factor expression on retinal ganglion cell genesis. Confocal images of retinal organoid sections immunolabeled with retinal ganglion cell (RGC) markers. (A) Co-labeling of viral markers with a pan-POU4F antibody in LV-NEP and LV-ATOH7f infected retinal organoids at day 43. (B) Merged images of viral marker GFP co-labeling with NF145 and NeuN at day 36, and DCX at day 46. Scale Bars, 20 μm.

FIGURE 4

Figure 4. Impacts of elevating neurogenic factor expression on retinal ganglion cell genesis. (A) Immunofluorescent images of attached retinal organoid cultures infected with LV-GFP or LV-ATOH7f after 4-day Dox induction at day 39. Co-labeling for viral vector markers with POU4F and NF145 at day 38, or POU4F1 and RBPMS at day 41 shows high correspondence of LV-ATOH7f infection and RGC marker expression. Scale bars, 20 μm. (B) Representative flow cytometry profiles of dissociated cells at day 39 from LV-GFP, LV-ATOH7f, and LV-NEP infected retinal organoids after 4-day Dox induction. Cells were co-labeled for RGC marker POU4F and viral marker Flag for LV-ATOH7f, or GFP for LV-GFP and LV-NEP, respectively. (C,D) Bar graphs show flow cytometry quantification of RGC markers POU4F and ISL1 among total cells (C) and viral infected cells (D) at day 39. (E) Quantification of cultured monolayer cells from LV-GFP, LV-AEP, and LV-NEP infected retinal organoids at day 47 after 7-day Dox induction. Percentages of RGC marker ISL1, DCX, and NeuN positive cells among viral infected GFP⁺ cells are shown. 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.

We also performed flow cytometry analyses of dissociated retinal organoid cells to quantify the effects of viral mediated ATOH7 or Neurog2 expression on RGC genesis (Figures 4B–D). After 4 days of Dox induction at day 39, ATOH7 expression led to a twofold increase of POU4F⁺ cells, from 6.6 ± 0.68% of total cells in LV-GFP infected organoids to 13.7 ± 0.59% of total cells in LV-ATOH7f infected organoids (Figure 4C). Among LV-ATOH7f infected cells, 57.1 ± 1.9% were POU4F⁺ compared to 6.1 ± 1.1% among LV-GFP infected cells (Figure 4D). However, LV-ATOH7f did not increase expression of the RGC marker ISLET1 at this stage (Figure 4C). In the parallel analysis, LV-NEP infection did not significantly promote POU4F⁺ cells, but instead increased ISLET⁺ cells from 5.9 ± 0.95% to 11.1 ± 0.89% of total cells (Figure 4C), and from 4.4 ± 1.1% to 18.4 ± 1.6% among viral infected cells (Figure 4D). Similar analyses at day 47 following a 7-day Dox induction showed significantly increased co-labeling with the RGC markers ISLET1 (>10-fold of 6.3%), DCX (>4-fold of 21%), and NeuN (>20-fold of 4%) of LV-AEP and LV-NEP infected cells compared with the control LV-GFP infected cells (Figure 4E). These results demonstrate that elevated ATOH7 or Neurog2 expression promoted human RGC production in retinal organoids.

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.

FIGURE 5

Figure 5. Characterization of lentivirus-infected human retinal organoids by single cell RNA-seq. (A) UMAP visualization of cell clustering based on sequencing of single cell cDNA libraries from lentivirus infected retinal organoids between culture days 45–48 after 8-day Dox induction. The numbers of single cell that passed quality control and used for downstream analysis are 3004 for LV-GFP, 2063 for LV-AEP, and 3909 for LV-NEP samples. (B) Dot plots show assignments of UMAP clusters into seven major cell categories based on expression of known genes for different cell states or types during early retinal development. PNP, pre-neurogenic progenitor; NP, neurogenic progenitor; EP, cell cycle-exiting progenitor; NB, neuroblast; RGC, retinal ganglion cell; PR, photoreceptor; HC/AC, horizontal cell and amacrine cell. Minor cell clusters in LV-GFP and LV-NEP infected retinal organoids that do not show characteristic retinal gene expression are designated as “unknown.” The average expression levels are represented by the intensity bar, and the percentages of cells expressing a given gene in a cluster are indicated as the dot size.

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 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.

FIGURE 6

Figure 6. Distinct cell categories based on single cell RNA-seq profiles and their developmental trajectory. (A) UMAP presentations based on cell category assignments for single cell clusters from LV-GFP, LV-AEP, and LV-NEP infected retinal organoids. Each cell category is represented by a single color. The side legend shows the color codes used to represent distinct cell categories. (B) Pseudotime trajectory analysis using Monocle of different cell categories from LV-GFP, LV-AEP, and LV-NEP infected retinal organoid. (C) Bar graph shows the percentages of different cell categories among total cells in LV-GFP, LV-AEP, and LV-NEP infected retinal organoid cells. PNP, pre-neurogenic progenitor; NP, neurogenic progenitor; EP, cell cycle-exiting progenitor; NB, neuroblast; RGC, retinal ganglion cell; PR, photoreceptor; HC/AC, horizontal cell and amacrine cell.

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 pre-neurogenic 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.

FIGURE 7

Figure 7. Effects of elevating neurogenic factor expression on developmental transitions. (A) Heatmaps of the top 10 DEGs in each cell type category from cells of LV-GFP, LV-AEP, and LV-NEP infected retinal organoids. The side legend indicates the color codes for cell category identities atop the heatmaps. (B) Volcano plots show significant DEGs in LV-GFP infected retinal organoids during the three transitions: from the neural stem cell to the neurogenic progenitor state, from the cycling progenitor to the exiting progenitor state, and from the exiting progenitor to the neuroblast states. Genes with p < 0.01 and fold change > Log2 0.75 are shown as red dots. (C) Bar graphs show the predominant biological processes using GO pathway analysis for each cell category of LV-GFP infected retinal organoids. The false discovery rates (FDR) and the gene numbers among the top 25 DEGs associated with a given GO term are shown. PNP, pre-neurogenic progenitor; NP, neurogenic progenitor; EP, cell cycle-exiting progenitor; NB, neuroblast; RGC, retinal ganglion cell; PR, photoreceptor; HC/AC, horizontal cell and amacrine cell.

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 cytokinesis (Lee et al., 2008; Supplementary Figure 6). Slingshot analysis performed for the combined dataset illustrated the developmental trajectory of cell clusters from the pre-neurogenic state through the different phases of the cell cycle toward cell cycle exit (Figure 8A). Although cell clusters 5 and 6 both expressed M phase genes, they clearly had distinct transcript 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).

FIGURE 8

Figure 8. Influence of neurogenic factors on cell cycle progression and retinal ganglion cell differentiation. (A) UMAP of the combined sc RNA-seq dataset from LV-GFP, LV-AEP, and LV-NEP infected retinal organoids. Each cell cluster is assigned to an identity as shown on the side color legend based of feature plots of known genes. Line edges connecting individual clusters are generated by SLINGSHOT analysis to show the predicted developmental trajectory of the entire cell cohort. PNP, pre-neurogenic progenitor; NP, neurogenic progenitor; EP, cell cycle-exiting progenitor; NB, neuroblast; RGC, retinal ganglion cell; PR, photoreceptor; HC/AC, horizontal cell and amacrine cell. (B) UMAPs show cell clusters for individual samples from LV-GFP, LV-AEP, and LV-NEP infected retinal organoids with the same cell identity assignments as in the combined UMAP. (C) Distribution of cells in different phases of the cell cycle as percentages of total proliferative progenitor cells. The M phase cells include the exiting progenitors (cluster #6). (D) Heatmap shows differential expression levels of genes in the two RGC clusters 8 and 9. (E) Volcano plot shows DEGs between the two RGC clusters 8 and 9. Genes with p< 0.01 and fold change > Log2 0.5 are shown as red dots.

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).

FIGURE 9

Figure 9. Regulatory relationship among genes playing roles in early human retinogenesis. (A) Feature plots show distribution of ASCL1, ATOH7, and NEUROG2 expressing cells among cell clusters in LV-GFP, LV-AEP, and LV-NEP infected retinal organoids. (B,D) Violin plots show comparisons of LV-GFP (green) versus LV-AEP (red) or LV-NEP (blue) induction on expression levels of endogenous genes in the retinal organoids. The box within the violin represents the middle 50% of the data, and the yellow line within the box indicates the median expression level. Statistical analysis generated P-values were based on both the cell counts and expression levels. ***p < 0.001, **p < 0.01, *p < 0.05. (C,E) Bar graphs show percentages of cells expressing endogenous genes among total viral infected cells in LV-GFP, LV-AEP, and LV-NEP samples. (F) Schematic model based on STRING analysis summarizes a gene network involved in cell cycle exit and early retinogenesis. Gray edges represent protein-protein associations. Positive and negative regulatory relationships are indicated as arrows and short bars, respectively. Previously known molecular interactions from curated databases are shown as blue arrows and short bars. Effects of lentiviral induced ATOH7f and Neurog2 on endogenous human gene expression in retinal organoids are indicated as green arrows and red bars. Looping arrows indicate positive auto-regulation. (G) Summary of sequential expression of bHLH genes (red) and other key genes (blue) during early human retinal organoid development. PNP, pre-neurogenic progenitor; NP, neurogenic progenitor; EP, cell cycle-exiting progenitor; NB, neuroblast; RGC, retinal ganglion cell; PR, photoreceptor; HC/AC, horizontal cell and amacrine cell.

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 pre-neurogenic 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-to-cell 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 (Wang et al., 2020). 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 (Prasov et al., 2012; Chiodini et al., 2013; Zhang et al., 2018). We therefore propose that integrated cell extrinsic signals and interacting cell-intrinsic 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.

Materials and Methods

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 (Dull et al., 1998; 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).

Human Retinal Organoid Derivation

Human H9 ES cells were cultured and passaged on Matrigel (Corning, 356231) coated dishes in mTeSR1 medium (Stemcell Technologies, 05850). Retinal organoids were generated based on a previously described protocol (Ohlemacher et al., 2016) with modifications. At the start of the culture (day 0), H9 ES cells (at 80–90% confluency) were enzymatic detached using dispase (1 mg/ml, Stemcell Technologies, 07923). Detached cells were transferred into medium at 3:1 ratio of mTeSR1 to Neural Induction Medium (NIM) that consists of DMEM/F12 with 1x N2 supplement (Thermo Fisher, 17502048), 1x non-essential amino acid (NEAA; Thermo Fisher, 11140050), 2 μg/ml Heparin (Thermo Fisher, H7482), and 1x Antibiotic-Antimycotic (Anti-Anti; Thermo Fisher, 15240112) in low-attachment plates (Corning, 3471) to allow the formation of embryonic bodies (EBs). During the next 3 days, the medium was replaced daily with the ratio of NIM to mTeSR1 increased to 100%. At day 6, human BMP4 (R&D Systems, 314-BP-010) was added to the NIM medium to a final concentration of 55 ng/ml (Kuwahara et al., 2015). At day 7, EBs were collected and seeded in NIM as adherent cultures in 6-well dishes (Corning, 3516) till day 16. At day 16, the visible neural rosettes formed from attached EBs were manually lifted, collected, and further cultured as suspensions in Retinal Differentiation Medium (RDM) consisting of DMEM to F12 at 3:1, 1x B27 supplement (Thermo Fisher,1754044), 1x NEAA, and 1x Anti-Anti. During the 7 days after lifting neural rosettes, 5 μM SU-5402 (Sigma-Millipore, SML0443) and 3 μM GSK inhibitor CHIR99021 (Stemgent, 04-0004) were added to RDM. At day 20, the translucent optic vesicle-like structures were manually separated from the rest of the suspension culture, collected, and cultured as retinal organoids in RDM. From day 24, 10% FBS (Sigma-Aldrich, 12103C), 100 μM Taurine (Sigma-Aldrich, 0625), and 500 μM retinoic acid (Sigma-Aldrich, R2625) were added, and the medium was changed twice a week.

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.

Attached and Dissociated Retinal Organoid Cultures

Retinal organoids between days 30 and 35 were cut into small pieces (0.1–0.5 mm) and plated on Matrigel coated culture dishes, or dissociated with Trypsin (Sigma-Aldrich, T9935) to single cells and plated on poly-D-lysine and laminin coated glass coverslips (Corning, 354087). After attachment, retinal organoid cells were cultured in RDM or BrainPhys neuronal medium (Stemcell Technologies, 05790) with SM1 (Stemcell Technologies, 05711) and N2 supplements (Thermo Fisher, 17502048), 20 ng/ml BDNF (PeproTech, 450-02), 20 ng/ml GDNF (Stemcell Technologies, 78058), 1 mM dibutyryl cyclic-AMP (Stemcell Technologies, 73882), and 200 nM ascorbic acid (Stemcell Technologies,72132) till desired time, followed by immunofluorescent labeling or electrophysiological recordings.

Electrophysiological Recording

Electrophysiological recordings were performed using dissociated retinal organoid cells cultured as a monolayer on glass coverslips between days 40 and 45 using previously described methodologies (Hirooka et al., 2002; Hartwick et al., 2004; Lalonde et al., 2006; Sargoy et al., 2014). Whole cell patch clamp was performed at room temperature using an Axopatch 200B amplifier controlled by pClamp 11 data acquisition software (Molecular Devices). The pipette solution contained 20 mM KCl, 120 mM K-gluconate 0.1 mM CaCl₂, 1 mM EGTA, 10 mM HEPES, 3 mM Mg-ATP, 0.2 mM Li-GTP, and 8 mM phosphocreatine, at pH7.2. The bathing solution (for recording K⁺ currents, Na⁺ currents and for current clamp) contained 125 mM NaCl, 3 mM KCl, 2 mM CaCl₂, 1.25 mM NaH₂PO₄, 1 mM MgCl₂, 25 mM NaHCO₃ and 10 mM glucose bubbled continuously with 95% O₂—5% CO₂. A high barium external solution for Ca²⁺ channel current recordings containing 110 mM choline chloride, 5 mM KCl, 1 mM MgCl₂ 7 mM BaCl₂, 15 mM TEACl, 0.1 mM 4-aminopyridine, 20 mM glucose, 10 mM HEPES and 1 μM tetrodotoxin, adjusted to pH 7.4, was used with a CsCl intracellular solution containing 140 mM CsCl, 1 mM CaCl₂, 11 mM EGTA, 2 mM MgCl₂, and 10 mM HEPES, at pH7.2.

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).

The following primary antibodies were used: mouse anti-GFP (1:200; Millipore, MAB#3580); goat anti-GFP (1:200; Rockland Inc., 600-101-215); rabbit anti-GFP (1:200; Rockland Inc., 400-401-215); mouse anti-FLAG (Clone M2) (1:500; Sigma-Aldrich, F3165); rabbit anti-ATOH7 (1:100; NovusBio, NBP1-88639); goat anti-Ngn2 (1:200; Santa Cruz Biotechnology, sc-19233); mouse anti-BrdU (1:1; GE Healthcare, RPN202); mouse anti-PCNA (1:500; Sigma-Aldrich, P3825); rabbit anti-phospho-histone 3 (ser10) (1:3,000; Upstate Biotechnology, 06-570); rabbit anti-Pax6 (1:200; Chemicon, ab5409); goat anti-CHX10/VSX2 (N-18) (1:50; Santa Cruz Biotechnology, sc-21690); goat anti-BRN3a (1:100; Santa Cruz Biotechnology, sc-31984); goat anti-pan BRN3 (1:50; Santa Cruz Biotechnology, sc-6026); mouse anti-Iselt1 (1:10; Developmental Study Hybridoma Bank, 39.4D5); mouse anti-doublecortin (E-6) (1:50; Santa Cruz Biotechnology, sc-271390); rabbit anti-NeuN (1:200; Abcam, ab177487); rabbit anti-NF145 (1:750; Millipore, AB1987); rabbit anti-RBPMS (1:200; Ye et al., 2018). Secondary antibodies used were: Alexa Fluor 488 donkey anti-mouse (1:500; Thermo Fisher, A32766); Alexa Fluor 488 donkey anti-rabbit (1:500; Thermo Fisher, A32790); Alexa Fluor 488 donkey anti-goat (1:500; Thermo Fisher, A32814); Alexa Fluor 594 donkey anti-mouse (1:500; Thermo Fisher, A32744); Alexa Fluor 594 donkey anti-rabbit (1:500; Thermo Fisher, A32754); Alexa Fluor 594 donkey anti-goat (1:500; Thermo Fisher, A32758); Alexa Fluor 647 donkey anti-mouse (1:500; Thermo Fisher, A32787); Alexa Fluor 647 donkey anti-rabbit (1:500; Thermo Fisher, A32795); Alexa Fluor 647 donkey anti-goat (1:500; Thermo Fisher, A32849).

Slides were mounted with Fluoro-Gel (Electron Microscopy Sciences, 17985-10) after staining with DAPI (Sigma-Aldrich, D9542). Confocal images were acquired using an Olympus Flowview FV1000-BX61 (upright) scanning laser microscope with Plan-APO objectives. Images were arranged using Adobe Photoshop.

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). Non-infected retinal organoid cells were used to set thresholds for selecting EGFP-positive cells. Sorted EGFP-positive cells were collected in HBSS without Ca²⁺ and Mg²⁺ (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¹ (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 progressions 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² (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³ (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.

Funding

This work was supported by NIH grant R01EY026319 to X-JY, NIH core grant P30EY000331, and an unrestricted grant from the Research to Prevent Blindness to the Department of Ophthalmology at University of California Los Angeles.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fcell.2021.653305/full#supplementary-material

Supplementary Figure 1 | Characterization of human ES cell-derived 3D retinal organoids. (a) Bright field image shows morphology of a group of H9 ES cell-derived 3D retinal organoid at Day 33.Scale bar, 1 mm. (b) RT-PCR assay detects expression of eye field and neural retina genes at Day 16 andDay 29 in H9 ES cell-derived cultures. (c) Whole mount images show a 3D retinal organoid coimmunolabeledfor PAX6 and VSX2 at Day 24. The top panels show low magnification images of theentire retinal organoid (scale bar, 100 μm), and the bottom panels show confocal images with nuclearlabeling of PAX6 and VSX2 in retinal progenitor cells (scale bar, 20 μm).

Supplementary Figure 2 | Effect of Dox induction time course on retinal organoid development. Whole mount images show effects of different Dox induction durations on retinal organoiddevelopment. After 24-hour of Dox induction, both LV-AEP and LV-NEP infected cells showed atendency toward localizing to the inner layer, compared to the control LV-GFP virus infected retinalorganoids. This trend became more pronounced after 48-hour Dox induction. By 72-hour after theonset of Dox treatments, the majority of GFP₊ cells were concentrated in the inner layer of the retinalorganoids. In contrast, most LV-GFP infected cells remained a ventricular zone distribution pattern after48- and 72-hour induction. Scale bars, 200 μm for the lower magnification; 100 μm for the insets.

Supplementary Figure 3 | Electrophysiological properties and functionality of Human H9 ES cell-derivedretinal neurons. (a) Bright field view of the dissociated cell culture derived from retinal organoids at Day 40. Scale bar, 20μm. (b, c) A neuron being patch recorded is shown in bright field (b) and after filling with lucifer yellowdye (c). (d-i) Whole cell patch clamp recording of dissociated neurons (Day 37-40) derived from H9 ESretinal organoids cultured as a monolayer. d Whole cell voltage clamp of a cell with multipolar neuritesstepped from −60 mV to +30 mV in 10 mV steps of 40 ms duration. (e) Example of well-clamped I_Naisolated by digital substraction following block with 100 nM TTX. Steps from −90 to +20 mV are shown. (f)Outward K₊ currents isolated by digital substraction following 10 mM TEA application. Steps from −70 to+30 mV are shown. (g) Train of action potentials elicited with depolarizing current in current clampedcells having large I_Na. (h) Phasic action potential generation in cells having smaller I_Na. (i) Summary ofblock of peak INa and I_K at +40 mV by TTX and TEA, respectively. TTX blocked 92% of the transient inwardcurrent (n=7) and TEA blocked 62% of the sustained outward current at +40 mV (n = 6).

Supplementary Figure 4 | Fluorescent activated cell sorting of lentivirus infected retinal organoids for single cell RNA-sequencing. (a) FACS profiles of dissociated retinal organoids infected with LV-GFP, LV-AEP, and LV-NEP between Day45-48 after 8 day ox treatment. Cells from non-infected retinal organoids were used as negativecontrols to set the thresholds for GFP₊ cells. FACS enriched GFP₊ cells were used for single cell RNA-seqanalyses. (b) Violin plots show the numbers of genes detected in single cell RNA-seq analyses using 10XGenomics Chromium and NovaSeq work flow. The cutoff used in this study was 2500 UMI per cell,resulting in mean gene per cell ranging from 2935-3079. For downstream analysis 3004 cells for LVGFP,2063 cells for LV-AEP, and 3909 cells for LV-NEP were used. (c) Violin plots show percentage ofmitochondrial encoded genes detected in single cell RNA-seq analyses using 10X Genomics Chromiumand NovaSeq work flow. The low rates (<0.03%) of transcripts from the mitochondrial genome indicatethat the transcripts analyzed in this study are from the nuclear genome.

Supplementary Figure 5 | Expression of known genes used to assign cell categories of cell clusters. Feature plots of known genes in LV-GFP, LV-AEP, or LV-NEP infected retinal organoids shown as UMAPs.Genes representing different cell category or states are used to assign cell cluster identities in Figure 5.

Supplementary Figure 6 | Expression of featured known genes in combined sample clusters. Feature plots of known genes in the combined LV-GFP, LV-AEP, LV-NEP sample clusters shown asUMAPs. Genes representing different cell cycle phases and cell type categories are used to assigncluster identities in Figure 8.

Supplementary Figure 7 | Expression patterns of Notch signaling components and selected bHLH genes inlentiviral infected retinal organoid cells. Feature plots of selected genes quantified in Figure 9 are shown in UMAPs in LV-GFP, LV-AEP, and LVNEPinfected retinal organoid cells. HES1 Is predominantly expressed among neural stem cells andprogenitors. HES6, DLL1, and NEUROD1 are upregulated in exiting progenitors and neuroblasts. NHLH1is expressed among postmitotic cells.

Supplementary Figure 8 | Expression patterns of genes involved in cell cycle exit and early retinogenesis. Feature plots of selected genes involved in cell cycle exit (PLK1, CDKN2C) and differentiation of earlyretinal neurons (SOX4, ONECUT1, OTX2) are shown in UMAPs in LV-GFP, LV-AEP, and LV-NEP infectedretinal organoid cells.

Supplementary Table 1 | DEGs in Different Cell Categories of LV-GFP infected retinal organoids.

Supplementary Table 2 | DEGs in Different Cell Categories of LV-AEP infected retinal organoids.

Supplementary Table 3 | DEGs in Different Cell Categories of LV-NEP infected retinal organoids.

Supplementary Table 4 | Top 10 DEGs in LV-GFP infected Different Cell Categories.

Supplementary Table 5 | Top 10 DEGs in LV-AEP infected Different Cell Categories.

Supplementary Table 6 | Top 10 DEGs in LV-NEP infected Different Cell Categories.

Footnotes

References

Ahlmann-Eltze, C. (2019). ggsignif: Significane Brackets for ‘ggplot2’, R Package Version 0.6.0 ed.

Google Scholar

Akagi, T., Inoue, T., Miyoshi, G., Bessho, Y., Takahashi, M., Lee, J. E., et al. (2004). Requirement of multiple basic helix-loop-helix genes for retinal neuronal subtype specification. J. Biol. Chem. 279, 28492–28498. doi: 10.1074/jbc.m400871200