Evolving Transcriptomic Signature of Human Pluripotent Stem Cell-Derived Retinal Pigment Epithelium Cells With Age

We differentiated the human embryonic stem cell line H9 into retinal pigment epithelium (RPE) cells to assess temporal changes in transcriptomic profiles of cells. We performed single cell RNA-Sequencing of a total of 16,576 cells, and analysed the resulting data to access the molecular changes of RPE cells across two culture time points (1 and 12 months). Our results indicate the stability of the RPE transcriptomic signature over time in culture, with results indicating maturing populations of RPE could be observed over time, with no evidence of an epithelial – mesenchymal transition. Assessment of gene ontology (GO) pathways reveals that as cell cultures age, RPE cells upregulate expression of genes involved in metal binding and antioxidant functions, which might reflect an increased ability to handle oxidative stress as cells mature. Comparison with the transcriptional profiles of native RPE identified a progression towards a maturing RPE profile. These results suggest that in vitro long-term culture of RPE cells allow the modelling of phenotypes observed in native mature tissue. Our work highlights the changing transcriptional landscape of hPSC-derived RPE as they age in culture, which provides a reference for native and patient- samples to be benchmarked against.


Introduction
The retinal pigment epithelium (RPE) is a polarized monolayer of post-mitotic, pigmented cells that is essential to the health and function of photoreceptors and functionality of the underlying vasculature. The RPE phagocytoses and recycles photoreceptor outer segments -a waste product of visual cycling -and protects the retina against photo-oxidation by effectively absorbing light.
It governs the exchange of fluid, nutrients, and waste products to and from the apical and basal surfaces. Laterally, tight junctions reinforce the cell network, and intercellular gap junctions couple neighboring cells, forming a structure of metabolically connected cells. This arrangement forms part of the outer blood retina barrier, a physical barrier that efficiently isolates the neural retina from systemic circulation in the underlying vascular choriocapillaris. Importantly, the RPE is key to the immune privilege of the eye, by its physical contribution to the blood retina barrier, and also its expression of molecules repressing immune cells from migrating within the retina (1).
In the human retina, aging is associated with vision decline and delayed dark adaptation, both of which are direct consequences of tissue stress and retinal damage (2). It is hypothesized that overtime, oxidative stress leads to the death of retinal neurons; a decrease in RPE numbers; an accumulation of the toxic waste lipofuscin within the RPE; and an accumulation of basal toxic deposits called drusen underneath the RPE (2). Together, these events contribute to a loss of homeostasis and low-grade inflammation within the retina (2). Although it is clear that the RPE is key to the health of the retina, the precise molecular mechanisms responsible for its aging are not well understood. We need to better understand the causal pathways behind aging -and we can only do so if we move past current methodological impasse, which highlights the fact that we cannot collect samples from living individuals, particularly from tissue such as the eye. 5 Human pluripotent stem cells (hPSCs), including human embryonic stem cells (hESCs) and induced pluripotent stem cells (iPSCs), have generated great expectations in the field of regenerative medicine due to their ability to propagate indefinitely in vitro and give rise to any cell type in the body, including cells that form the retina. A variety of protocols have been described to differentiate hPSCs into RPE cells (3)(4)(5)(6)(7)(8)(9)(10). RPE cells are generally assayed after a few weeks of differentiation, at which stage they demonstrate similarity to their human native counterparts, in terms of morphology/expression of key proteins, functions and expression profiles, however with a profile closer to a foetal stage than adult stage (4,11,12). Interestingly, the transcriptome profile of hPSC-derived RPE cells as they age in culture is unknown. To date, the majority of RNA-seq studies have been conducted on bulk samples, consisting of millions of individual cells -the result of which is that transcript quantification represents the average signal across the cell population being studied. Recent developments to isolate single cells and genetically barcode their expressed transcripts has enabled the transcriptomes of single cells to be sequenced in a high throughput manner. By sequencing a large number of cells, it is now possible to dissect the cellular composition of apparently homogeneous cell cultures, providing a powerful way to identify and dissect molecular pathways involved in cellular events of interest.
Here, we performed single cell RNA-sequencing (scRNA-seq) on hESC-derived RPE cells differentiated for 1 and 12 months to assess the impact of time in culture on the RPE transcriptome and whether genetic hallmarks of maturation can be observed overtime.

