Single-Cell RNA Sequencing Reveals Molecular Features of Heterogeneity in the Murine Retinal Pigment Epithelium

Transcriptomic analysis of the mammalian retinal pigment epithelium (RPE) aims to identify cellular networks that influence ocular development, maintenance, function, and disease. However, available evidence points to RPE cell heterogeneity within native tissue, which adds complexity to global transcriptomic analysis. Here, to assess cell heterogeneity, we performed single-cell RNA sequencing of RPE cells from two young adult male C57BL/6J mice. Following quality control to ensure robust transcript identification limited to cell singlets, we detected 13,858 transcripts among 2667 and 2846 RPE cells. Dimensional reduction by principal component analysis and uniform manifold approximation and projection revealed six distinct cell populations. All clusters expressed transcripts typical of RPE cells; the smallest (C1, containing 1–2% of total cells) exhibited the hallmarks of stem and/or progenitor (SP) cells. Placing C1–6 along a pseudotime axis suggested a relative decrease in melanogenesis and SP gene expression and a corresponding increase in visual cycle gene expression upon RPE maturation. K-means clustering of all detected transcripts identified additional expression patterns that may advance the understanding of RPE SP cell maintenance and the evolution of cellular metabolic networks during development. This work provides new insights into the transcriptome of the mouse RPE and a baseline for identifying experimentally induced transcriptional changes in future studies of this tissue.


Introduction
Cells of the retinal pigment epithelium (RPE), an epithelial monolayer located between the neurosensory retina and the choriocapillaris, perform activities that are critical to ocular development and visual function [1][2][3]. As part of the outer blood-retinal barrier, RPE cells control the flow of electrolytes, water, gases, nutrients, and waste products between the retina and circulation that is essential for retinal development and homeostasis [1][2][3][4]. RPE cells phagocytose the tips of photoreceptor outer segments and thereby contribute to the daily turnover of the phototransduction machinery that maintains visual function [5,6]. These cells contribute directly to vision by regulating the concentration of ions in the subretinal space, which influence light-dependent electrophysiological responses in photoreceptor cells [7]. They also participate in the visual cycle, in which all-trans retinaldehyde released in photoreceptor cells upon light stimulation is reisomerized to the 11-cis configuration required for detecting additional stimuli [8]. RPE cells are heavily pigmented with melanin, which absorbs light to improve visual contrast and scavenges reactive oxygen species to maintain tissue homeostasis [9][10][11]. Although many genes and gene products that contribute to these functions are known, further studies are needed for a molecular understanding of how the RPE contributes to vision, retinal homeostasis, posterior eye development, and ocular disease. Analysis of the native RPE transcriptome represents a primary approach to this end. 2

of 23
The RPE cell population in the mammalian eye is heterogeneous [12], which may add complexity to transcriptomic analysis. Heterogeneity of morphological features among native RPE cells, such as cell area, shape, melanin pigmentation, and the number of nuclei per cell, is well documented in human [13][14][15][16][17][18][19][20][21][22] and other mammalian species [13,23,24], including mice [25,26]. Depending on species, morphological differences among RPE cells are accentuated with age and vary topographically with respect to the ocular region (central-peripheral, dorsal-ventral, nasal-temporal) and proximity to ocular specializations, including the macula, area centralis, visual streak, and tapetum lucidum [23]. Additionally, morphological differences have been observed in adjacent RPE cells or small patches of cells independent of topographical location, resulting in cellular mosaicism [12]. RPE cells are also functionally heterogeneous. For example, although the adult RPE is largely post-mitotic, studies of human eyes have identified a small population of stem cells that proliferate and differentiate when subsequently cultured in vitro [27]. Similarly, rare cells containing mitotic figures have been identified in the adult albino rat RPE [13], and a small population of mitotically active cells has been reported in the peripheral RPE of the adult rat [28,29], which may be related to human RPE stem cells. Analysis of another functional readout, differential indocyanine green dye uptake, revealed cellular mosaicism in the human and mouse RPE [30,31]. Further evidence of heterogeneity has come from the histological analysis of cellular components in ocular sections or RPE-choroid-sclera flatmounts or from the biochemical analysis of dissected regions of the posterior eye [12]. Finally, focal RPE changes have been observed in individuals affected with inherited macular diseases, such as butterfly-shaped pigment (or pattern) dystrophy [32] and Best vitelliform macular dystrophy [33], and in animal models of these diseases [32,34]. The non-uniform distribution of pathological changes in these diseases raises the possibility of a heterogeneous RPE response to the genetic and/or environmental conditions that induce disease. Overall, these studies provide compelling evidence for the topographic and cellular heterogeneity of the RPE. However, the underlying mechanisms that lead to this heterogeneity are poorly understood and may benefit from transcriptomic approaches that provide information at the single-cell level.
Transcriptomic studies have been described using RPE preparations from human donors and mouse models (Supplementary Table S1). The analysis of expressed sequences tags identified novel genes associated with the RPE [35,36]. Microarray analysis, based on the hybridization of known genes, extended these initial insights [37][38][39][40][41][42]. RNA sequencing (RNAseq) was used subsequently to improve the quality and depth of sequencing and to provide a view of absolute transcript abundance [43][44][45]. Most recently, single-cell RNA sequencing (scRNA-seq) has been used to examine cellular heterogeneity in human RPE and to provide clues to the development of this tissue [46][47][48].
Here, we apply scRNA-seq to RPE cells isolated directly from the mouse posterior eyecup by enzymatic and mechanical disaggregation and by further selection for viability based on fluorescence-activated cell sorting. A transcriptome of about 2000 highly variable genes was documented in each of roughly 2700 cells, enabling the use of cluster analysis to identify distinct but related RPE cell populations that appear to be distributed along a maturation time course. Bioinformatic analysis identified major known RPE pathways, including those related to visual cycle and melanogenesis genes, as well as additional transport and metabolic pathways that appear to be coordinately regulated upon RPE maturation. The approach described may benefit future efforts to understand the molecular basis of the RPE function in vision, development, homeostasis, and ocular disease.

