Dissecting the transcriptome landscape of the human fetal neural retina and retinal pigment epithelium by single-cell RNA-seq analysis

1During the development of the vertebrate embryo, the inside of the optic cup forms the neural retina (NR), whereas the outside becomes the retinal pigment epithelium (RPE). Multipotent retinal progenitor cells (RPCs) are located in the inner layer of the optic cup and give rise to essentially all cell classes in the retina, including six classes of neurons—namely, the output neurons: retinal ganglion cells (RGCs); interneurons such as horizontal cells (HCs), amacrine cells (ACs), and bipolar cells (BCs); and rod and cone photoreceptor cells (PCs), as well as glial cells and Müller glia cells. In mice, RGCs are generated at embryonic day (E)8, with a peak at approximately E11. ACs appear at E8, with a peak at E16. BCs are born at E17.

2Many basic helix-loop-helix (bHLH) transcription factors are expressed in RPCs and are involved in determining retinal cell fate. For example, Pax6, Sox2, and Vsx2 are transcription factors known to be involved in RPC multipotency. Atoh7 induces the expression of Pou4f2 and Isl1, which are important for the specification and differentiation of RGCs. Foxn4 expression in RPCs is required for AC and HC specification. Foxn4 activates the expression of Neurod1, Neurod4, and Ptf1a to confirm the AC fate. Rax2 and Otx2 are involved in the generation of PCs. Vsx2 and Ascl1 contribute to the specification of BCs. The overexpression of Hes1, Hes5, and Hey2 promotes the Müller glial cell fate. Although functional studies of the transcription factors mentioned above revealed solid clues for retinal cell fate determination, it is difficult to perform functional experiments on human fetal retina. Thus, the transcription factor profiles of the developing human retina, especially active transcription factors, need to be investigated and could offer helpful clues for the study of human retina development in retinal organoids.

3In humans, the analysis of a histological collection of human eyes identified the sequence of retinal layer development. Immunohistochemical staining for transcription factor expression in human embryos at 6–12 weeks (W) showed the division of the early NRs into the KI67 (a marker of proliferating cells)-positive outer layer and the Tuj1 (a neuronal marker)-positive inner layer. Additionally, the expression patterns of known eye-field transcription factors and related transcription factors in these two layers have also been revealed. Although the histological analysis of the developing human retina has helped researchers understand the structures and morphology of retinal cells, detailed molecular evidence for the development of human retina is lacking. Numerous studies of the human NR and RPE transcriptome revealed the gene expression in these two tissues. Hoshino and colleagues and Aldiri and colleagues analyzed bulk RNA sequencing (RNA-seq) data of human fetal retina. More specific expression profiles of each human fetal retinal cell class need to be investigated.

4Single-cell RNA-seq analysis could satisfy the need of higher specificity. Macosko and colleagues and Shekar and colleagues performed single-cell RNA-seq on mouse retina after birth using Drop-seq. They offered rich resources for mouse retina studies. There are several studies on the human retina preprinted at bioRxiv. Liang, Qingnan, and colleagues and Lukowski, Samuel, and colleagues completed the adult retina single-cell transcriptome atlas, and Lu, Yufeng, and colleagues provided a resource for the fetal human retina at the stages after 12 W. However, none of them discussed the development of fetal retinal cells at early stages (before 12 W) or the interaction of NR and RPE. Our work is the first to dissect the transcriptome landscape of NR and RPE before 12 W.

5We performed single-cell RNA-seq using a previously described modified single-cell tagged reverse transcription (STRT) method on individual cells collected from the human fetal NR and RPE at 5–24 W, following receipt of ethical approval, to address these questions. We identified the major cell classes in the NR and RPE, providing the single-cell transcriptome landscape of each cell class. Importantly, we dissected the temporal order of the generation of human retinal cells. Through a systematic bioinformatics analysis, we revealed important transcription factors with high activity in regulatory networks required for the development of the NR and RPE. Moreover, we revealed dynamic changes in the transcriptome features of RGCs and RPE cells. We explored the developmental period at which the visual cycle might occur based on the single-cell transcriptome features of PCs and RPE cells, as well as other interactions of these two cell classes. Furthermore, we mapped disease-related genes in our data set to explore potential targets of inherited retinal diseases. Thus, we offer a rich resource for studying the human NR and PRE. The expression patterns of single-cell RNA-seq data in the present study are available at http://49.4.93.68:30004/.

6We performed single-cell RNA-seq analysis of individual cells isolated from the human fetal NR and RPE ranging from 5 to 24 W (10 developmental stages) to obtain transcriptional maps of these two layers (Fig 1A). After quality control (gene number ≥ 1,000, unique molecular identifier [UMI] ≥ 10,000), 2,421 cells were retained, and an average of 4,631 expressed genes (transcripts per million [TPM] > 0) (S1 Fig) and 129,105 transcripts were detected in each individual cell (S1 Fig). We did not detect significant batch effects in these experiments (S1 Fig).

