Mapping developmental trajectories and subtype diversity of normal and glaucomatous human retinal ganglion cells by single-cell transcriptome analysis.

Glaucoma is characterized by a progressive degeneration of retinal ganglion cells (RGCs), leading to irreversible vision loss. Currently, there is no effective treatment for RGC degeneration. We used a disease-in-a-dish stem cell model to examine the developmental susceptibility of RGCs to glaucomatous degeneration, which may inform on the formulation of therapeutic approaches. Here, we used single-cell transcriptome analysis of SIX6 risk allele (SIX6risk allele ) primary open angle glaucoma patient-specific and control hRGCs to compare developmental trajectories in terms of lineage- and stage-specific transcriptional signature to identify dysregulated stages/genes, and subtype composition to estimate the relative vulnerability of RGCs to degeneration because their ability to regenerate axons are subtype-specific. The developmental trajectories, beginning from neural stem cells to RGCs, were similar between SIX6risk allele and control RGCs. However, the differentiation of SIX6risk allele RGCs was relatively stalled at the retinal progenitor cell stage, compromising the acquisition of mature phenotype and subtype composition, compared with controls, which was likely due to dysregulated mTOR and Notch signaling pathways. Furthermore, SIX6risk allele RGCs, as compared with controls, expressed fewer genes corresponding to RGC subtypes that are preferentially resistant to degeneration. The immature phenotype of SIX6risk allele RGCs with underrepresented degeneration-resistant subtypes may make them vulnerable to glaucomatous degeneration.

Recent advances in single‐cell transcriptomics are paving the way to a comprehensive understanding of disease modeling in terms of cellular complexity, and dysregulated genes and signaling pathways. Application of this approach to the generation of retinal ganglion cells (RGCs) from SIX6 glaucoma patient‐specific and healthy control induced pluripotent stem cells revealed a flawed developmental trajectory in the former with immature and deficient subtype specification, likely due to dysregulated mTOR and Notch signaling pathways. The observations of this study shed light on the fidelity of RGC generation in vitro and influence of the primary open angle glaucoma risk allele on RGC development and subtype specification that may make RGCs susceptible to glaucomatous degeneration.

INTRODUCTION

Glaucoma is a complex group of diseases with multiple risk factors and genetic variants, in which a selective degeneration of the output retinal neurons, the retinal ganglion cells (RGCs), leads to irreversible blindness. , The mechanism underlying RGC degeneration is poorly understood, thus its treatment options remain limited to pharmacological or surgical mitigation of intraocular pressure, associated with primary open angle glaucoma (POAG). Given this intractable situation, stem cell modeling of glaucomatous degeneration may shed light on underlying pathology for the formulation of therapeutic approaches. In the last decade, significant progress has been made toward modeling glaucoma using pluripotent stem cell technology. For example, RGCs have been directly generated from human embryonic stem/iPS cells by default , or by stage‐specific recruitment of development mechanisms in two‐dimensional (2D) culture. The reproducible generation of hRGCs from iPS cells led to the development of a (a) disease model for POAG associated with the missense variant (rs33912345; C > A; His141Asn) in the exon of SIX6 (SIX6 ), , , for testing the hypothesis that RGCs are developmentally susceptible to glaucomatous degeneration, and (b) optic nerve regeneration model that demonstrated the supportive role of the mTOR signaling in the regeneration of hRGC axons following chemical axotomy. However, both these models can be further improved by characterizing (a) the developmental trajectories of control and disease‐specific RGCs in terms of lineage‐ and stage‐specific transcriptional signatures to give insight into dysregulated genes and pathways , that can be targeted therapeutically and (b) the range of subtypes generated to estimate the relative vulnerability of RGCs to degeneration, because their ability to regenerate optic nerve is subtype‐specific. , Here, to examine further the fidelity of modeling glaucoma and optic nerve regeneration in terms of cellular and molecular complexity, we have used a comparative single‐cell transcriptome analysis of hRGCs, generated from normal (controls) and SIX6 iPS cells. We observed that the developmental trajectories, defined by lineage‐ and stage‐specific transcripts, were similar for normal and SIX6 hRGCs. However, the development of SIX6 hRGCs appeared relatively stalled at the postmitotic precursor stage, resulting into fewer RGCs. These RGCs were immature compared with controls, as demonstrated by reduced expression of genes involved in RGC development and maturation. Additionally, SIX6 RGCs demonstrated expression of fewer RGC subtype‐specific genes, compared with controls, particularly of those that confer resistance to RGC degeneration. A comparative analysis of differentially expressed genes (DEGs) mapped on signaling pathways suggested that the immature phenotype of SIX6 RGCs, in which subtypes resistant to degeneration are underrepresented, is due to dysregulated mTOR and Notch signaling pathways in SIX6 RGCs. In summary, the developmental trajectories of SIX6 and control RGCs are similar, but the effect of the risk allele is manifested during the differentiation of retinal progenitor cells (RPCs) and/or/retinal precursors in the former resulting in RGCs, which are fewer, immature, and deficient in subtypes that are resistant to degeneration. This immature phenotype, if carried into adulthood, may make RGCs vulnerable to glaucomatous degeneration.