Preparation of Single RPE Cells
To assess the cellular and molecular heterogeneity of RPE cells in an unbiased manner, we performed scRNA-seq on cells from two young adult male C57BL/6J (B6) mice at postnatal day 36 (P36). This age is three weeks past the last major wave of RPE cell division, which completes at about P15 [25]. To obtain single RPE cells, RPE-choroid-sclera eyecups were incubated with a concentrated trypsin solution and agitated gently to release RPE sheets, which were then disrupted mechanically to obtain single RPE cells ( Figure 1A). Cells from both eyes of each mouse were pooled to yield replicate samples (R1 and R2) and were incubated with calcein acetoxymethyl ester (calcein AM) to mark viable cells. In the experiment described here, fluorescence-activated cell sorting (FACS; Figure 1B, Supplementary Figure S1) of R1 and R2 yielded 27,657 and 12,540 single viable cells, respectively. These yields correspond to 26-28% and 12-13%, respectively, of the total population based on estimates of 5.4 × 10 4 and 4.9 × 10 4 RPE cells per adult B6 eye [25,26]. Single cells obtained by this approach were heavily pigmented and often exhibited two lobes of variable size ( Figure 1C). This shape is consistent with an apical and basal cellular domain separated by a junctional actin band as observed in other studies, in which single RPE cells were isolated from native or cultured sheets [49][50][51]. FACS-purified cells were concentrated and applied to the wells of an scRNA-seq chip to create a bar-coded cDNA library for sequencing ( Figure 1A). The elapsed time between enucleation and loading RPE cells onto the scRNA-seq chip was about 2 h; cells were kept on ice during this period except for the 30-min trypsin incubation at 37 • C and the 15-min calcein-AM incubation at room temperature.

Preparation of Single RPE Cells
To assess the cellular and molecular heterogeneity of RPE cells in an unbiased manner, we performed scRNA-seq on cells from two young adult male C57BL/6J (B6) mice at postnatal day 36 (P36). This age is three weeks past the last major wave of RPE cell division, which completes at about P15 [25]. To obtain single RPE cells, RPE-choroid-sclera eyecups were incubated with a concentrated trypsin solution and agitated gently to release RPE sheets, which were then disrupted mechanically to obtain single RPE cells (Figure 1A). Cells from both eyes of each mouse were pooled to yield replicate samples (R1 and R2) and were incubated with calcein acetoxymethyl ester (calcein AM) to mark viable cells. In the experiment described here, fluorescence-activated cell sorting (FACS; Figure  1B, Supplementary Figure S1) of R1 and R2 yielded 27,657 and 12,540 single viable cells, respectively. These yields correspond to 26-28% and 12-13%, respectively, of the total population based on estimates of 5.4 × 10 4 and 4.9 × 10 4 RPE cells per adult B6 eye [25,26]. Single cells obtained by this approach were heavily pigmented and often exhibited two lobes of variable size ( Figure 1C). This shape is consistent with an apical and basal cellular domain separated by a junctional actin band as observed in other studies, in which single RPE cells were isolated from native or cultured sheets [49][50][51]. FACS-purified cells were concentrated and applied to the wells of an scRNA-seq chip to create a bar-coded cDNA library for sequencing ( Figure 1A). The elapsed time between enucleation and loading RPE cells onto the scRNA-seq chip was about 2 h; cells were kept on ice during this period except for the 30-min trypsin incubation at 37 °C and the 15-min calcein-AM incubation at room temperature.

Characterization of RPE scRNA-seq Datasets
A total of 2667 and 2846 cells from R1 and R2, respectively, that passed all quality checks (see Methods) were analyzed. Unsupervised clustering of individual cell transcriptomes using Louvain community detection revealed six transcriptionally distinct clusters in both samples (Figure 2A,B, Table 1). By default, Seurat software labeled these as clusters 0-5 based on the population size of each cluster ( Figure 2C). We renumbered the clusters

Characterization of RPE scRNA-seq Datasets
A total of 2667 and 2846 cells from R1 and R2, respectively, that passed all quality checks (see Methods) were analyzed. Unsupervised clustering of individual cell transcriptomes using Louvain community detection revealed six transcriptionally distinct clusters in both samples (Figure 2A,B, Table 1). By default, Seurat software labeled these as clusters 0-5 based on the population size of each cluster ( Figure 2C). We renumbered the clusters produced by Seurat as C1-C6, based in part on the clustering tree map produced from the clustering trees tool [52], which indicates how clusters split and cells partition between clusters as resolution increases. The clustering tree map for single cell data from R1 revealed that cells in Seurat cluster 5 at a final chosen resolution of 0.6 split from the remaining cells in the first step (clustering resolution = 0.1), suggesting a distinct gene expression profile in these cells compared to other cells. Therefore, we renumbered Seurat cluster 5 as cluster C1 ( Figure 2C). On the other hand, Seurat clusters 1-3 appeared at later steps of this analysis (clustering resolution = 0.3, 0.6), so we renumbered these as clusters C4-C6 ( Figure 2C). Finally, Seurat clusters 0 and 4 were labelled as C2 and C3, respectively. The rationale for the final adjustment of the cluster order is given in Section 2.4. produced by Seurat as C1-C6, based in part on the clustering tree map produced from the clustering trees tool [52], which indicates how clusters split and cells partition between clusters as resolution increases. The clustering tree map for single cell data from R1 revealed that cells in Seurat cluster 5 at a final chosen resolution of 0.6 split from the remaining cells in the first step (clustering resolution = 0.1), suggesting a distinct gene expression profile in these cells compared to other cells. Therefore, we renumbered Seurat cluster 5 as cluster C1 ( Figure 2C). On the other hand, Seurat clusters 1-3 appeared at later steps of this analysis (clustering resolution = 0.3, 0.6), so we renumbered these as clusters C4-C6 ( Figure 2C). Finally, Seurat clusters 0 and 4 were labelled as C2 and C3, respectively. The rationale for the final adjustment of the cluster order is given in Section 2.4.   We performed correlation analysis to assess whether the clusters identified in each replicate were distinct and to test whether the replicates were similar. In each replicate, correlation analysis between clusters using average gene expression as a metric indicated strong positive correlation between Clusters C2-6 but weak correlation of these clusters with C1 (Supplementary Figure S2A,C). However, correlation analysis using the logarithm (base 2) of the fold-change in gene expression (log 2 FC) as a relative expression metric revealed negative correlations among the clusters (Supplementary Figure S2B,D), suggesting that cluster analysis identified distinct but related RPE cell populations. To compare clustering in R1 and R2, we measured the correlation between the average expression of genes in each cluster from both replicates ( Figure 2D) as well as the correlation between log 2 FC in each cluster from both replicates ( Figure 2E). We noticed substantial similarity between clusters in R1 and R2. C1 from R1 and R2 correlated strongly with each other, and clusters C2-6 correlated with each other ( Figure 2D,E). Taken together, these results suggest that the gene expression profiles in R1 and R2 are generally similar and confirm that the RPE cells isolated by our approach represent a robust heterogeneous population. Results from R1 are presented below, and those for R2 are summarized at the end of the Results section.