7We performed an unbiased principal component analysis (PCA) of these 2,421 cells using highly variable genes to examine the global transcriptome patterns in the single-cell RNA-seq data (Fig 1B). The first segregating factor was different tissues, as cells were divided into two groups consistent with the tissue layers from which they were collected. The positive genes in PC1 were enriched for gene ontology (GO) terms such as nervous system development (P value = 6.2E-09), and the negative genes in PC1 were enriched for GO terms such as retinol metabolic process (P value = 1.9E-04). The second segregating factor was the developmental stage. Cells were ordered in a time course from the early to the late stage. The positive genes in PC2 were related to visual perception (P value = 2.3E-12), and the negative genes were related to cell division (P value = 2.4E-10) (S1 Table).

8We reduced the 18 statistically significant principal components in the PCA to two dimensions using t-distributed stochastic neighbor embedding (t-SNE) to define the cell identities. These 2,421 cells were divided into 21 distinct clusters (Fig 1E). Cells from different layers were clearly segregated and ordered in a time-dependent manner (Fig 1C and 1D). Based on the specific expression of known cell type–specific markers for the NR and RPE (Fig 1F), we identified major cell classes of human fetal retina—namely, RPE cells expressing SERPINF1, TYR, MITF, RPE65, BEST1, and TTR; RPCs expressing SFRP2, MKI67, SOX2, HES5, and FZD5; RGCs expressing gamma synuclein (SNCG), NEFL, ATOH7, and EBF3; HCs expressing ONECUT1/2/3; ACs expressing MEIS2, GAD1, and GAD2; BCs expressing VSX1 and VSX2; PCs expressing PDC, PDE6G, SAG, CRX, and NRL; Müller glia expressing PLP1, TF, and RLBP1; microglia expressing CX3CR1, C1QA, C1QB, and C1QC; fibroblasts expressing COL3A1 and COL1A1; and a few blood cells expressing ALAS2 and HBM. Markers of each cluster and their GO analysis results were provided in S2 and S3 Tables.

9We further analyzed clusters containing various cell classes to obtain a more accurate identity for each single cell (S4 Table). For example, cluster 16 contained both VSX1 (BC marker)-positive and NRL (a rod transcription factor)-positive cells (Fig 1F). We performed t-SNE analysis on the cells in cluster 16 to define BCs and PCs in cluster 16 (S2 Fig). Cells in cluster 16 could be divided into two subclasses, one expressing BC-enriched genes CHN2 and SCN3A and the other expressing PC markers PDC and AIPL1. The same strategies were performed for cluster 10 (containing RPCs and Müller glia) (S2 Fig), cluster 4 (containing HCs and ACs), and cluster 3 (containing RPC subclasses) (S2 Fig). Subsequently, each single cell was assigned a potential identity for further analyses (Fig 2B, S5 Table). Cell class–specific markers were identified using the FindMarker function in Seurat (Fig 2A and S4 Table). According to the cell class–specific markers, we inferred cells in the PC group as rod photoreceptors because they expressed earlier markers of rod photoreceptors, such as NRL and ROM1. Besides, we found that NTRK1, which was reported to be expressed in injured RGCs, was one of the markers in human fetal HCs. We detected the expression of NTRK1 in ONECUT2 (a marker of HCs)-positive cells by RNAscope assay (S3 Fig).

10We detected four subclasses of RPCs (S2 Fig). The differentially expressed gene (DEG) analysis showed that one RPC subclass had relatively higher expression levels of RPC-enriched genes SFRP2, SHISA2, and CPAMD8. One subclass was highly proliferative, expressing TOP2A and MKI67. RGC markers, such as NEFM, PRPH, NEFL, ISL1, and GAP43, were expressed in one RPC subclass. We detected a small group of RPCs expressing C1QA, CD74, and HLA-DRA.

11The single-cell RNA-seq analysis was used to identify the different classes of cells in the developing NR and markers for each retinal cell class. A GO analysis of markers for each cell cluster indicated the potential states of each cell class (S2 and S3 Tables). For example, RPCs were grouped in cluster 3, and the GO analysis showed that these cells were enriched for cell division processes (P value = 1.8E-20). Clusters 2, 5, 8, and 12 were RGCs that primarily participated in nervous system and axon development. Clusters 13 and 4 were HCs. Markers of cluster 13 were related to neuron projection development (P value = 4.1E-05) and retinal layer formation (P value = 5.1E-05). Cluster 6 contained ACs that participate in neurotransmitter secretion (P value = 5.0E-12) and chemical synaptic transmission (P value = 4.3E-09). Markers of cluster 9 were significantly related to PC maintenance (P value = 4.9E-07). BCs in cluster 14 were involved in visual perception (P value = 2.8E-09), response to stimulus (P value = 1.1E-05), and retinal homeostasis (P value = 3.8E-05).