MATERIALS AND METHODS

Ethics statement

The study protocols were approved by the University of Nebraska Medical Center Institutional Review Board/Embryonic Stem Cell Research Oversight committees (and 397‐14‐EP/263‐10). Participants for this study were recruited after signing informed consent at Duke Eye Center, North Carolina, under the Duke University Health System protocols (CR006_Pro00009605).

Control and cells

We used previously established human iPS cell lines reprogrammed from peripheral blood mononuclear cells (PBMCs) isolated from blood samples from SIX6 POAG patients and age‐ and sex‐matched control donors by ectopic expression of KLF4, OCT4, SOX2, and C‐MYC, using a nonintegrating approach. Individual colonies of prospective induced pluripotent stem cells (iPSCs) were manually picked and clonally expanded to yield lines of control and patient‐specific iPSCs. Individual colonies were maintained in feeder‐free conditions using mTeSR1 medium (STEMCELL Technologies, Canada, https://www.stemcell.com/) in Matrigel‐coated dishes. Cells were passaged at 70%‐80% confluence at a 1:3 to 1:4 split ratio every 5 days using cell dissociation reagent (STEMCELL Technologies). All experiments were carried out using cells from passages 12 to 30. Cells were characterized for expression of pluripotency markers before RGC differentiation.

Retinal induction and RGC differentiation

Differentiation of human induced pluripotent stem cells (hiPSCs) along the retinal lineage and subsequent RGC differentiation was performed using the previously published protocol. Briefly, approximately 30 individual clumps of undifferentiated hiPSCs were plated onto Matrigel‐coated dishes in retinal differentiation medium consisting of DMEM/F12 B27 supplement (2%), N2 supplement (1%) in the presence of Noggin (1 ng/mL), human recombinant Dkk‐1 (1 ng/mL) and human recombinant insulin‐like growth factor‐1 (1 ng/mL) for 3 days. The concentration of factors was then increased to 10 ng/mL for up to 3 weeks. The resulting neural rosettes (NRs) were manually picked and grown as cellular aggregates for differentiation into RGCs. For RGC differentiation, NRs were manually picked and plated onto Matrigel‐coated dishes and RGC differentiation was initiated by treating cells for 2 days with Shh (250 ng/mL), FGF8 (100 ng/mL), and DAPT (3 μM). RGC differentiation was facilitated by treatment with Follistatin (100 ng/mL), Shh (250 ng/mL), and DAPT (3 μM) for 3 days. Finally, RGC maturation and survival were promoted by supplementing the medium with BDNF (100 ng/mL), Forskolin (10 μM), NT4 (5 ng/mL), CNTF (10 ng/mL) cAMP (400 μM), Y27632 (10 μM), and DAPT (3 μM) for the next 10 days. Medium was changed every 2‐3 days. All reagents were purchased from R&D systems (R&D Systems, Minneapolis, Minnesota, https://www.rndsystems.com).

Immunocytochemical analysis

For immunocytochemical analysis, cells at day 15 of differentiation were fixed in 4% paraformaldehyde for 15 minutes. After two washes with phosphate‐buffered saline (PBS), fixed cells were exposed to 5% normal goat/donkey serum in PBS for 30 minutes at room temperature, were permeablized with Triton X‐100 (0.4% or 0.2% for nuclear or cytoplasmic staining, respectively), followed by an overnight incubation with primary antibody at 4°C. Cells were incubated with fluorescence (Cy3/FITC)‐tagged secondary antibodies (Life Technologies, Carlsbad, California), diluted in PBS with 10% normal goat serum, for 2 hours at room temperature. DAPI (4′,6‐diamidino‐2‐phenylindole) was used to counterstain nuclei. Samples were mounted using VectaShield (Vector Laboratories, Burlingame, California) and fluorescent images were acquired with the Zeiss ApoTome Imager M2 upright microscope (Axiovert 200M) and the Axiovision 4.8 software (Carl Zeiss, GmbH, Germany, http://www.zeiss.co.in). A list of antibodies and working dilutions is provided in supplemental online Table 1.

3′ Droplet‐based scRNA sequencing

Single‐cell RNAseq libraries were prepared using the 10x Genomics Single Cell 3' Gel Bead and Library Kit following the manufacturer's instruction. Briefly, the single‐cell suspensions were diluted to the working concentration required to target 8000 number of cells captured. Using the Chromium Controller (10x Genomics), single cells were partitioned into gel beads emulsion containing unique barcoded primers with unique molecular identifier (UMI), followed by lysis of cells and barcoded reverse transcription of RNA, amplification of barcoded cDNA, fragmentation of cDNA to ∼200 bp, 5' adapter attachment and sample indexing. Libraries were quantified using the KAPA Library Quantification Kits for Illumina platform (KAPA Biosystems). The processed libraries were sequenced on an Illumina NextSeq550 instrument using 26 × 8 × 101 sequencing.

3′ Droplet‐based scRNA‐sequencing data analysis

Cellranger (10x Genomics) analysis pipeline was used for demultiplexing the scRNA‐seq output for aligning reads and gene expression analysis. Briefly, Cellranger mkfastq was used to convert the raw base call (BCL) files of scRNA‐seq into FASTQ files, which were then used for alignment, filtering, and processing barcode data for clustering and gene expression analysis through the cellranger count. Specifically, the first 10 principal components analysis was applied to the data to change the dimensionality to a user‐selectable number of principal components, followed by t‐stochastic neighbor embedding (t‐SNE) to visualize the data in a 2D space and graph‐based clustering method to cluster the cells. We used Loupe cell browser (Loupe cell browser; https://support.10xgenomics.com/single‐cell‐gene‐expression/software/visualization/latest/what‐is‐loupe‐cell‐browser), R packages including Cellranger R‐kit (http://cf.10xgenomics.com/supp/cell‐exp/cellrangerrkit‐PBMC‐vignette‐knitr‐2.0.0.pdf), complex heatmap, and Geom_violin for more in‐depth analysis. , The DEGs between groups of cells were identified using the Cell Ranger default differential express gene function algorithm. For each cluster, the algorithm is run on that cluster vs all other cells, yielding a list of genes that are differentially expressed in that cluster relative to the rest of the sample. A fold change cutoff of 0.5 and P‐value of .05 were applied to identify DEGs in each cluster for Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis. KEGG enrichment analysis was done using g:Profiler, and the Bonferroni correction method was used for multiple test correction.

Quality control pre‐ and postsequencing

For quality control, before sequencing the cell suspension was filtered through a fine‐mesh cell strainer to remove cell aggregates and then live/dead assay (Thermo Fisher Scientific) was performed to measure the cell viability. Postsequencing, during the analysis the captured cells were filtered by their UMI count. Using the Barcode Rank Plot in 10x genomics, which represents the distribution of barcode counts and which barcodes were inferred to be associated with cells.

RESULTS

Transcriptional identity of subtypes of cells (cell clusters)

We directly differentiated human iPS cells, reprogrammed from healthy nonglaucomatous donor blood, using the 2D chemically defined method as previously described. , To characterize developmental trajectory and molecularly distinct cell types, scRNA‐seq (10x Genomics) was performed on iPSC‐derived neural progenitors at day 15 of RGC differentiation (Figure 1A). Cells were collected, dissociated, and processed through 10x genomics chromium controllers to create single‐cell libraries from approximately 8000 cells, followed by sequencing barcoded libraries on Illumina platform. Sequencing data were analyzed using the single‐cell software, Cell Ranger (10x Genomics) to align reads, perform cluster and gene expression analyses, and for t‐SNE projection. An average of 70 000 reads were generated for each cell, identifying on average 3000 genes per cell. Next, we performed the t‐SNE analysis to explore the diversity of these cells, which identified nine distinct cell clusters (C1–C9), clearly separated in the t‐SNE plots (Figure 1B). Individual cell clusters possessed variable number of cells, ranging from 1241 to 419 cells per cluster (Figure 1B). Each cluster was defined by a unique set of DEGs, as illustrated in the heat map (Figure 1C). Next, in order to determine the identity of clusters along different developmental stages, we performed functional enrichment analysis on DEGs. We observed 27 significantly enriched gene ontology (GO) terms across the clusters. Based on the common GO terms, the nine clusters were categorized into three discrete groups corresponding to specific biological processes, namely, “Mitotic Cell Cycle” (C1 and C4), “Regulation of Neurogenesis” (C2, C3, C8) and “Eye Development” (C5, C7, C9), depicted in the heat map of the DEG (Figure 1C), t‐SNE plots of cluster‐specific genes (supplemental online Figure 1), and violin plots of the levels of expression of selected genes across nine clusters (Figure 2A). In these analyses, cell cluster 6 (C6) emerged as a separate group, characterized by DEGs corresponding to biological processes, which included Cerebral Cortex Neuron Differentiation, Forebrain Neuron Differentiation, and Cellular Developmental processes. This was a nonretinal group, expressing nonretinal genes such as HOXB4 minus the expression of retinal cell type‐specific genes. Of all cell clusters, C3 and C8 represented differentiating RGCs, given the expression of POU4F2, encoding a homeodomain TF regulating the differentiation of RGCs , (see below).

FIGURE 1

Cell type classification of human induced pluripotent stem cell (hiPSC)‐derived retinal ganglion cells (RGCs) by scRNA‐seq. A, Schematic diagram displaying the differentiation of RGCs from hiPSC through a stage‐specific chemically defined approach. At day 15 of differentiation, hiPSC‐derived cells were subjected to single‐cell RNAseq analysis. B, t‐distributed stochastic neighbor embedding (t‐SNE) plot showing the overall gene expression relationship among the 6652 single cells with ∼3000 genes detected in each cell. Different cell clusters are color‐coded. C, Heatmap displays the top 30 differentially expressed genes that were cell type‐specific in different clusters (intensity displayed from red [high] to blue [low]). The enriched biological processes for each gene group were shown in the right column. Classical cell type‐specific genes were labeled

FIGURE 2

Generation of a developmental trajectory of human induced pluripotent stem cell (hiPSC)‐derived retinal ganglion cells (RGCs). A, Violin plot showing the expression of cell type‐specific genes across the nine cell clusters. Each cluster is color‐coded. The y‐axis indicates normalized expression value, Log2(expected count +1). NES and SOX1 = neural progenitor‐specific genes; CCNB1 and UBE2C = mitotic cell‐specific genes; PAX6, SFRP2, SIX6, RAX, and VSX2 = RPC‐specific genes; POU4F2, DLX2, NEUROD2, and HOXB4 = neuronal genes. B, A developmental trajectory was constructed consisting of a major (retinal) and a minor (nonretinal) branch, based on cell type‐specific differentially expressed genes. C, The proportion of specific cell types at different stages in the two branches of the developmental trajectory is depicted in the graph

Developmental relationship between cell clusters

To know whether or not these cell clusters represented discrete dynamic stages of RGC development we characterized them into neural progenitors, precursors, and neurons based on DEGs (Figure 2B,C). We observed that cell clusters could be broadly segregated into two groups: (a) progenitors/precursors (C1, C2, C4, C5, C7, C9), expressing neural progenitor‐specific genes such as NES (Nestin), PAX6, SFRP2, SIX6, RX, and VSX2, and (b) neurons (C3, C6, C8), expressing genes corresponding to cell type‐specific regulators such as POU4F2, DLX2, NEUROD2, and HOXB4. The cell clusters in the progenitor/precursor group were further segregated into proliferating progenitors (C1 and C4) expressing mitotic genes (eg, MKI67, CCNB1, NUSAP1, TOP2A, TPX2, CENPF, SMC4, HMGB2, HIST1H4C, and UBE2C) and postmitotic precursors (C5, C7, C9), which were not expressing mitotic genes listed above. In the proliferating progenitor group, cluster C1 cells appeared developmentally more primitive than cluster C4 cells because they expressed neural stem cell marker NES, whereas the latter appeared to have progressed along the retinal lineage, as suggested by the expression of the canonical RPC‐specific genes such as PAX6, SFRP2, SIX6, RX, and VSX2. Therefore, cluster C1 represented developmentally the most primitive of all cell clusters, and thus occupied top of the hierarchy in the developmental trajectory. From cluster C1 emerged two different developmental trajectories; a retinal trajectory (∼58% of total cells), consisting of cells expressing RPC and/or RGC markers, and a nonretinal trajectory (∼23% of total cells), with cells absent in RPC and retinal cell‐specific marker genes. Cluster C4, consisting of proliferating cells expressing RPC‐specific genes, led the retinal trajectory. Cell clusters (C2, C5, C7, C9) in postmitotic precursor group, with the exception of cluster C2, expressed canonical or recently identified (SFRP2/PENK/DKK3) RPC marker genes, represented the intermediate steps toward RGC differentiation in the retinal trajectory. However, cells in cluster C2 were enriched in the expression of NES and did not express RPC‐specific genes; therefore, they were assigned to the nonretinal trajectory. Cell cluster C6, expressing HOXB4, EMX2, and ASCL1, belonged to the nonretinal trajectory.

Characteristics of RGC‐specific cell clusters and subtype diversity

Cell clusters C3 and C8, which belonged to the neuron group in the retinal trajectory, were characterized by the shared expression of RGC‐specific regulators, ATOH7, POU4F2, and ISL1 and other RGC‐specific genes , , , , (Figure 3A,B). However, a closer examination of DEGs shed light on the developmental relationship between these cell clusters. For example, cluster C8 appeared to be developmentally primitive than cluster C3 because, although it expressed RGC‐specific genes, it maintained robust expression of RPC markers, RX, SFRP2, and SIX6, compared with cluster C3 (Figure 3B). Additionally, the expression levels of the majority of RGC‐specific genes, particularly those involved in RGC maturation, were significantly lower in cluster C8 vs cluster C3, further demonstrating the relative immaturity of the former (Figure 3A). Twice as more cells in cluster C3 (∼15%) belonged to the category of matured RGCs than in cluster C8 (∼7%). To determine the fidelity of directed RGC differentiation, we examined cluster C3 cells for the DEGs, corresponding to early and late born retinal cells. Expressions of other cell type‐specific genes such as PROX1 (amacrine/horizontal cells), VSX1 (bipolar cells), and SLC1A3 (Müller glia) were observed at barely detectable levels in cluster C8 and cluster C3, demonstrating the dominance of RGC‐specific transcriptional activity in cells belonging to the RGC sublineage (supplemental online Figure 2). Next, given the emerging evidence that the function of RGCs, their susceptibility to degeneration, and their capacity to regenerate axons are subtype‐dependent, , we examined cluster C3 cells for their RGC subtype characteristics, primarily following the Sanes and Masland classification, based on physiological, morphological and molecular criteria, supplemented with subtype‐associated genes, recently identified through the single‐cell transcriptome analysis , (Figure 4A). Differential expression of subtype‐specific genes suggested representation of all four broad functional subtypes of RGCs, namely the ON RGCs, OFF RGCs, ON/OFF DS RGCs, and IP‐RGCs (Figure 4B). These genes were expressed in cluster C3 cells and not in the nonretinal cluster C6 (data not shown), demonstrating the association of subtype‐specific gene expression specifically with RGCs. The expression of these genes was significantly less in cluster C8 cells vs cluster C3 cells, further attesting to a relative immature state of the former and demonstrating that acquisition of the subtype specificity may be a part of the maturation process (Figure 4B). Next, we carried out immunocytochemical analysis of hRGCs and observed POU4F2+ cells expressing immunoreactivities corresponding to four basic subtypes of RGCs, demonstrating that the subtypes classified on the basis of DEGs may represent different subpopulations of RGCs (Figure 4C‐G).

FIGURE 3

Expression profile of key genes related to retinal ganglion cell (RGC) specification, differentiation, and maturation in RGC‐specific cell clusters. A, Heatmap revealing the scaled expression of RGC‐specific genes in RGC‐specific cell clusters. Green to red represents low to high gene expression. B, Violin plot showing the expression profile of representative genes in the RGC‐specific cell clusters. Different clusters are color‐coded

FIGURE 4

Identification of different subtypes in retinal ganglion cell (RGC)‐specific cell clusters. A, The schematic representation of RGC subtypes classification based on function and representative genes. B, Heatmap showing expression profile of genes corresponding to ON/OFF DS RGCs, ON DS RGCs, alpha‐RGCs, and IP‐RGCs subtypes in control clusters C3 and C8. C, Representative merged immunofluorescence images demonstrating cells coexpressing POU4F2 and immunoreactivities corresponding to DCX (ON/OFF RGCs), D, FSTL4 (ON DS RGCs), E, SMI32 (alpha‐RGCs), F, CALB2 (transient alpha‐RGCs), and G, TBR (IP‐RGCs) when human induced pluripotent stem cells were directly differentiated along the RGC lineage. Scale bar: 50 μm

Developmental trajectories of cell clusters obtained from cells

We have recently demonstrated that RGCs generated from SIX6 POAG patient‐specific iPS cells have developmental abnormalities that result in shorter and less complex neurites and immature electrophysiological signatures, accompanied by dysregulated genes and pathways, compared with those generated from age‐ and sex‐matched healthy donor iPS cells, used above (supplemental online Figure 3A). Here, we used single‐cell transcriptome analysis to determine the effect of the SIX6 on the developmental trajectory of RGCs and whether or not it influenced the RGC subtype diversity, which might influence RGC survival and regeneration. Single‐cell RNAseq analysis on SIX6 RGCs was carried out as described above for the age‐ and sex‐matched healthy donor control. Fourteen cell clusters (C1‐C14) were revealed, clearly separated in the t‐SNE plots (supplemental online Figure 3B). Each cluster was defined by unique DEGs and classified into specific groups, corresponding biological processes and subsequently into specific developmental stages as described above (supplemental online Figure 3C). Cluster C13 did not pass the P‐value threshold and failed to enrich DEGs relative to every other cluster, and was excluded from the analysis. The developmental trajectory of SIX6 RGCs consisted of retinal and nonretinal branches and could be superimposed on the developmental stages determined for nine clusters of cells derived from the control iPSCs (Figure 5A). However, the following differences were observed (Figure 5B): (a) the majority of cells (61%) in SIX6 RGC trajectory belonged to the proliferating RPC and postmitotic precursor groups as compared with 36% of such cells in controls; (b) only 4% of cells (C11) constituted the RGC group, compared with 22% in controls (C3 + C8), suggesting that the efficiency of RGC generation was rendered deficient in the SIX6 iPS cells, and the road block in the generation of RGCs appeared to be at the RPC/retinal precursor stages, particularly at the postmitotic precursor stage which consisted of five clusters (C1, C5, C7, C12, and C14) and (c) cells in the nonretinal branch (C2and C8), although not expressing RPC or retinal cell type‐specific genes, were diverse in the expression of genes compared with controls, such as LMO3, encoding a LIM only family of protein, , and GRID2, encoding a glutamate ionotropic receptor.

FIGURE 5

Generation of a developmental trajectory of SIX6 induced pluripotent stem cell (iPSC)‐derived retinal ganglion cells (RGCs). A, A developmental trajectory was constructed consisting of a major (retinal) and a minor (nonretinal) branch, based on cell type‐specific differentially expressed genes. B, Graph depicting a comparative profile of the proportion of cells in each cell types across all the clusters in both control (Figure 3) and SIX6 RGCs

A comparison of DEGs in cluster C11 and control cluster C3 cells revealed that although the former expressed key regulators of RGC specification (eg, ATOH7, POU4F2, and ISL1) it was deficient in the expression of regulators of RGC maturation (eg, DLX1 and DLX2), axonogenesis (eg, KLF6, SOX4, and SOX11), and axonal guidance (eg, ROBO2, GAP43, and DCC), compared with the latter (Figure 6), suggesting a risk allele‐associated molecular barrier toward maturity. The relative immaturity of cluster C11 cells suggested that their maturation into different subtypes would be compromised. To test this premise, we carried out subtype characterization of cluster C11 cells as described for control cluster C3 cells. DEGs in cluster C11 included those that were associated with ON/OFF DSCS, ON‐DSGCs, alpha RGC, and ipRGCs. However, as compared with control C3 cluster, C11 cells expressed fewer genes characteristic of different RGC subtypes (Figure 7A,B). Since DEGs in each cluster represented the cumulative expression of genes in all cells in that cluster, we examined the expression of subtype specific genes in subclusters to know if their expression remained significant at the single‐cell levels and thus identified a cell or a group of cells belonging to a specific subtype. This analysis revealed six subclusters in control cluster C3 (supplemental online Figure 4), and five subclusters in the disease‐specific cluster C11 (supplemental online Figure 5). Four subclusters in C3 (SC1, SC3, SC4, and SC6) expressed different classical subtype‐specific genes, in which two subclusters (SC2 and SC5) remained unclassified but expressed recently identified subtype‐specific transcription factor genes (supplemental online Figure 4A,B). In contrast, only three C11 subclusters (SC1, SC2, and SC4 ) could be subtyped and each expressed fewer subtype‐specific genes, compared with controls (supplemental online Figure 5A,B). SC3 and SC1 remained unclassified for not expressing classical subtype specific genes. Next, we tested the premise that SIX6 RGCs may have molecular and cellular attributes that may shed light on glaucomatous susceptibility, which in the SIX6 POAG patient is demonstrated by thinned retinal nerve fiber layer (RNFL). , Examination of the genes recently identified with the RGC subtypes, which are preferentially resistant (eg, OFF‐Sa, ON‐a, IP‐RGCs) or susceptible (eg, OFF‐Ta and ON‐OFF DSGCs) to optic nerve injury, revealed that these were expressed in both C3 and C11 cells (Figure 7C). However, we observed that the number and expression of resistant‐ and susceptibility‐associated genes were decreased and increased in C11 cells vs C3 control cells, respectively. Next, we wanted to understand the putative mechanism underlying SIX6 ‐associated delay in the developmental trajectory and subtype‐specific differentiation of RGCs. We mapped DEGs on the KEGG pathway to identify developmentally relevant signaling pathways, adversely affected by the SIX6 , compared with controls (Figure 7D, E). We observed that besides the cell cycle checkpoint pathway, active in proliferating cells, the Notch pathway, whose inhibition is necessary for RPCs to differentiate into RGCs, , was increased in SIX6 RGCs, compared with controls. In contrast, the mTOR pathway, which facilitates RGC differentiation, was decreased in SIX RGCs vs controls. Together, these observations suggested that the in vitro generated SIX6 RGCs were developmentally compromised due to the stalling of the development at the levels of RPC/retinal precursors, presumably due to dysregulation of developmentally relevant signaling pathways. In addition, SIX6 RGCs were likely to be more susceptible to injury compared with controls, due to deficient expression of degeneration resistant subtype‐specific genes.

FIGURE 6

Comparison of differentially expressed genes between control and SIX6 retinal ganglion cells (RGCs). Heatmap representing a comparative expression profile of key genes corresponding to eye field transcription factors (EFTFs), RGC specification, differentiation, RGC maturation and axon guidance, synaptogenesis, and axonogenesis in the RGC‐specific clusters in control vs SIX6 groups

FIGURE 7

Comparative expression profile of retinal ganglion cell (RGC) subtype, and resilient/susceptible RGCs specific genes in control and SIX6 RGCs. A, Heatmap showing expression of fewer RGC subtype‐specific genes in SIX6 RGCs compared with controls. B, Representative merged immunofluorescence images demonstrating SIX6 RGCs coexpressing POU4F2 and immunoreactivities corresponding to DCX (ON/OFF RGCs), FSTL4 (ON DS RGCs), SMI32 (alpha‐RGCs), CALB2 (transient alpha‐RGCs), and TBR2 (IP‐RGCs). C, Heatmap demonstrating high expression level of susceptible and low expression level of resilient RGC‐specific markers in SIX6 compared with control RGCs. D, Kyoto Encyclopedia of Genes and Genomes pathway analysis demonstrating pathways mapped on up regulated (mTOR signaling pathway, Axon guidance) and down regulated (Notch signaling pathway, Cell cycle checkpoints) genes in control RGC‐specific cell cluster, and E, representative pathways mapped on upregulated (cell cycle checkpoints, Notch signaling pathway) and downregulated (PI3K‐Akt‐mTOR signaling, G2/M transition) genes in SIX6 RGC‐specific cell cluster. Scale bar: 50 μm

DISCUSSION

The stem cell approach to glaucoma, whether to replace degenerating RGCs or establish a disease‐in‐a‐dish model of glaucomatous RGC degeneration, requires ex vivo generation of hRGCs from pluripotent stem cells with high efficiency and fidelity. This includes generation of hRGCs through normal developmental trajectory containing complementation of different subtypes, against which the developmental aspects of RGC abnormality in a disease model can be evaluated. Information about subtype characteristics is not only important from functional viewpoints but also for understanding the underlying mechanism of glaucomatous degeneration given the emerging evidence that the susceptibility and resistance of RGCs are subtype‐dependent. ,

Developmental trajectories

Segregation of populations of cells based on unique combinations of DEGs revealed cell clusters that identified two dynamic branches, beginning from the generic neural stem cell population; (a) retinal branch consisting of RPCs and RGCs, and (b) nonretinal branch giving rise to nonretinal neurons. The majority of control cells (∼58%) were on the trajectory of retinal lineage (Figure 2B) demonstrating the preferential influence of the culture conditions to promote retinal differentiation. However, of those cells on the retinal trajectory, ∼24% and ∼20% were in postmitotic retinal precursor stage and mature RGC stages, respectively. The developmental trajectory suggests that the efficiency of RGC generation could be increased by facilitating the transition of cells from the postmitotic retinal precursor stage to RGCs either by changing the culture conditions or increasing the culture time or both. Preliminary observations that extending the culture time substantially increased the number of POU4F2+ cells (Teotia and Ahmad, unpublished results) suggest that the majority of cells along the retinal lineage when culture was terminated were in transition toward RGC differentiation. The developmental trajectory of SIX6 RGCs was similar to that of control RGCs, however, twice as many cells as in controls were held at the RPC stage (Figure 5B). Cells in the postmitotic retinal precursor stage were more than 1.5‐fold higher in SIX6 RGCs culture than in controls. In contrast, cells in the RGC stage were 5‐fold less in the SIX6 RGCs culture than in controls. Together, these observations suggest that the presence of the risk allele impeded the differentiation of RPCs along RGC lineage and those that differentiated into RGCs did not express genes that are important for RGC maturation, including those that regulate neuritogenesis and axon guidance. This corroborates the earlier finding of immature phenotype of SIX6 RGCs in terms of simpler neurites and electrophysiological profiles, compared with controls. The risk allele may adversely influence the mechanisms that underpin RGC differentiation. Mapping of DEGs on KEGG pathways identified two signaling mechanisms, the Notch and mTOR pathways, whose dysregulation can profoundly influence the transition of RPCs from one stage to another and their differentiation into RGCs. Notch signaling is a gate‐way to differentiation; its inhibition is necessary for the progression of uncommitted progenitors along a specific lineage. , , , Inhibition of Notch pathway leads to the activation of ATOH7, a first essential step for committing the progenitors on the RGC lineage. The mTOR pathway positively regulates RGC differentiation and evidence suggests that this may involve interactions with the Notch pathway. The increase and decrease in Notch and mTOR pathways in C11 SIX6 cluster, compared with control C3 cluster, respectively, suggest that in the condition of persistent Notch signaling the proper transitions of RPCs will be impeded, and with reduced mTOR signaling, differentiation and maturation of committed precursors will be hampered. This would lead to inefficient and immature differentiation of RGCs.

RGC subtype classifications

RGCs are heterogeneous in terms of their response to direction and intensity of light, size, mosaic distribution, pattern of dendritic arborization, and gene expression profile. Recent information of cell type‐specific gene expression through promoter‐reporter transgenes , and scRNA‐seq analysis, , superimposed over aforementioned classic criteria for RGC subtyping, aided by serial section electron microscopy and optical imaging of electrical activities suggests that there are in excess of 40 RGC subtypes. However, the majority of these studies were carried out in mice, therefore the RGC subtype classification in human remains poorly understood. A recent study carried out on RGCs derived from human iPSCs, using single‐cell qPCR for known RGC‐ and RGC subtype‐specific markers, demonstrated expression of transcripts corresponding to DS RGCs, alpha‐RGCs, and IP‐RGCs. However, this study was limited in capturing the subtype diversity because of limited sets of genes that were examined. The unbiased analysis of DEGs carried out here revealed that POU4F2 cluster (C3) reflected diversity corresponding to ON RGCs (ON DSRGCs and ON alpha‐RGCs), OFF RGCs, ON‐OFF DSRGCs and IP‐RGCs. This cluster was diverse in the expression of additional subtype‐specific genes recently identified by single‐cell transcription profiling of purified mouse RGCs. For example, analysis of these genes in subclusters revealed the association of TBR1, ZIC1, and ZFXH3 with SC1 (IP‐RGCs); ZFHX3 and TAGLN2 with SC3 (ON‐OFF DSRGCs and alpha‐RGCs); and NEUROD1, CRAB1, and TAGLN2 with SC4 (ON DS RGCS and transient OFF alpha‐RGCs). The two unclassified subclusters, SC2 and SC4, based on the absence of known subtype‐specific genes, were characterized by the expression of CRABP1 + ZIC1 and NEUROD1 + CRAB1 + ZHX3, respectively. That the acquisition of this diversity is progressive was revealed by low expression of fewer subtype‐specific genes in immature RGC cluster, C8. A similar immaturity in the expression of the range of subtype‐specific genes was observed in SIX6 cluster (C11Six6‐RA). Besides the range of responses to intensity and directions of incident light, the diversity of RGC subtypes also underpins the response of RGCs to their survival and optic nerve regeneration. Recent studies using the optic nerve crush (ONC) model in mice have revealed that alpha RGC/IP‐RGCs and ON‐OFF DSRGCs are relatively more resistant and susceptible, respectively, to ONC and that alpha‐RGCs preferentially regenerate following injury. , Resistant RGCs maintain normal dendritic morphology along with subtype‐specific light response following ONC. Examination of human RGC response to injury in a microfluidic model of chemical axotomy revealed that the regenerating axons belong to RGCs expressing markers corresponding to alpha‐RGCs, suggesting that the RGC subtype‐specific response to injury and regeneration is evolutionarily conserved. Recently, a single‐cell transcription profile of RGCs following ONC has identified DEGs preferentially associated with resistant and susceptible RGCs. Nine of these resistant associated genes (RAMP1, KBTBD11, EOMES, IFITM10, CHL1, SPP1, TBR1, ESRRG, and GLDN) were expressed in control RGCs vs two (ESRRG and GLDN) in SIX6 RGCs. This may reflect the developmental abnormality or immaturity of SIX6 RGCs vs controls. However, if the expression of genes that play a role in RGC survival and optic nerve regeneration remain deficient in the adult it may lead to glaucomatous pathology seen in SIX6 POAG patients with thinned RNFL. ,

CONCLUSION

In summary, single‐cell transcriptome analysis revealed the recapitulation of normal developmental trajectories during the generation of hRGCs, deviation from which, due to SIX6 ‐associated dysregulation of developmentally relevant signaling pathways, may lead to immature phenotype including fewer RGC subtypes. The deficient generation of those RGC subtypes that confer resistance to RGC degeneration may make SIX6 RGCs vulnerable to glaucomatous degeneration.

CONFLICT OF INTEREST

The authors declared no potential conflicts of interest.

AUTHOR CONTRIBUTIONS

P.T.: collection and/or assembly of data, data analysis and interpretation, and manuscript writing; M.N.: data analysis and interpretation; I.A.: conception and design, manuscript writing, financial support, and final approval of manuscript.

Supplementary Figure 1 Characterization of progenitor and neuronal cell types by scRNA‐seq analysis. A) t‐SNE plots showing expression of cell type‐specific genes in distinct cell clusters. The gene expression level is color‐coded (gray, no expression; red, increased expression)

Click here for additional data file.

Supplementary Figure 2 Specificity of RGC generation from hiPSCs. A) Heatmap showing barely detectable expression of select marker genes of early (amacrine, horizontal and photoreceptors) and late born (bipolar cells and Müller glia) retinal cell types, demonstrating the specificity of RGC differentiation, ascertained by scRNA‐seq analysis

Click here for additional data file.

Supplementary Figure 3 Cell type characterization of SIX6 RGCs by scRNA‐seq analyssis. A) Representative merged immunofluorescence images showing morphological differences between control RGCs and SIX6 RGCs coexpressing immunoreactivities corresponding to POU4F2 and βIII tubulin. B) t‐SNE plot showing the overall gene expression relationship among the 7177 single cells with ∼2982 genes detected in each cell derived from the SIX6 iPSCs. Different cell clusters are color‐coded. C) Heatmap displaying the DEGs that were cell type‐specific in different clusters. The enriched biological processes for each gene group were shown in the right. Classical cell type marker genes were labeled. The color key from blue to red indicates low to high gene expression respectively