Functional Analysis of Clusters
To establish the molecular differences among clusters, we identified marker genes unique to each cluster by comparing the gene expression in each cluster against all other clusters. We found ≥25 cluster-specific marker genes for each cluster in R1 (adjusted p value [padj] < 0.05). A heat map comparing the fold-change in expression for the top 20 marker genes that distinguished each cluster from the others is shown in Figure 3A, and several marker genes are indicated for each cluster (the full list of marker genes is provided in Data S1). Overall, C1 showed the greatest differences in gene expression compared to C2-C6, which were relatively less distinguished from each other ( Figure 3A). The top 20 differentially upregulated genes in C1 included Mlana, Dct, Trpm1, and Gpnmb (Data S1), which participate in melanogenesis [53]. Upregulation of these genes indicates that C1 cells are positioned along the melanocytic developmental pathway. C1 also showed a higher expression of genes implicated in stem/progenitor maintenance and/or stemness, such as Aldoc [54], Dkk3 [55,56], and Id3 [54,57,58], compared to C2-C6. By contrast, the top 20 upregulated genes in clusters C2-C6 included RPE-specific marker genes, such as the visual cycle genes Rpe65, Lrat, and Rrh [8] ( Figure 3A; Data S1). These results suggest that C2-C6 consists of heterogeneous but closely related differentiated RPE cell populations.
We considered two possible origins for C1 cells. Melanocytes of the posterior eye include RPE cells, which originate in the neuroectoderm [59], and additional pigmented cells of the choroid, ciliary body, and iris, which ultimately derive from the neural crest [60]. The ciliary margin, ciliary body, and iris were removed from posterior eyecups in our studies by dissecting below the limbus; thus, C1 cells are either RPE cells or choroidal melanocytes. The top 20 genes upregulated genes in C1 include Pax6 ( Figure 3A; Data S1), a key transcription factor that orchestrates developmental processes in the eye [61], such as RPE development [62] and RPE melanogenesis [63]. By contrast, it is likely that choroidal melanocytes rely on Pax3 rather than Pax6 for development and/or melanogenesis, similar to neural crest-derived melanocytes in the skin and hair follicles [64][65][66]. This premise is supported by the labeling of choroidal but not RPE cells using a Pax3-green fluorescent protein reporter in mice [67], by the identification of Pax3 but not Pax6 in the transcriptional signature of murine choroidal melanocytes [68,69], and by reports that PAX3 mutations in patients with Waardenburg syndrome cause choroidal hypopigmentation without affecting RPE pigmentation [70][71][72]. Pax3 was not identified among C1 marker genes (Data S1). Together, these observations support the hypothesis that C1 comprises melanocytic RPE cells. Rpe65, Lrat, and Rrh [8] ( Figure 3A; Data S1). These results suggest that C2-C6 consists of heterogeneous but closely related differentiated RPE cell populations. We considered two possible origins for C1 cells. Melanocytes of the posterior eye include RPE cells, which originate in the neuroectoderm [59], and additional pigmented cells of the choroid, ciliary body, and iris, which ultimately derive from the neural crest In support of the above evidence for upregulation of five melanogenesis genes, GO analysis confirmed the enrichment for melanin biosynthetic processes in C1 ( Figure 3B; Data S2). GO analysis also indicated substantial differences in C1 energy metabolism compared to C2-C6, including an upregulation of genes associated with both oxidative phosphorylation ( Figure 3B; Data S2) and glycolysis (Data S2). Gene ontology (GO) analysis of differentially upregulated genes (padj < 0.05) in C2-C6 were enriched for multiple biological processes as indicated by GO terms related to WNT signaling, cellular response to metal ions, detoxification, and others ( Figure 3B; Data S2).
Melanogenesis genes in the RPE are abundantly expressed during early embryogenesis, prior to melanin accumulation, and are downregulated at later developmental stages [46,73]. Thus, the increased expression of melanogenesis genes in C1 may indicate that this population consists of melanocytes at an earlier developmental stage than those in C2-C6, possibly corresponding to RPE stem/progenitor (SP) cells. GO analysis did not identify SP processes among the differentially upregulated genes in C1 (Data S2). Nevertheless, many of the top 20 upregulated genes in C1 ( Figure 3A; Data S1) or their human orthologs appear to have important roles in SP cells. Six of these genes considered SP cell markers were upregulated, including Aldoc in quiescent neural stem cells (qNSCs) from the adult mouse hippocampus [74,75] and ventricular-subventricular zone of the brain [76], and in embryonic radial precursors of neural stem cells from the cortex [77]; Id3 in hippocampal qNSCs and cortex radial precursors [75,77], in quiescent limbal epithelial stem cells of the cornea [78], and (as its human ortholog ID3) in human embryonic stem cells [57]; Dkk3 in stem cells of the interfollicular epidermis [55] and DKK3 in human bone marrowderived mesenchymal stem cells [79]; Ifitm3 in stem cells and committed progenitors of the interfollicular epidermis [80], and in quiescent limbal epithelial stem cells [54,57,58,78]; MAP1B in human mesenchymal stem cells [81]; and Folr1 in ventral mesencephalic floor plate neural progenitors [82]. Four of the top 20 genes are upregulated in cancer stem cells, including Id3 in cancer stem cells of a mouse mammary tumor model [58] and ID3 in intrahepatic cholangiocarcinoma tumors [83]; MGST1 in human pancreatic cancer stem cells [84]; and GPX3 and GSN in quiescent colon cancer stem cells [85]. Importantly, several of the top 20 genes influence stemness, defined as the ability of SP cells to maintain an undifferentiated state capable of self-renewal and to differentiate into multiple cell types. Three of the top 20 genes promote the self-renewal of SP cells (DKK3 [79], Id3 [58,83,86], and Ifitm3 [78]) and with one other gene have been shown to influence the differentiation of SP cells into other cell types (DKK3 [79], Id3 [57], Ifitm3 [78], and Tmsbx4 [87]). In summary, there is strong literature support for the possibility that C1 cells are RPE SP cells expressing early melanocyte markers.