12We summarized the number of cells of each class at each stage based on the data in S5 Table to determine the time at which each cell class in the NR was generated (Fig 2C and 2D). RPCs and RGCs were already observed at 5 W. The proportion of RPCs peaked at 6 W, and the highest proportion of RGCs was observed at 8 W. Subsequently, HCs developed (reached a peak at 9 W), followed by ACs (reached a peak at 17 W), PCs, Müller glia cells, and BCs.

13Next, we analyzed the cell cycle features of each cell class we identified (S6 Table). At 5 and 6 W, fibroblasts, microglia, and RPE cells were actively proliferating (Fig 3A–3C, S4 Fig). RPCs in the NR were actively proliferating (Fig 3A, 3B and 3D, S4 Fig). We used KI67, the marker of proliferating cells, and the RPC-enriched gene SOX2 to evaluate the location and state of RPCs in the developing human retina (Fig 3E and 3F). KI67+ and SOX2+ cells were located in the outer layer of the retina. The density of KI67+ cells in the NR decreased as retinal development proceeded. We detected two classes of SOX2+ cells in the NR at 13 W. One SOX2+ cell class had a spindle-shaped nucleus resembling RPCs at 8 W and 9 W; this class of cell was located at the outer nuclear layer (ONL). The other SOX2+ cell class had a round-shaped nucleus and was located in the inner nuclear layer (INL). We postulated that the round SOX2+ cells in the INL were Müller glia cells.

14Transcription factors play important roles in eye development. To identify transcription factors with high activity, we applied a method called single-cell regulatory network inference and clustering (SCENIC), which evaluates the activity of the gene regulatory networks (GRNs) in each cell. The t-SNE analysis of the binary regulon activity matrix (S5 Fig) was consistent with the t-SNE analysis based on expression matrices in Fig 2B.

15In the present study, transcription factors that were active among more than 50% of cells in a particular cell class were retained. The transcription factors with high activities were selected by the AUCell algorithm (Fig 4A and 4B), and they were defined as active transcription factors. In the NR, RPCs expressed differentiation-related genes, including SOX2, RAX, and VSX2. POU4F2, POU6F2, KLF7, NPDC1, and TFAP2D were active in RGCs. SOX4 and MEIS1 were active in ACs. ONECUT1 and ONECUT2 were active in developing HCs. NEUROD1, NRL, CRX, RAX2, and PRDM1 were primarily expressed in PCs. VSX1 and VSX2 were expressed at high levels in BCs. In the RPE, MITF, BHLHE40, OTX2, and SMAD3 were active.

16Using SCENIC, the target genes of each transcription factor could be identified (Fig 4C). RAX2 was active in PCs. During the development of the retina, the increased target genes of RAX2 were involved in visual perception and detection of light stimulus, which is consistent with the function of PCs. KLF7 was active in RGCs. A set of KLF7 target genes was constantly expressed and related to neuron projection development. MITF was uniquely active and expressed in RPE, which targeted amino acid transport and neural tube closure–related genes at early stages and cation transport genes at late stages. MITF also targeted genes that were related to pigmentation during development and regulation of cell motion.

17SOX2 was mainly expressed in RPCs and Müller glia cells. POU4F2 was an RGC-specific transcription factor. NEUROD1 was detected in various retinal cells but was mainly expressed in PCs (S5 Fig). We verified these representative transcription factors by immunofluorescence staining and determined their distribution and expression in structures of the human fetal NR (S5 Fig). At 9 and 13 W, RPCs were SOX2+ and NEUROD1−. RPCs were located in the outer layers of the NR (near the RPE) at 9 W and then were shifted away from near the RPE toward INL at 13 W. Few RPCs were detected at 23 W, whereas Müller glia cells were observed in the INL of the NR at 13 and 23 W. At 13 W, the NEUROD1+SOX2− cells adjacent to the RPE were likely PCs. We detected SOX2+NEUROD1+ cells at 23 W, and they might be newborn neurons. The NEUROD1+SOX2− cells located next to Müller glia cells might be BCs. A majority of the BRN3B+ (gene POU4F2) (S5 Fig) cells were located farthest away from the RPE, and these cells were probably RGCs. At 9 W, the BRN3B−NEUROD1+ TFAP2B+ cells distributed among the RGCs might be ACs.