Results
Quality control 6 The H9 hESC line was differentiated to RPE following the protocol described in the methods section. RPE cells originating from the same cell culture of origin and the same original passaging of RPE cells were isolated after 1 and 12 months of differentiation, dissociated to single cells and processed to generate libraries for scRNA-Seq analysis, in order to generate a transcriptional map of the aging RPE cells (Figure 1A). The capture of the 1-month single cell library detected 12,873 cells at the mean read depth of 40,499 reads per cell, while the 12-month capture detected 4,850 cells at the mean read depth of 114,503 reads per cell (Supplementary Table 1). To account for the disparities between the sequencing depth of the samples, the 12-  Table 1). The remaining 16,576 cells were retained for further analysis.

Identification and characterisation of subpopulations
Cluster analysis was performed independently and identified 12 subpopulations in each sample (Supplementary Table 2). After the data were integrated using anchors identified using a method described by (16), MetaNeighbor was used to match common subpopulations across both samples (17), denoted as "Common". Clusters unique to 1-month and 12-month samples were respectively denoted as "Young" and "Aged". In total, 18 subpopulations were identified with six common subpopulations ("Common", 8,484 cells), six subpopulations exclusive to the 1-month dataset ("Young", 5,758 cells) and six subpopulations exclusive to the 12-month dataset ("Aged", 2,334 cells). Conserved cluster markers for each distinct cluster (Supplementary   Tables 3-5) and cell counts per cluster ( Figure 1B) were identified, and clusters visualised by Uniform Manifold Approximation and Projection (UMAP) plot ( Figure 1C). 3,070 cells could not be assigned to a subpopulation. We compared variations in the transcriptomic profiles between the 1-month-old and 12-month-old samples as to identify potential changes in phenotypes upon aging of RPE cells in vitro. Cell type specific markers were identified using the Find Conserved Markers function for each cluster and Gene Ontology (GO) analysis (13,14) was performed using a PANTHER overrepresentation test against the Homo sapien genome   Figure 1D). As these transcripts are associated with stages of RPE maturity, our data suggests that all subpopulations are of RPE lineage, potentially at various stages of differentiation and maturation. Differential gene expression analysis and pathway enrichment were performed to characterise the molecular signature of these subpopulations.

Assessment of the Common Subpopulation
Of the cells examined, more than half (8,484 cells) were clustered into "Common" subpopulations (1-6) that intersect the 1-month-old and 12-month-old cell cultures (Figures 2-5 indicating a large shared transcriptional profile between the two conditions. A range of RPE markers were observed as conserved between conditions (1-month-old and 12-month-old) (Figures 2, 3A). In particular, RPE markers associated with melanin biosynthesis (MITF  Table 3). Examples of genes characteristics of neural differentiation, extracellular matrix and RPE are illustrated in Figure 3B-D respectively. Some intrasubpopulation variations in the pattern of gene expressions were observed between cells identified from the 1-month-old or 12-month-old cultures (Figure 3B-D). Process Overrepresentation test (Fisher exact test, FDR <0.05) analysis found an overrepresentation of pathways involved in mitochondrial and ribosomal function; protein biogenesis, transport, assembly and function; ATP biosynthesis and metabolism; and response to oxidative stress. Expression of mitochondrially-encoded genes and ribosomal genes has been correlated with development and maturation (19), including of the retina (20,21). Hence, the presence of RPE markers, together with mitochondrially-encoded genes and ribosomal genes suggests this subpopulation comprises of a highly metabolically active maturing RPE phenotype. abundantly produced by RPE cells (22) and its secretion diminishes with aging (23). DCT is expressed in the developing retina (24), is important to the production of melanin and to the RPE homeostasis (25,26) and its downregulation is associated with mature native RPE (27). We thus suggest this subpopulation is a mature functional RPE population.  (38)).
Only the marker CRABP1 was highly conserved and upregulated, whilst other RPE markers were downregulated (such as SERPINF1, BEST1, RLBP1 or RPE65). This thus suggests a subpopulation of immature cells. Altogether our data suggests that the common population to both time points is heterogenous, with subpopulations representing different stages of RPE cell differentiation.