Assessment of Possible Retinal or Choroidal Cell Contamination
To identify additional possible sources of cellular heterogeneity present in R1 and R2, we first compared our cell clusters with mouse retinal cell clusters from a previous study [88]. C1 exhibited a significant positive correlation (p < 0.05) with most of the mouse retinal cell clusters, while clusters C2 and C4 did not correlate significantly with mouse retinal cell clusters ( Figure 4A). Clusters C5 and C6 exhibited significant negative correlation (p < 0.05) with mouse retinal cell clusters ( Figure 4A). The significant positive correlation of cluster C1 with mouse retinal cell clusters from 2-week-old mice [88] suggests that cells in C1 either are not RPE cells, despite the expression of Pax6 as discussed above, or possibly represent a multipotent RPE cell type that expresses retina-associated genes.
retinal cell clusters, while clusters C2 and C4 did not correlate significantly with mouse retinal cell clusters ( Figure 4A). Clusters C5 and C6 exhibited significant negative correlation (p < 0.05) with mouse retinal cell clusters ( Figure 4A). The significant positive correlation of cluster C1 with mouse retinal cell clusters from 2-week-old mice [88] suggests that cells in C1 either are not RPE cells, despite the expression of Pax6 as discussed above, or possibly represent a multipotent RPE cell type that expresses retina-associated genes.  To distinguish between these possibilities, we analyzed the expression of genes associated with specific RPE pathways, which we identified from the full list of differentially expressed genes (Data S1). Human RPE stem cells have been associated with low levels of visual cycle genes, high levels of melanin pigment biosynthetic genes and SP/stem cell markers [27]. Similarly, low expression of visual cycle genes and high levels of melanogenesis and SP cell marker genes were observed in C1 compared to C2-6 ( Figure 4B). These results provide further evidence for heterogeneity of the RPE cell population and for the assignment of C1 as an RPE SP cell type.
As an additional test of whether C1 and C2-6 represent bona fide RPE cells, rather than possible contaminating cell types, we investigated the expression of other genes in the cell clusters. Microglia are sometimes observed at the interface between the retina and RPE of B6 mice [89], and therefore are a possible contaminating cell type. We identified very low levels of key microglial cell type marker genes in all clusters ( Figure 4B). We also implemented the CELL-ID method [90] to verify the identity of the cell clusters using marker genes of microglial cell type as well as marker genes reported for the microglial cluster in the mouse retina [88] as a reference. CELL-ID annotated only five cells as microglia in R1 ( Table 2), indicating that there was little contamination from these cells in this dataset (<0.2% of total, <3% of any cluster). In R2, 26 microglia were identified, but remain a small percentage of the cell population (<1% of total, <3% of any cluster). Next, we examined the possible presence of choroidal melanocytes cells, as melanogenesis related genes, such as Pmel and Mlana, which were highly expressed in C1 ( Figure 4B; Supplementary Figure S3C), are also expressed in both mouse and human choroidal melanocytes [47,69]. We compared the top 100 gene signatures of each cell in C1-C6 with those in choroidal melanocytes from a mouse study [69] using the CELL-ID approach [90]. This approach identified eight cells in C2 with a significant choroidal melanocyte gene signature in R1, and one cell in C2 of R2 (Table 2; <0.4% of total, <1% of any cluster). None of the cells in cluster C1 were identified with this signature. Taken together, these results indicate that most of the cells in all clusters, including C1, correspond to RPE cells, and that contamination with microglia or choroidal melanocytes is rare.

A Proposed Cluster Maturation Timeline
Based on the expression of genes related to RPE-specific pathways such as visual cycle and melanogenesis and correlation with mouse retinal cell clusters ( Figure 4A,B), we propose C1 consists of immature RPE cells and C2-6 contain mature RPE cell types. We staged a possible maturation timeline from C1 to C6 based on well-known attributes of RPE cells (Figure 4C,D). Melanogenesis-associated genes exhibited significantly increased expression (log 2 FC > 1; padj < 0.05) in C1 relative to other clusters (Table 3, Figure 4C) and their expression declined progressively with maturation from C2-6 ( Figure 4C). Similarly, the expression of visual cycle and retinoid uptake genes was significantly reduced (log 2 FC < −1; padj < 0.05) in cluster C1 relative to other clusters (Table 4, Figure 4D). We also examined genes known to contribute to SP cell maintenance and renewal or are differentially upregulated in SP cells, including Aldoc, Dkk3, Id3, Tmsb4x [91], Anxa2 [92], Nbl1 [93], Rax [94], and Rarres2 [75]. Transcripts from these genes were more abundant in cluster C1 ( Figure 4B, Table 5) and declined in C2-6. Overall, these results support the identification of C1 as a SP cell population and the proposed maturation of RPE cells from C1 to C6.