18We obtained 585 RGCs from 5 to 23 W. We performed pseudotemporal ordering of the single cells using Monocle2 to reveal the dynamic changes in the transcriptome of the developing human RGCs. We applied the extreme gradient boosting (xgboost) R package v0.7.0 training a model to identify any different RGCs based on the expression levels of all genes. The genes with the top 10 highest importance scores in this model were selected to obtain a pseudo–developmental time using the discriminative dimensionality reduction tree (DDRTree) algorithm of the Monocle package combined with the trajectory information (S7 Table). Based on our data, the pseudotemporal order of the cells matched their actual developmental stages (Fig 5A). Using the genSmoothCurves function, we analyzed the dynamic expression levels of RGC-enriched genes.

19During the RGC development, the expression levels of the bHLH genes ATOH7 and NEUROD1, the cell adhesion molecule CNTN2, and progenitor markers SFRP2 and ID1 decreased. We detected increased expression levels of RGC master genes, such as FXYD7, MEG3, NEFL, NEFM, SNCA, MAPT, NTM, and RXRG, in developing RGCs. The expression levels of ELAVL4, GAP43, PFN1, PRPH, RTN1, and TUBB2A remained constant in RGCs at different stages (S6 Fig).

20According to the results of the GO analysis, genes specific to early RGCs were involved in the response to BMP (P value = 2.0E-07). Genes specific to late RGCs participated in the energy derivation by oxidation of organic compounds (P value = 4.9E-11) and fatty acid beta-oxidation (P value = 3.1E-05) (S6 Fig).

21We stained fixed human fetal retinas for the RGC marker SNCG to obtain spatial information about RGC generation in the human retina. Immunostaining revealed the gradual generation of RGCs from the optic nerve head toward the periphery, and RGCs were essentially spread throughout the ganglion cell layer (GCL) at approximately 9 W (Fig 5D). The proportion and density of RGCs gradually decreased during the development of the human retina (Fig 5E). We also detected SNCG proteins in axons located in the optic nerve (Fig 5D1–5D3).

22We dissected the dynamic changes in the transcriptome of RPE cells using Monocle2. RPE cells were ordered in pseudotime consistent with the actual developmental stages (Fig 6A). The states of RPE cells could be divided into early, middle, and late states. We applied the R package (cummeRbund) to group RPE high-variable genes (S7 Fig, S8 Table). Nine clusters of RPE high-variable genes were identified. The down-regulated genes were related to cell division (Midkine [MDK], HSPA1A, and HSPA1B), cell cycle, regulation of neuron differentiation (ID3), and neural precursor cell proliferation (DCT, PAX6, SOX11, and WNT2B). The up-regulated genes were related to extracellular structure organization (CST3, EFEMP1, ITGAV, CRISPLD1, and ITGB8), lipid biosynthetic process (HACD3, PLA2G16, PLCE1, PTGDS, ABHD2, CYP27A1, and INPP5K), and the canonical retinoid cycle in rods (LRAT, PLTP, RGR, PLBP1, RPE65, and TTR) (Fig 6C and 6D). The expression of RPE65 was validated by immunofluorescence staining (Fig 6E and 6F). In addition, Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway mapping of RPE-enriched genes showed that amino acid signaling pathway–related genes were highly expressed at early stages (Fig 6G–6I).

23To identify the interactions between RPE cells and PCs in fetal human retina, we explored the expression levels of visual cycle–and phototransduction-related genes. First, we analyzed the states of PCs using PCA, and these cells were divided into three subgroups (S7 Fig). Groups 2 and 3, comprising cells from retinas at 23 and 24 W, were more mature PCs, as they expressed not only PC-enriched genes, such as NRL, AIPL, RP1, GNGT2, NEUROD1, PDC, RXRG, CRX, and PDE6G, but also genes related to phototransduction, such as RHO, SAG, GNGT1, and CNGA1 (S7 Fig). Group 1 seemed to be less mature PCs because of their lower expression levels of visual perception genes, and most of them were collected from earlier stages. SAG, GNGT1, and CNGA1 are involved in the rhodopsin-mediated signaling pathway. Moreover, we detected ABCA4, RBP3, and RHO in PCs (Fig 7B). At 24 W, more PCs expressed RHO, a rhodopsin expressed in rod cells that is also involved in the rhodopsin-signaling pathway. Light transduction was reported to begin with photon absorption by rhodopsin. Based on our results, light transduction might occur at 24 W. In RPE cells, RPE65, RDH5, LRAT, RGR, and RLBP1 were expressed at 13 W. In the visual cycle, RGR (also known as RPE-retinal G-protein-coupled receptor) is used to isomerize all-trans-retinal to 11-cis-retinal, and then the 11-cis-retinal is delivered back to PCs. RPE65, RDH5, and LRAT enzymes catalyze the reaction from all-trans-retinol to 11-cis-retinal. RLBP1 (also known as the cellular retinaldehyde binding protein [CRALBP]) speeds up the reaction (Fig 7A). Taken together, the visual cycle was probably initiated as early as 13 W, long before visual perception began.