Assessment of the Young Subpopulation
Around a third of all cells (5,625 cells) clustered into the "Young" subpopulation (Supplementary Table 4). Common RPE markers were found to be conserved markers within a few "Young" subpopulations (Figures 2, 3A, Supplementary Table 4). (Heterogeneous Nuclear Ribonucleoprotein H1). All play fundamental roles for the homeostasis and functions of the RPE: ATP1B1 encodes an apical Na + /K + ATPase which expression reduces with age and in AMD (39); defects in CTNNB1 are linked to abnormalities in RPE development and pigmentation (40); GSTP1 is a survival factor for RPE cells, which levels increase as cells mature (41); CD63 is a late endosome/exosome marker known to be released by RPE cells (42) and HNRNPH1 levels are associated with improved survival of RPE cells in culture (43). Hence, those markers indicate cellular functions suggestive of functional RPE cells. It is interesting to note that known canonical RPE markers, such as SERPINF1, RLBP1, TTR, PMEL, CRYAB (Crystallin Alpha B) were less expressed in this subpopulation than in all other subpopulations.
This could indicate a stage of differentiation/ maturation of the RPE cells. Taken together, our results suggest that all subpopulations within the "Young" cohort are immature cells developing to RPE cells.

Assessment of the Aged Subpopulation
Less than 10% of all cells (2,334 cells) clustered into the "Aged" subpopulation (Supplementary Table 5). All six identified "Aged" subpopulations were subjected to the same analysis; however, a statistical overrepresentation test returned no positive results for "Aged" Subpopulations Zero, One and Five (most likely owing to the low number of conserved genes identified). Only a few common RPE markers associated with melanin biosynthesis (TTR, DCT), extracellular structure organization (CST3, CRISPLD1), secretion (SERPINF1, VEGFA), visual cycle (RLBP1) and lipid biosynthesis (INPP5K) were conserved in some of the "Aged" subpopulations (Figures 2, 3A, Supplementary Table 5).  Table 5). This data suggests that the RPE cells of these "Aged" subpopulations are increasing their handling of metals and antioxidant abilities, which could reflect a further maturation of the RPE cells.