k-Means Clustering and Functional Profiling
To identify other biological pathways that exhibit similar trends across the maturation timeline as the melanogenesis, visual cycle, and retinoid uptake pathways, we performed k-means clustering on differentially expressed genes (padj < 0.05) across all RPE clusters. Differentially expressed genes were classified into 10 different groups (Gp1-Gp10) based on their gene expression profiles along the proposed maturation timeline (Supplementary Data S3). Of these, gene sets in Gp5 and to a lesser extent Gp2 exhibited an almost identical expression profile as in our initial analysis of visual cycle and retinoid uptake genes along the maturation timeline (that is, downregulated in progenitor C1). In contrast, genes in Gp4 and Gp10 exhibited similar expression profile as melanogenesis genes (that is, upregulated in progenitor C1) along the maturation timeline ( Figure 5A). We then performed GO analysis on these gene sets to identify significant enrichment (FDR < 0.05) of multiple biological processes. Genes in Gp2 were enriched for biological processes such as: "regulation of lipid localization", "transport", "tissue migration" and "transforming growth factor" ( Figure 5B, Supplementary Data S4). Genes in Gp5 were enriched for biological processes such as "morphogenesis of an epithelial fold", "lipid localization", "homeostasis" and "ERK1 and ERK2 cascade" (Figure 5B, Supplementary Data S4). Gene sets in Gp4 were enriched for "metabolic process", "pigmentation, "establishment of cell polarity", and "neuron differentiation". Gp10 genes were enriched for biological processes, such as "oxidative phosphorylation", "negative regulation of immune system process", and "epithelial cell proliferation" (Figure 5B, Supplementary Data S4). Overall, we identified multiple biological processes that were differentially regulated in proposed progenitor/stem cluster C1.

Replication of Results in R2
To assess the reproducibility of these findings, we examined results from R2. As described above, unsupervised clustering revealed six transcriptionally distinct clusters in R2, similar to the clustering results in R1 ( Figure 2B). Average gene expression was highly correlated in respective clusters from R1 and R2 ( Figure 2C). Heatmaps of the top 20 marker genes distinguished C1 from C2-6 (Supplementary Figure S3A) and indicated that clusters C2-6 were relatively less distinguished from each other, suggesting heterogeneous but related RPE cell populations in R2. C1 showed higher expression of SP marker genes such as Dkk3, Id3, and Aldoc. Correlation analysis with mouse retinal cell clusters [88] identified significant positive correlation (p < 0.05) between C1 to most of the mouse retinal cell clusters (Supplementary Figure S3B), while C2 and C4 showed significant negative correlation (p < 0.05) with some of the mouse retinal cell clusters (Supplementary Figure S3B).
Moreover, in R2, C1 expressed high levels of melanogenesis and SP cell marker genes and low levels of visual cycle genes, while other clusters exhibited higher levels of visual cycle genes and reduced levels of progenitor markers (Supplementary Figure S3C). A maturation timeline along C1-6 was noted for R2, primarily based on expression profile of melanogenesis genes (Supplementary Figure S3D). However, the expression profile of visual cycle genes was not in complete agreement with this proposed maturation timeline (Supplementary Figure S3E), possibly indicating a lower sample quality in R2. Nonetheless, an overall maturation timeline representing C1 as an immature/progenitor RPE cell type and C2-6 as mature RPE cells was obtained. Importantly, C1 in R2 was highly correlated with C1 in R1 ( Figure 2D,E). Thus, we were able to identify SP cells in R2 as in R1.
Finally, we also performed k-means clustering on differentially expressed genes in each cluster in R2 followed by GO analysis of group of genes with similar expression profile as melanogenesis and visual cycle genes (Supplementary Figure S4). Groups of genes with reduced expression in C1 but increased levels in mature RPE cells (Gp1 and Gp2) were enriched for "regulation of lipid localization and transport", "response to toxic substances", and "transforming growth factor beta production" biological processes (Supplementary Figure S4B, Supplementary Data S4). Group of genes with increased expression in C1 but reduced in mature RPE cells (Gp5 and Gp8) were significantly enriched for "metabolic process" and "regulation of transport activity" and "establishment of cell polarity" (Supplementary Figure S4B, Supplementary Data S4). Overall, we observed similar cell cluster profiles in R2 and R1, reinforcing the presence of a SP population and confirming evidence for RPE cell heterogeneity due to differences in RPE maturation in young adult mice. Int. J. Mol. Sci. 2022, 23, x FOR PEER REVIEW 12 of 23

A Mature RPE Transcriptome
As indicated above, cluster analysis and k-means testing support the identification of C2-6 as a mature RPE cell population. Importantly, this population is free of contaminating cell types, and its transcriptome may therefore provide an opportunity to characterize the RPE transcriptome with high confidence. To allow such analysis, we determined the average transcript count for each gene in the combined C2-6 "mature" population that passed quality control criteria (Supplementary Data S5). Counts from C1 and from individual clusters C2-6 were provided in parallel. A useful (though arbitrary) inclusion threshold is an average expression value of one transcript per cell. As expected, signature RPE genes were expressed in the mature population at high levels above this threshold (Rgr, 186.9; Rpe65, 47.8), whereas photoreceptor genes were expressed below the threshold (Gnat1, 0.008; Rho, 0.02). Surprisingly, several genes associated with photoreceptors were detected at levels above the threshold (Rom1, 4.2; Slc24a1, 3.7; Abca4, 1.9; Gnb1, 1.9), raising the possibility that they are also expressed in the RPE (see Discussion). These data can be queried to assess the presence and abundance of specific gene transcripts in the mature RPE cell population.

Discussion
Transcriptional profiling of the mammalian RPE has been pursued to provide insights into RPE function in ocular health and disease [42,95]. In this report, we have demonstrated the use of scRNA-seq to assess the transcriptional profile of individual cells obtained from the native RPE of young adult mice. Our results may aid future efforts to explain morphological and functional RPE heterogeneity at the molecular level.