24In addition, using CellPhoneDB (www.CellPhoneDB.org), we found more clues about interactions between RPE cells and PCs (Fig 7C and 7D). For example, we detected two neurotrophic factors MDK and PTN (pleiotrophin) in RPE cells. PTN and MDK are involved in inducing and stimulating neuronal differentiation. Receptor protein tyrosine phosphatase type Z (PTPRZ1), a shared receptor of MDK and PTN, was highly expressed in PCs. The binding of MDK to PTPRZ1 is important for MDK-dependent survival of embryonic neurons. The interaction of PTN with PTPRZ1 is known to play an important role in cell–cell adhesion, cell motility and migration, cell division, and epithelial–mesenchymal transition. The interactions of PTPRZ1 with both MDK and PTN might be key for the development of PCs and RPE cells.

25RetNet (https://sph.uth.edu/retnet/) has provided genes whose mutations cause inherited retinal diseases. We wondered whether human fetal retinal cells could be disease targets. In the human fetal retina, most inherited retinal disease–related genes were enriched in PCs, BCs, and RPE cells, such as autosomal dominant/recessive retinitis pigmentosa, X-linked macular degeneration, and other autosomal dominant/recessive retinopathy (Fig 8A and S8 Fig). A previous study showed that USH2A, EYS, and CRB1 were the top three genes contributing to inherited retinal dystrophy in Chinese patients. In our data set, USH2A was mainly enriched in PCs. EYS was highly expressed in BCs and PCs. Interestingly, CRB1 was enriched not only in Müller glia cells but also in RPCs. In addition, RPCs were also targets of CLRN1, KCNJ13, and KIF11 (Fig 8B). As RPCs are important for retinogenesis, mutations in disease-related genes in RPCs could be key for the study of disease processes.

26Microglia played a role in the pathogenesis of human age-related macular degeneration (AMD). In mice, microglia were reported to express genes associated with AMD. In addition to TLR4, which was also expressed in mouse microglia, we found that C2 and C3 were enriched in human microglia (Fig 8A and 8B).

27To our knowledge, the present study is the first to highlight the transcriptome landscape of the human fetal NR and RPE in parallel at a single-cell resolution. In general, the gene expression networks of the NR and the RPE showed substantial differences. NR-specific genes were related to nervous system development, whereas RPE-specific genes were related to retinol metabolism. However, both tissues underwent processes ranging from active cell proliferation to the termination of proliferation and the stepwise maturation of visual perception throughout development; these results indicated their functional interactions at late stages. The NR and RPE have been reported to codifferentiate together. Once the connection is lost, the single-layered RPE forms a multilayered retina-like structure.

28In the present study, all 2,421 cells were clustered into well-defined cell classes based on known markers identified through mouse models and adult human retina studies, indicating the conserved feature of these markers among mice and humans. Then, we revealed markers of human retinal cells and dissected the temporal order of the generation of retinal cells in vivo. In humans, RGCs were already present at 5 W and peaked at 8 W. HCs were detected at 7 W and peaked at 9 W. Subsequently, ACs peaked at 17 W, followed by PCs. BCs and Müller glia cells were observed last.

29We compared the markers of fetal NR cells with the markers of adult NR cells (S9 Fig). We found that most top 15 markers (10 out of 15) in adult human BCs were also detected in fetal human BCs, such as TRPM1, LRTM1, TMEM215, PLXDC1, NETO1, CA10, ST18, CHN2, SLC4A10, and KCNMA1. The top 30 RGC markers in adult retina showed less specificity in fetal RGCs. However, some other RGC markers were shared in adult and fetal retina, such as GAP43, PRPH, SYT4, INA, and NEFL. This finding might due to the fact that BCs detected in our study were from late embryos, whereas RGCs were from early embryos. Cell states in late embryos might be more similar to mature cells than those in early embryos. Our findings highlighted the differences between retinal cells at early stages and those at late embryonic stages or in adult retina. In addition, fetal PCs in our study highly expressed adult rod genes. The adult cone genes were seldom detected.

30We compared our single-cell RNA-seq data to the bulk transcriptomic studies of Hoshino and colleagues and revealed more details on transcriptome dynamics of fetal human retina (S9 Fig). In general, the DEGs between retinas at early and late time points found by bulk RNA-seq analysis could mapped on specific cell classes in our single-cell data. However, we found that some DEGs up-regulated at later time points might be expressed in the same cell class at early time points. For example, PROX1, an important transcription factor for generating HC, was reported to begin expressing after fetal day (D) 67 (9.5 W), whereas in our study, PROX1 could be detected in HCs as early as 8 W. Although our study could dissect more transcriptome features at single-cell resolution, it has limitations. As mentioned above, fetal PCs in our study highly expressed rod genes, and seldom were cone genes detected. In the bulk RNA-seq data, the cone gene GNAT2 was up-regulated starting at D80 (11.4 W). Cone genes might be detected if we expanded the cell numbers. We evaluated the specificity of each human fetal retinal cell class–enriched gene in retinas of macaque and mouse (S10 Fig). We found that more human fetal retinal cell class–enriched genes demonstrated enrichment in mouse retina than did those in macaque retina.

31One important issue in retinal development is the determination of retinal cell fate by transcription factors. We provided a key profile of transcription factors that were active and cell class–specific in developing human retinal and RPE cells. The functions of transcription factors in the developing human retina are difficult to verify; thus, an evaluation of the activities of these transcription factors is helpful in inferring their importance. We used a computational method called SCENIC to analyze the activity of the GRNs in each cell. TFAP2D has been detected in the developing mouse retina from E13.5 to E16.5. In the present study, TFAP2D and other transcription factors important for RGCs were detected in the developing human retina as early as 5 W. In human retina, ONECUT1 and ONECUT2 were active in developing HCs. In mice, Onecut1 and Onecut2 are expressed in both the neuroblast layer (NBL) and GCL at early stages (E11.5 and E12.5). Because of the high quality of our data set, the target genes of active transcription factors could be well defined. The functional annotations and dynamic expression patterns of active transcription factor target genes revealed major events during the development of retinal cells. Moreover, we verified transcription factors by immunofluorescence staining to explore the spatiotemporal gene expression network in the developing human retina. In conclusion, we identified various retinal cells in the human fetal NR. By exploring the expression patterns of transcription factors at both the transcriptome and protein levels, we provided important information about the spatiotemporal gene expression network in the developing human retina.

32We labeled RGC axons with SNCG (which might be involved in modulating axonal architecture during development) at different stages to trace the states of RGCs during retinal development. RGCs were concentrated in the inner layer of the NR, and the density decreased during NR development. RGCs were globally generated in the direction from the optic nerves head toward the periphery and almost spread throughout the inner layer at approximately 9 W. However, in the present study, the sample size was too small to do a serious cell type analysis.

33In the present study, RPE cells underwent sequential development with the decreased expression of PAX6 and increased expression of visual cycle genes, such as RPE65 and LRAT. We detected proliferating RPE cells at 5 and 6 W. In the present study, rhodopsin expression was first detected at 24 W, whereas genes related to visual cycle were already expressed in RPE cells as early as 13 W. In addition, the expression levels of genes related to photoreception, such as RHO, SAG, GNGT1, and CNGA1, were increased as PC differentiation proceeded. Moreover, using CellPhoneDB, we identified receptor–ligand pairs of RPE cells and PCs, such as PTPRZ1 and PTN/MDK, NGFR and TTR, and FLT1/NRP1/NRP2 and VEGFA, which also provided clues for RPE and PC interactions.

34The ability to sense light requires interactions between PCs and RPE cells. The present study revealed the development of the NR and RPE, which will improve our understanding of neural development and provide clues about neural regeneration. The RPE is a single-layered epithelium and participates in the visual cycle, outer-segment phagocytosis process, and retinal-blood barrier. The RPE is also associated with retinal degenerative diseases. Thus, studies of the development of both the human NR and RPE are important. We dissected the transcriptome landscape of the developing human NR and RPE, as well as the inherited retinal disease map on human fetal retina. This study provides key information about human fetal retinal development, which is likely to be important for understanding and curing pathological conditions.

35The present study was approved by the Reproductive Study Ethics Committee of Peking University Third Hospital (Research License 2012SZ-013 and 2017SZ-043). All donors signed the written informed consent and voluntarily donated aborted fetuses for the present study.

36We collected fetal eyes from voluntarily donated aborted fetuses at 5 to 24 W of gestation. The NR and RPE were separately dissected under microscopes. Then, the NR or RPE tissues were digested into single cells using 2 mg/mL Collagenase/Dispase (Roche) by incubating them at 37 °C for approximately 5 min.

37The modified STRT protocol for multiplexed single-cell RNA-seq was described in previous studies. Briefly, after digestion, single retinal cells were randomly selected and transferred into 2 μL of cell lysis buffer using a mouth pipette. The cell lysis buffer contained 0.1 U/μL RNase Inhibitor (TaKaRa, 2313B), 0.0475% Triton X-100 (Sigma-Aldrich, X100), 2.5 μM dNTP mixture (Thermo, R0193), and 2.5 μM barcode-reverse transcriptase (RT) primer (TCAGACGTGTGCTCTTCCGATCT-XXXXXXXX-NNNNNNNN-T25, with X representing the nucleotide of cell-specific barcodes and N representing the UMIs) (S9 Table). We thoroughly vortexed the tubes containing a single cell and the lysis buffer and incubated the tubes at 72 °C for 3 min to release the linearized RNA molecules and immediately placed the tubes on ice. The RT reaction was performed by adding 2.85 μL of RT mixture (40 U SuperScript II RT (Invitrogen, 18064071), 5 U of RNase Inhibitor, 5X Superscript II first-strand buffer, 25 mM DTT, 5 M betaine (Sigma-Aldrich, B0300), 30 mM MgCl2 (Sigma-Aldrich, 63020), and 1.75 μM TSO primer (AAGCAGTGGTATCAACGCAGAGTACATrGrG+G, rG represents riboguanosines one [+G] and +G indicates the LNA-modified guanosine) to each tube, and then we placed the tubes in a thermo cycler and incubated them at 25 °C for 5 min, 42 °C for 60 min, 50 °C for 30 min, and 70 °C for 10 min.

38After the RT reaction, we added 7.5 μL of the PCR mixture (6.25 μL of 2× KAPA HiFi HotStart Ready Mix [KK2602] and 300 nM of ISPCR oligo [AAGCAGTGGTATCAACGCAGAGT], and 1 μM of 3’Anchor oligo [GTGACTGGAGTTCAGACGTGTGCTCTTCCGATC]) to each tube. The amplification was performed using the following program: 4 cycles at 98 °C for 20 s, 65 °C for 30 s, and 72 °C for 5 min, followed by 15 cycles at 98 °C for 20 s, 67 °C for 15 s, and 72 °C for 5 min, with a final cycle at 72 °C for 5 min.

39Subsequently, we pooled the amplified samples with different cell barcodes together and purified these samples with DNA Clean & Concentrator-5 (Vistech, DC2005) and twice with 0.8X Ampure XP beads (Beckman, A63882).

40Purified cDNAs were amplified with biotinylated index primer (Biotin/CAAGCAGAAGACGGCATACGAGATindexGTGACTGGAGTTCAGACGTGTGCTCTTCCGATC) and ISPCR oligo for 4–5 cycles. Biotinylated cDNAs were first purified and subsequently sheared into fragments of approximately 300 bp using Covaris. Next, we enriched the fragmented cDNAs, which were primed by biotinylation and contained barcodes and UMI sequences, using Dynabeads MyOne Streptavidin C1 (Invitrogen, 65002).

41We constructed a library based on the enriched cDNA fragments, which were attached to the C1 beads, using KAPA Hyper Prep Kits (KK8505). We used the NEB U-shape adaptor for ligation. The final amplification was performed with the Illumina QP2 primer (CAAGCAGAAGACGGCATACGA) and the short universal primer (AATGATACGGCGACCACCGAGATCTACACTCTTTCCCTACACGAC) for seven cycles. Libraries were sequenced to generate 150-bp paired-end reads on an Illumina Hiseq4000 platform (sequenced by Novogene).

42A custom script was used to trim the primer, TSO, and polyA sequences for the raw input reads. Additionally, the cell barcode was also extracted and input into the read.id information of the fastq file. Subsequently, the extracted reads were mapped to the hg38 genome using Hisat2 v2.0.5. The unique mapped reads were then counted by FeatureCounts v1.5.2 using the Ensembl release 84 GTF transcriptome annotation. Finally, we obtained the following table (Table 1):

43Spark v1.6.2 was used to identify cells with at least 500 genes and then to cast the table into a sample (2,553)-gene (58,830) count matrix.

44We used the Seurat R/Bioconductor toolkit to perform the dimension reduction process. Gene expression levels were quantified as TPM, and the copy numbers of each transcript were quantified based on the number of distinct UMIs. Only cells expressing more than 1,000 genes and at least 10,000 transcripts were retained. Then, the variation for each gene was calculated, and only genes with a proper level of variation (fxn.x = expMean, fxn.y = logVarDivMean, x.low.cutoff = 0.5, x.high.cutoff = 10, y.cutoff = 0.5, and do.contour = F) were considered in subsequent analyses.

45After two steps, 2,421 samples with 1,230 high-variable genes remained for the PCA. The top 30 principal components were used for the JackStraw analysis, and 18 principal components were ultimately used in the RunTSNE analysis. Finally, the cell group information was manually annotated for each cell based on the t-SNE clustering group information.

46We used the Seurat function find_all_markers (test.use = "roc") to identify cluster-specific markers.

47We used the R package Seurat to identify subclasses in one cluster.

48First, we overlapped the DEGs found by Seurat with transcription factors identified by SCENIC. Second, to ensure transcription factors active, we retained transcription factors that were active among more than 50% of cells in a particular cell class.

49A detailed description of the SCENIC method is available at https://github.com/aertslab/SCENIC.

50Monocle v2.3.2 was used to construct the trajectory landscape of RGC/RPE development. As for RGCs, we applied the extreme gradient boosting (xgboost) R package v0.7.0 to modify the cell trajectory. First, RGCs were labeled with the group number identified by t-SNE in the previous steps to select representative genes. Then, xgboost R package v0.7.0 was used to train a model that could identify these different t-SNE labels based on the expression level of all genes. After the training process, genes with the top 10 highest importance scores in this model were selected. Subsequently, the trajectory was built using the DDRTree algorithm of the Monocle package, and we obtained a pseudo–developmental time from the trajectory data. Finally, the expression levels of the selected genes were plotted using the genSmoothCurves function.

51We applied the R package (cummeRbund) to group target genes of cell class–specific transcription factors using variable genesCluster (k = 3) and group RPE high-variable genes using variable genesCluster (k = 9).

52We used a software called CellPhoneDB—which has a database of ligands, receptors, and their interaction—to systematically predict cell–cell interaction molecules with default settings (https://github.com/Teichlab/cellphonedb). For the top receptor–ligand pairs, we applied R package ggalluvial to visualize the interaction modules.

53The eyes were fixed with 1–10 mL of 4% paraformaldehyde overnight at 4 °C. Eyes were dehydrated with 1–5 ml of a 20% sucrose solution overnight at 4 °C and embedded in Tissue-Tek O.C.T. Compound (Sakura #4583).

54For immunofluorescence staining, 12-μm sections were dried at 55 °C for 30 min and washed three times (10 min each) with 0.1% Triton in PBS (PBT). Sections were then permeabilized with 0.5% Triton in PBS for 30 min at room temperature and incubated with blocking buffer (0.1% Triton and 10% donkey serum in PBS) for 90 min at room temperature. Primary antibodies were diluted in blocking buffer. Subsequently, sections were incubated with diluted primary antibodies overnight at 4 °C. After three washes (20 min each), sections were incubated with diluted secondary antibodies (1:400) at room temperature for 2 h. Finally, sections were washed three times (20 min each) and incubated with Prolong Gold Antifade Reagent with DAPI (Invitrogen #1846939). Images were acquired using a confocal laser scanning microscope (Leica TCS SP8, Leica Microsystems, Wetzlar, Germany).

55The primary antibodies used in the present study were rabbit anti-SNCG (Abcam #ab55424, 1:500 dilution), mouse anti-RPE65 (Abcam #ab78036, 1:100 dilution), mouse anti-NEUROD1 (Abcam, #ab60704, 1:100 dilution), mouse anti-SOX2 (Abcam #ab79351, 1:100 dilution), rabbit anti-SOX2 (Abcam #ab59776, 1:100 dilution), rabbit anti-KI67 (Abcam #ab15580, 1:100 dilution), rabbit anti-BRN3B/POU4F2 (Abcam #ab56026, 1:100 dilution), and rabbit anti-AP2 beta (Abcam #ab221094, 1:100 dilution).

56We used the DAVID Functional Annotation tool (version 6.8, https://david.ncifcrf.gov/) to do gene-enrichment analysis. The P value listed after GO term in this paper is the EASE Score, a modified Fisher exact P value.

57We applied RNAscope assay to validate RNA detected in our single-cell RNA-seq data. Briefly, fresh human fetal retinas were fixed for at least 24 h in 4% paraformaldehyde (with DEPC) and proceeded to dehydrate with graded concentration of sucrose. After embedding with O.C.T. compound, 10-μm retina frozen sections were obtained for the following mRNA detection. According to the instructions provided by commercial RNAscope Multiplex Fluorescent Reagent Kit (Advanced Cell Diagnostics), the target probes ONECUT2 and NTRK1, which were synthesized by Advanced Cell Diagnostics, were hybridized to the retina sections. After hybridization, mRNA signal was amplified and then labeled with TSA Plus Fluorescein. Images were acquired through Nikon A1RSi+ confocal microscope.

58Codes for building the single-cell RNA-seq results browser (http://49.4.93.68:30004/) were deposed at https://github.com/huboqiang/singlecellbrowser. The needed files are S5 Table of this paper and GSE107618_Merge.TPM.csv.gz, which has been uploaded in GEO: GSE107618.