RPE cells express many complement components in various retinal diseases, inflammation
and/or aging (53). The complement regulators C1s (Complement C1s), C1r (Complement C1r) and C1QBP which form the C1 complex were conserved markers in a few "Young" and "Common" subpopulations (all in "Common" Subpopulation Three; C1s and C1r in "Young" Subpopulation Five and "Common" Subpopulations One and Five; C1r and C1QBP in "Common" Subpopulation Two; C1s "Common" Subpopulation Four; C1r only "Common" Subpopulations Two and Six; C1QBP only in "Common" Subpopulations Two, Three and Six) (Supplementary Tables 3, 4). Other genes associated with the complement response (such as APOE (Apolipoprotein E) plays a role in lipid metabolism, drusen content and is also modified with complement activation in the RPE (54). We observed that APOE is a conserved marker of the "Young" subpopulations One and Five, as well as all "Common" Subpopulations (Supplementary Tables 3, 4). The levels of expression per cell is however different between subpopulation, with an upregulation in the "Young" Subpopulation Five, and a downregulation in the "Young" Subpopulation One as well as in the "Common" subpopulations One, Two, Three and Six. Interestingly, within the "Common" subpopulation Four, cells arising from the 1-monthold culture expressed less APOE than cells arising from 1-month culture of all other "Common" subpopulations whilst the opposite was observed with the 12-month-old cells within this subpopulation. The opposite pattern was observed in the "Common" Subpopulation Five (Supplementary Tables 3, 4).

Comparison to known signatures of native RPE cells.
We compared the hPSC-derived RPE signature to that of fetal native RPE cells, which was obtained by scRNA-Seq (24). As observed foetal native tissue, the hPSC-derived RPE cells expressed SERPINF1, TYR, MITF, RPE65, BEST1, and TTR, and more immature cells expressed SFRP2, MKI67 and DCT (Figures 1D-3 TYRP1, SLC6A15 is associated with mature native RPE (27). We thus compared the expression profiles of those genes in the RPE cultures overtime, in order to assess the maturity of the cultured cells (27) (Figure 5). We observed that the "Common" Subpopulation (both times of culture) presented high expression of RPE65, TTR, BEST1, RBP1, PTGDS, SERPINF1, DUSP4 and GEM, associated with a downregulation of DCT and TYRP1 (Supplementary Table 3). This thus suggests that this population correlates more closely with the adult native RPE.

Discussion
Here, we provide a dynamic profile of the transcriptome of hPSC-derived RPE cells over 12 months. Our data confirms expression of markers of RPE homeostasis and functions in hESCderived RPE (24,27,55,56) and provides novel information on the timing of expression of those markers. At both early and late time points, we observed that hPSC-derived RPE cells expressed genes associated with secretion, visual cycle, melanin synthesis, phagocytic activity, metal binding and oxidoreductase activity. Based on expression of genes known to be associated with levels of RPE maturity and on GO analysis, the identified 18 subpopulations regroup into three main populations: immature and progenitor cells, maturing RPE cells and functionally mature RPE cells. Within those, some subpopulations were comprised of highly metabolically active cells.
An essential function of the RPE is photoprotection of the retina, which is accomplished by different mechanisms. Those include absorbing radiation, binding and sequestering redox-active metals such as iron, and scavenging free radicals and reactive oxygen species (57).
Metallothioneins are metal-binding proteins that are protective against oxidative stress.
Compared to the 12-month-old RPE, the 1-month-old RPE cells express less mRNA for the metallothioneins MT1E, MT1F, MT1G, MT2A, MT1X and higher levels of the Dopachrome Tautomerase DCT. This data suggests a variation in the handling of metals and in the antioxidant abilities of RPE cells, which could be reflective of either a necessity to handle more oxidative stress in an aging in vitro environment or a maturation of RPE cells towards a more mature and protective phenotype (57). The assessment of variations in other gene expression between the two time points indicates a maturation profile of cells rather than an increased stress. Indeed, DCT is known to be expressed in human retinal progenitor cells (24) and its expression is regulated by the early RPE marker MITF. It is thus not surprising that as RPE cells mature, MITF expression reduces (58) and subsequently reduces DCT expression, as is observed in human fetal retina (24). Similarly, CTNNB1 regulates MITF and OTX2 expression and subsequently RPE differentiation (40). Finally, SOX11 is known to be expressed in early retinal progenitor and early in differentiating RPE cells (24,48), hence its downregulation as cell culture ages further supports a maturation of RPE cells in culture. The "Young" Subpopulations Zero and One are characterised by the expression of CST3, which encodes for the cysteine proteinase inhibitor Cystasin. Interestingly, this protein is known to decrease in native RPE cells with aging (23). Its presence in the cell population further strengthens the suggestion of a maturing RPE population over time.
The high level of expression of mitochondrial and ribosomal -related genes in some cell populations also indicates that cells are very active, necessitating energy and ribosomal activities for protein synthesis. It could also suggest that ribosomes potentially contribute to extraribosomal functions in the RPE cells, such as cell development and maturation (19), as already reported for melanocyte development (59), retinal development (20,21) and retinal degeneration (60). 22 Our analysis also revealed that cells in culture can develop a transcriptomic profile more closely related to the adult native RPE. In particular one subpopulation found in both 1-monthold and 12-month-old cultures showed a profile of common RPE gene expression more closely related to the adult than to the fetal cells (27). This population could potentially be selected for assessment of more mature characteristics and phenotypes.

Conclusion
The

Materials and methods
Ethical approval 23 The experimental work was approved by the Human Research Ethics Committee of the

University of Melbourne (1545484) with the requirements of the National Health & Medical
Research Council of Australia (NHMRC) and conformed with the Declaration of Helsinki.

Cell culture and differentiation of hESCs to RPE cells
The hESC line H9 (Wicell) was maintained on vitronectin-coated plates using StemFlex (Thermo Fisher), with medium changed every second day (61). Cells were differentiated into RPE cells as previously described (3) with the following modifications. Briefly, hESCs were maintained in culture until 70-80% confluent at which stage StemFlex was replaced with Complete E6 (Stem Cell Technologies) supplemented with N2 (Thermo Fisher) to induce retinal differentiation, with thrice weekly media changes for 33 days. On Day 33, medium was replaced with RPEM (α-MEM, N1 supplement, 5% fetal bovine serum, non-essential amino acids (all Thermo Fisher), penicillin-streptomycin-glutamine (Sigma), taurine-hydrocortisone-triiodothyronin (THT) (in-house preparation) to promote RPE differentiation, with media changed every second day. Cells were cultured for 32 days, after which point maximal pigmentation is routinely observed. Cells were harvested with an 8-minute exposure to 0.25% Trypsin-EDTA (Thermo Fisher) and inactivated with RPEM. Non-RPE contaminants (visible as unpigmented cells), were manually removed from the culture, which begin shedding off the culture plate after ~ 2 minutes. Cells were seeded at a density of 75,000 cells / cm 2 onto growth-factor-reduced Matrigel-coated tissue culture plates (Corning). Media was changed every second day, with the first sample of cells harvested after 30 days (D30) and the second sample harvested on day 367 (D367) for scRNA-Seq analysis. (Figure 1A). 24 RPE cells were dissociated to single cells using 0.25% Trypsin-EDTA (Thermo Fisher) for 8 min and inactivated with RPEM. Cells were centrifuged at 300g for 1 minute to pellet and resuspended in a small volume of RPEM containing 0.1% v/v propidium iodide (PI, Sigma Aldrich) to exclude non-viable cells. Single cell suspensions were passed through a 35 µm filter prior to sorting. A minimum of 60,000 live cells (PI-negative) were sorted on a BD FACSAria IIU (100 µm, 20psi) into culture medium. Cells were centrifuged at 300g for 5 min and resuspended in PBS containing 0.04% BSA to a concentration of ~800-1,000 cells/µl. Approximately 17,400 cells were loaded onto a 10X chip for a target recovery of 10,000 cells.

Bioinformatics mapping of reads to original genes and cells
Reads underwent initial quality control, sample demultiplexing and quantification with Cell Ranger 3.0.2 by 10x Genomics (http://10xgenomics.com). Raw base calls were demultiplexed into individual samples and converted to FASTQ using the cellranger mkfastq pipeline. The reads were then mapped to the Homo sapiens genome (GRCh38, Annotation: Gencode v29) using STAR (v2.5.1b) called via the cellranger count pipeline. Resulting data for each sample were then aggregated and depth-equalized via the cellranger aggr pipeline.

Quality control and normalisation
Each sample underwent quality control independently. First, the following values were calculated for each cell: total number of Unique Molecular Identifiers (UMIs), number of detected genes and proportion of mitochondrial and ribosomal-related genes to total expression. Cells that had any of the 4 parameter measurements lower than 3x median absolute deviation (MAD) of all cells were removed from subsequent analysis. Cells were also removed if mitochondrial-associated genes accounted for more than 25% of total expression, and/or ribosomal-associated genes accounted for more than 60% of total expression. Cell-cell normalization was then performed using the SCTransform function from Seurat (62).
Confounding sources of variation -specifically, mitochondrial and ribosomal mapping percentage, were also regressed out with this function.

Dimensionality reduction, clustering and integration
The dimensionality of the data was first reduced with Principal Component Analysis (PCA).
Subsequently, the 30 most statistically significant principal components (PCs) were reduced to two dimensions by Uniform Manifold Approximation Projection (UMAP). These values were used to construct a Shared Nearest Neighbour (SNN) graph for each cell. The Louvain method 26 for community detection was then used to identify clusters in each dataset at the resolution of 0.6. The 1-month and 12-month datasets were then integrated using Seurat's SCTransform integration method (16). The unsupervised version of MetaNeighbor was used to evaluate the similarities between the 1-month clusters and 12-month clusters. Cluster pairs that received an Area Under the Receiver Operating Characteristic (AUROC) score greater than 0.9 were merged into one cluster (Supplementary Figure 2).

Clustering characterization and analysis
Network analysis was then performed on significant differentially-expressed genes using Reactome functional interaction analysis (63,64). Gene Ontology (GO) analysis was performed using (13,14).

Code availability
All code and usage notes are available at: https://github.com/powellgenomicslab/RPE_scRNA_AgedStudy. This repository consists of code used to process raw sequencing data in FASTQ format to cell-gene expression tables via the Cell Ranger pipeline, and code used to perform the following analysis: quality control, normalization, dimensionality reduction, clustering, differential expression and integration.

Data Records
Data is available at ArrayExpress (accession number E-MTAB-851). Files consist of raw FASTQ files, and a tab separated matrix of UMIs per gene for each cell passing quality control filtering. BAM files can be generated by using the supplied repository to process the FASTQ files via Cell Ranger.   Figure 1C). The intensity of gene expression is indicated by colour variation.