Cluster Analysis Identifies Two Major RPE Populations
Our results at P36 reveal a major population (>98% of total) of mature RPE cells with overlapping but distinct transcriptomic signatures (C2-6). The clusters are closely related but retain differences in genes associated with known pathways that reflect RPE development, such as the WNT signaling pathway [62], the visual cycle [8], and the cellular accumulation of copper and other metal ions, which contribute to melanogenesis and ultimately accumulate in melanosomes [96]. C2-6 may represent different stages in the RPE maturation process, which is not synchronized across the full posterior eye. Alternatively, the RPE at P36 may be fully mature, and the observed heterogeneity in gene expression may arise instead from differences in topographic location or from cellular mosaicism.
Our analysis also identified a small cluster of RPE cells (C1) with possible SP cell properties. Stem cells are characterized by self-renewal (the ability to proliferate in an undifferentiated state) and potency (the capacity to yield diverse differentiated states in response to suitable growth and differentiation stimuli) [97,98]. Progenitor cells are related to stem cells, but their self-renewal is limited to a small number of cell divisions, and their potency is limited to fewer cell types determined by commitment to a specific differentiation pathway [97,98]. Our results indicate a high correlation of gene expression between C1 and retinal cell clusters, possibly indicating a capacity to differentiate into multiple retinal cell types. High expression of selected stemness and SP cell maintenance genes was also observed among C1 cells. These are likely to be RPE cells, as they express melanogenesis genes at high levels, a known attribute of embryonic RPE [46,73]. Our analysis excluded another melanin-producing cell type, the choroidal melanocyte, as a major constituent of C1. These results suggest that C1 consists of multipotent RPE SP cells, which supports prior evidence for such cells in human and rodent eyes [13,[27][28][29].

Insights into RPE Maturation
GO analysis of clusters and k-means clustering reinforces the identification of C1 and C2-6 as SP and mature RPE cells, respectively, and provides new insights into possible RPE SP cell function. For interpreting these data, a relevant concept is the existence of a stemcell niche [99], which protects stem cells from injury due to the surrounding environment and immune system and which regulates the participation of stem cells in tissue growth, maintenance, and repair. Gp2 (R1) exhibits an overall trajectory of slightly increasing gene expression from C1 to C2-6. A major class of Gp2 GO terms involve development (kidney, glomerulus, cartilage, neural retina, nephron) or morphogenesis (ureteric bud, branching structure). Most of these annotations include the RPE genes Bmp4, Nog, and/or Sox9, which encode, respectively, the secreted growth factor and extracellular matrix (ECM) protein BMP4, its inhibitor NOG, which produces morphogen gradients based on its distribution relative to BMP proteins [100], and the developmental transcription factor SOX9, which regulates ECM production in diverse cell types [101][102][103]. Other ECM genes are variably associated with Gp2 GO terms, such as collagen genes Col4a3, Col8a1, and Col8a2. These results suggest that tissue growth and ECM production are important activities of the mature, differentiated RPE that are downregulated in the RPE SP niche. A second major class of Gp2 GO terms involves the transport of nutrients across the plasma membrane, dominated by solute carriers for amino acids, lipids, and energy metabolites. These results indicate an altered transport activity in C1 compared to the mature RPE, as might be expected as cells leave the stem cell niche and encounter a new environment.
Genes in the GO terms associated with Gp4 and Gp5 are predominantly those in the melanogenesis pathway and visual cycle, exhibiting decreased and increased expression, respectively, upon RPE maturation from C1 to C2-6. An additional gene of interest associated with diverse Gp4 GO terms describing cell polarity, proliferation, migration, and mitotic spindle orientation is Gja1. This gene encodes gap junction protein GJA1 (connexin 43), which participates in multiple pathways in the stem cell niche [104]. Its decrease in expression as cells mature from C1 to C2-6 is consistent with the departure of RPE cells from the stem cell niche. As in Gp2, several GO terms in Gp5 are associated with solute transport across membranes, consistent with an altered transport activity in the RPE SP cell population.
Genes associated with Gp10 GO terms indicate a shift in energy metabolism between C1 and C2-6. Gp10 transcripts from genes encoding glycolytic/gluconeogenic enzymes (Tpi1, Pkm) as well as mitochondrial components of the tricarboxylic acid cycle (Sdhb, Idh2) and oxidative phosphorylation pathway (Cox4i1, Atp5a1, Atp5c1, Atp5b, Atp5o) were relatively more abundant in C1 than in C2-6. Interestingly, increased isocitrate dehydrogenase IDH2 promotes conversion of α-ketoglutarate to citrate by reductive carboxylation, identified as a major RPE metabolic pathway [105]. Reductive carboxylation influences redox homeostasis, which is central to stem cell self-renewal [106,107], and contributes to RPE fatty acid synthesis [105], which is considered essential for human pluripotent stem cell survival [108]. The altered expression of energy metabolism genes among C1 cells may reflect differences in the local nutrient and redox status of the RPE SP cell niche.
Other interesting genes associated with Gp10 GO annotations include Id1, Tpm1, Pdlim4, and Cd47, which are common to terms involving the assembly of the actin cytoskeleton. TPM1 (tropomyosin 1 α) regulates actomyosin contraction in muscle and non-muscle cells [109], PDLIM4 promotes the formation of contractile actin bundles (stress fibers) [110], and CD47 regulates actin reorganization during induced cell death [111]. ID1 is a transcription factor that mediates cell stemness through negative regulation of basic helix-loop-helix transcription factors, thereby promoting self-renewal [112,113]. This protein promotes stress fiber formation in prostate epithelial cells upon treatment with transforming growth factor β1 [114] and regulates cytoskeleton remodeling during endothelial cell tubulogenesis [115]. These studies, together with our results, raise the intriguing possibility that RPE stemness may arise in part from an increased abundance of factors that regulate the actin cytoskeleton.

Possible Photoreceptor Gene Expression in the RPE
The use of scRNA-seq provides greater confidence in assigning transcripts to cell types that have traditionally been difficult to isolate or characterize in pure form. In our study, several transcripts were detected in the mature RPE cell population from genes that are thought to be expressed mainly in photoreceptor cells, including Abca4, Gnb1, Rom1, and Slc24a1. These transcripts are unlikely to be due to cell contamination but may arise from the phagocytic uptake of outer segments. However, it is also conceivable that some or all these genes are functionally expressed in the RPE. A prominent example is Abca4, which encodes a membrane flippase that transports vitamin A retinal-lipid adducts in the photoreceptor outer segment and is clinically associated with a form of Stargardt macular dystrophy [116]. The Abca4 gene and corresponding protein have recently been shown to be expressed both in photoreceptors and in the RPE [117], which may yield new avenues for understanding ocular vitamin A metabolism and interpreting the effect of ABCA4 mutations in Stargardt disease. Similar investigative opportunities may await other genes expressed in the RPE that have been predominantly characterized as photoreceptor specific. Interestingly, ROM1 variants cause an ABCA4-like macular dystrophy, raising the possibility that ROM1 dysfunction in the RPE contributes to this phenotype [118].

Relevance to Human RPE Heterogeneity and RPE SP Cells
Our studies parallel efforts to understand human RPE (hRPE) molecular heterogeneity. Early bulk transcriptional studies as well as more recent scRNA-seq experiments have reported transcriptional differences between the macular and peripheral hRPE [46][47][48], which may be important for understanding cellular processes that differentiate macular function and/or are targeted in macular disease. Although initial scRNA-seq studies yielded relatively few hRPE cells (289 cells from three donors, 54-92 years of age [47]; 5-185 cells at each embryonic time point [46]), which may limit the ability to detect heterogeneity, a later study yielded more cells (9302 from three donors, 29-64 years of age [48]) and identified multiple hRPE subpopulations. Similar to the multiple clusters identified in our studies, peripheral hRPE was found to include eight clusters differing in GO processes, such as extracellular matrix organization and nutrient transport. Macular RPE yielded two clusters differing in cell adhesion, ER stress response, and other GO processes [48]. Functional similarities between the human and mouse clusters were detected; for example, cellular response to the zinc ion was common to cluster P1 in hRPE and C5 in our study, and visual perception was common to subcluster P2-1 and C4. In addition, differentially expressed genes in early and late fetal hRPE [46] were also differentially expressed in our analysis comparing C1 to C2-6, further supporting the argument that C1 retains features of early embryonic RPE. Representative upregulated genes in early hRPE and C1 include DCT/Dct, PAX6/Pax6, ID3/Id3, and MDK/Mdk, and downregulated genes include TTR/Ttr, RPE65/Rpe65, and LRAT/Lrat [46]. Of interest, two scRNA-seq studies identified ID3 as a prominently upregulated gene in the macular RPE [47,48]. These findings, together with our result that the mouse ortholog Id3 is upregulated in possible RPE SP cells, lead to the interesting hypothesis that such cells may be present in the human macular RPE.
In addition to existing evidence for topographical differences in hRPE cellular and molecular features (for example, comparing peripheral with macular tissue; see Table S1), a recent study revealed regional heterogeneity and cellular mosaicism based on surveying hRPE cell morphology over the entire posterior eye [119]. The multiple RPE subpopulations identified by scRNA-seq in both the human and mouse raises the possibility that transcriptomic differences underlie this heterogeneity. Detailed in situ transcript analysis will be needed to tease apart regional and cellular differences in RPE function, which may lead to a deeper understanding of RPE activities in health and disease.
Previous reports indicate approximately 10% of cells isolated from the adult hRPE can be activated in vitro to yield self-renewing and multipotent cells capable of differentiation into neural and mesenchymal progeny [27,120]. One interpretation of these studies is that the isolated cells correspond to adult RPE stem cells [27,120]. Alternatively, SP cell properties may emerge due to transdifferentiation of tissue-derived mature RPE cells cultured in supplemented media [121][122][123]. It has been noted that there is no direct evidence for the existence of stem cells in the adult hRPE [124]. Our evidence for the existence of a small (1-2%) population of RPE SP cells in the native mouse tissue supports the hypothesis that hRPE stem cells may also exist in the human tissue, although additional studies are needed to confirm this hypothesis.

Limitations and Future Studies
The limitations of this study include the lengthy RPE cell isolation procedure, which may allow changes in RNA levels that prevent accurate determination of in vivo RNA abundance. Further, because replicate samples used in FACS analysis are stored for different lengths of time prior to loading on the single-cell processing system, RNA quality may vary among samples. Future refinements that shorten and/or synchronize the isolation procedure to minimize RNA degradation may benefit the approach. In addition, the recovery of a relatively small portion of the total RPE cell population may contribute to bias in assessing RPE heterogeneity. Methods that improve RPE sheet disruption may improve the proportion of the total population collected and thereby obtain a more representative cell population. Future efforts should include immunohistochemical or in situ RNA hybridization studies to validate the assessment of topographic and cellular mosaicism. Finally, longitudinal analysis using scRNA-seq may provide a better understanding of the sequential changes in cellular function that occur as the RPE population matures.

Mice and Colony Management
C57BL/6J (B6) mice were produced from an animal colony bred at The Jackson Laboratory (JAX, stock #000664). Mice were housed in the Research Animal Facility at JAX under a 12-h light/12-h dark cycle and provided an NIH31 (6% fat) diet and acidified water ad libitum. All animal procedures were approved by the Institutional Animal Care and Use Committees of the JAX and adhered to the ARVO Statement for the Use of Animals in Ophthalmic and Vision Research.

RPE Cell Isolation and scRNA-seq
Mice at P36 were sacrificed by asphyxiation with carbon dioxide and cervical dislocation. Eyes were enucleated and placed immediately in a dish of 1x phosphate-buffered saline (PBS) on wet ice. Eyes were punctured below the limbus with a 20 g needle, gripped at the cornea with straight forceps, and cut circumferentially below the limbus with angled Vannas microdissection scissors. After removal of the cornea, iris, and lens and trimming of the optic nerve flush with the sclera, the retina was immediately peeled from the posterior eyecup. The peeled eyecup was then placed in 1.0 mL prewarmed 0.5% trypsin-EDTA (Thermo Fisher, Waltham, MA, USA, 15400-054) and incubated for 30 min at 37 • C under an atmosphere containing 5% CO 2 . Following incubation, the interior of the eyecup was grasped at the optic nerve with forceps, and the entire eyecup was drawn out of the solution and resubmerged repeatedly to release RPE sheets. When no additional sheets were released, the remaining tissue containing the choroid and sclera was discarded. The solutions containing RPE sheets from both eyes were pooled with 8.0 mL of ice-cold collection buffer in a gentleMACS C tube (Miltenyi Biotec, Bergisch Gladbach, North Rhine-Westphalia, Germany) and disrupted using protocol C of a gentleMACS Dissociator (Miltenyi Biotec). Samples were centrifuged at 300× g for 10 min at 4 • C and resuspended by triturating the pellet in 250 µL DMEM (Thermo Fisher Scientific, 11054-020) containing 2% fetal bovine serum, 2 mM EDTA. Calcein-AM (0.1 µM) was added from a 40 µM working stock, and samples were incubated for 15 min. FACS was performed using a 100 µm nozzle on a FACSAria II cell sorter (BD Biosciences, San Jose, CA, USA) to isolate single cells based on side-and forward-scattering that were further gated for viability (DAPI-negative and calcein-positive). Sorted cells were collected into DMEM containing 20% fetal bovine serum, 1x N1 medium supplement (MilliporeSigma, Burlington, MA, USA, N6530-5BL), 1x MEM non-essential amino acids (Thermo Fisher Scientific, 11140050), 2 mM GlutaMAX-I (Thermo Fisher Scientific, A128601), 0.25 mg/mL taurine (MilliporeSigma T0625), 20 ng/mL hydrocortisone (MilliporeSigma, H0396), and 13 ng/ml triiodothyronine (MilliporeSigma, T5516). Collected samples were centrifuged at 300× g for 10 min at 4 • C and resuspended in 0.04% bovine serum albumin. A volume of the suspension containing 8000-12,000 cells was loaded in a single channel of a Chromium Single Cell Instrument (10x Genomics, Pleasanton, CA, USA), and barcoded cDNA libraries were prepared using a Chromium Single Cell 3 Chip Kit v2 (10x Genomics) according to the manufacturer's protocol. Amplified cDNA from each channel was used to construct a sequencing library (Illumina, San Diego, CA, USA). cDNA and libraries were checked for quality on a 4200 TapeStation (Agilent Technologies, Santa Clara, CA, USA) and quantified by KAPA qPCR. Sequencing was performed on an NextSeq500 System (Illumina, San Diego, CA, USA) using 150-cycle sequencing to an average depth of 50,000 reads per cell. The scRNA-seq data reported in this publication have been deposited in the National Center for Biotechnology Information Gene Expression Omnibus (GEO) with the accession number GSE203138.

scRNA-seq Data Processing
scRNA-seq data were processed using Cell Ranger v1.2.0 (10x Genomics; RRID: SCR_017344). Sequencing libraries were demultiplexed to individual cell FASTQ files utilizing the cellranger mkfastq function. Each library was aligned to an indexed GRCm38 (mouse) RefSeq genome with default parameters, followed by barcode counting, and UMI counting using Cell Ranger. In R1, 2844 single cells were sequenced with 191,744 mean reads and a median of 2111 detected genes per cell. In R2, 2944 single cells were sequenced with 128,744 mean reads and a median of 1813 detected genes per cell.
Downstream analysis was performed on filtered feature counts generated by Cell Ranger. SoupX v1.5.2 [125] was implemented to estimate and remove ambient mRNA contamination. We identified potential single-cell doublets using DoubletFinder v2.0.3 [126], with an expectation of a 4% doublet rate assuming Poisson statistics, as per the developer's code on GitHub. Further low-quality single cells containing <500 expressed genes or >5% mitochondrial transcripts were excluded from the analysis. Additionally, genes expressed in fewer than three single cells were also removed. Following the removal of lowquality and doublet cells, single cells were normalized and clustered using Seurat V4.0.0 (RRID:SCR_016341) [127]. Single-cell gene expression counts were normalized following a global-scaling normalization method with a scale factor of 10,000 and were log transformed using the Seurat NormalizeData function.
We applied principal component analyses to reduce the dimensionality of the data using the top 2000 most variable genes in the dataset. The top 10 principal components selected using the JackStraw and Elbow plot method were used in the RunUMAP analysis. Resolution parameters 0.6 and 0.5 were selected using the clustering tree method [52] to identify clusters from the R1 and R2 single cell datasets, respectively. Seurat V4.0.0 was used to identify cluster-specific marker genes, and visualization was performed with dot and feature plots. The genes specifically expressed in each cluster were examined to identify the cell types.
Separately, we also used the R package CELL-ID [90] to identify the cell types in our datasets. CELL-ID extracts unbiased per-cell gene signatures in a single-cell RNA-seq dataset and matches cells from the same cell type across independent datasets or reference datasets. In this study, we matched the top 100 gene signatures from each cell in our single-cell RNA-seq dataset with gene signatures from melanocytes cell types [69] and microglial cell types [88] from previous studies.

Correlation Analysis
We computed the Pearson correlation between our mouse RPE cell clusters and with mouse retinal cell clusters [88] using the cor.test function in R. Correlations that were significant at p < 0.05 are exhibited in correlation plots.

Functional Enrichment Analysis
Functional enrichment analysis was performed using the R package clusterProfiler [128]. Gene ontology analysis was performed using enrichGO functions from the clusterProfiler R package. The function compareCluster from this package was used to compare the enriched functional categories of each gene cluster. The significance threshold for all enrichment analyses was set to 0.05 using Benjamini-Hochberg corrected p-values.

Gene Expression Clustering of Marker Genes
Gene expression clustering of differentially expressed genes from each RPE cell cluster was performed using the K-means clustering approach. We used the relative expression (log 2 FC) of differentially expressed genes in each cluster to the group gene with similar expression profiles across the RPE cell clusters.