Click here for additional data file.

Supplementary Figure 4 Subtype‐specific genes in control RGC sub clusters. A) RGC‐specific cell cluster C3 was further analyzed by sub‐clustering. Cluster 3 was divided into 6 distinct sub‐clusters, marked by expression of POU2F2, ZCCHC12, TBR1, ZIC1, ZFHX3, EOMES (IP‐RGC; SC1, SC6), COL25A1 (ON/OFF DS RGCs; SC3), SPP1 (alpha RGCs; SC3), FSTL4, NEUROD1, CRAB1, TAGLN2 (ON‐DS RGCS; SC4), CALB2 (Transient alpha RGCS; SC4) and unclassified (SC2 and SC5). B) Table demonstrating the specific RGC subtypes present within the sub‐clusters based on the expression of genes characteristic for each subtype according to Sanes et al, 2015 and Rheaume et al, 2019

Click here for additional data file.

Supplementary Figure 5 Subtype‐specific genes in SIX6 RGC sub clusters. A) RGC‐specific cell cluster C11 was further analyzed by sub‐clustering. Cluster C11 was divided into 5 distinct sub‐clusters, marked by expression of ADCYAP1, POU4F2, POU6F2, TAGLN2 (IP‐RGC; SC1 ), TRH, CDH6 (ON/OFF DS RGCs; SC2 , SC4 ), and unclassified (SC3 , SC5 ). B) Table demonstrating the specific RGC subtypes present within the sub‐clusters based on the expression of genes characteristic for each subtype according to Sanes et al, 2015 and Rheaume et al, 2019

Click here for additional data file.

Supplementary Table 1 List of Primary antibodies used for immunocytochemistry and Flow cytometry.

Click here for additional data file.