Zinc Supplementation Induced Transcriptional Changes in Primary Human Retinal Pigment Epithelium: A Single-Cell RNA Sequencing Study to Understand Age-Related Macular Degeneration

Zinc supplementation has been shown to be beneficial to slow the progression of age-related macular degeneration (AMD). However, the molecular mechanism underpinning this benefit is not well understood. This study used single-cell RNA sequencing to identify transcriptomic changes induced by zinc supplementation. Human primary retinal pigment epithelial (RPE) cells could mature for up to 19 weeks. After 1 or 18 weeks in culture, we supplemented the culture medium with 125 µM added zinc for one week. RPE cells developed high transepithelial electrical resistance, extensive, but variable pigmentation, and deposited sub-RPE material similar to the hallmark lesions of AMD. Unsupervised cluster analysis of the combined transcriptome of the cells isolated after 2, 9, and 19 weeks in culture showed considerable heterogeneity. Clustering based on 234 pre-selected RPE-specific genes divided the cells into two distinct clusters, we defined as more and less differentiated cells. The proportion of more differentiated cells increased with time in culture, but appreciable numbers of cells remained less differentiated even at 19 weeks. Pseudotemporal ordering identified 537 genes that could be implicated in the dynamics of RPE cell differentiation (FDR < 0.05). Zinc treatment resulted in the differential expression of 281 of these genes (FDR < 0.05). These genes were associated with several biological pathways with modulation of ID1/ID3 transcriptional regulation. Overall, zinc had a multitude of effects on the RPE transcriptome, including several genes involved in pigmentation, complement regulation, mineralization, and cholesterol metabolism processes associated with AMD.


Introduction
The retinal pigment epithelium (RPE) is a highly polarized monolayer of cells lining the back of the eye which provides critical support for the functioning of the adjacent photoreceptors. It is part of the outer blood-retina barrier that regulates the transport of

Retinal Pigment Epithelial (RPE) Cell Culture
Primary human fetal RPE cells (ScienCell, Carlsbad, CA, USA) from one donor were purchased and used at passage three (P3) for the complete study in duplicates/triplicates with unknown clinical or genetic background. Cells were seeded onto Corning 6-well transwell inserts (10 µm thick polyester inserts with 0.4 µm pore size, 4 × 10 6 /cm 2 pore density, Corning, Wiesbaden, Germany) in 125.000/cm 2 of epithelial cell medium (EpiCM, ScienCell, Carlsbad, CA, USA). After one week in culture, cell culture media were replaced with Miller medium with 1% FBS [18,19] and cells were cultured for two, nine, and nineteen weeks in duplicates. Two types of short-term zinc treatment were also conducted, where one-one extra replicates of untreated controls were taken for the two types of zinc treatment experimental setup. After one week or eighteen weeks in culture, cell culture media were replaced with Miller medium with 1% FBS for an additional one week in the absence or presence of 125 µM externally added zinc (as zinc sulphate; Thermo Fisher Scientific, Waltham, MA, USA) both in the apical and basal chambers, resulting in~10 nM bioavailable or free zinc [15,20]. The resulting replicates were the following: duplicates of zinctreated samples, triplicates of untreated controls at the two-and nineteen-week time point, and duplicates of untreated controls at the nine-week time point. Cellular differentiation was monitored through the development of cobblestone cell morphology and increase in pigmentation using light microscopy. The increase in transepithelial resistance (TEER) was measured using the EVOM2 Epithelial Voltohmmeter and STX2 electrodes (World Precision Instruments, Sarasota, FL, USA).
At the sample collection time, as detailed above, cells were washed with PBS (Thermo Fisher Scientific, Waltham, MA, USA) two times for one minute. Cells were detached by incubation with 0.15 % Trypsin-EDTA for thirty minutes at 37 • C. The trypsinization was stopped using 100% FBS and trypsin neutralization solution (ScienCell, Carlsbad, CA, USA). The obtained single-cell suspensions were washed in PBS with 1% BSA (Thermo Fisher Scientific, Waltham, MA, USA) 2 times for 5 min at 1000 rpm. After automatic cell counting (EVE, Thermo Fisher Scientific, Waltham, MA, USA), 7 × 105 cells/mL were prepared, and the cells were kept on ice for a maximum of ten minutes before proceeding with single-cell RNA sequencing. In parallel to single-cell sequencing, adjacent samples were fixed for fifteen minutes in 4% PFA (Merck, Darmstadt, Germany) diluted in PBS (Thermo Fisher Scientific, Waltham, MA, USA) for immunofluorescence.

Experiment Overview
Our previous study showed individual differences in assaying primary hfRPE from different donors [17]. To overcome the variations introduced by variability in donor samples and to generate a reproducible zinc effect, in this manuscript, experiments were performed on primary hfRPE cells from a single donor. In the initial scRNA-Seq run, samples were obtained from RPE cells cultured for two weeks (2W), nine weeks (9W), and nineteen weeks (19W) (Supplementary Figure S1A) in duplicates. Cells were collected from two wells at these time points. A total of 7000 cells from each sample were loaded on 10× Genomics Chromium v1.3 with a target recovery of 4000. Libraries made from each sample were pooled and sequenced.
In the second run, samples originated from RPE cultures were treated with a zincsupplemented medium for one week either after: (1) one week in culture or (2) eighteen weeks in culture in duplicates. We also included one-one sample from untreated RPE culture in this run and the transcriptomic profiles were generated in a pooled fashion as described above. The actual cell recovery of both runs ranged from 3000 to 4000 in each well, resulting in a total recovery of~30,000 cells for the first run and~15,000 for the second run. The raw scRNA-Seq data were processed using CellRanger v3.0.0. and then Seurat v3.1 to determine the heterogeneity of our specimens using unsupervised clustering, followed by annotation based on hierarchical clustering of a pre-defined set of canonical RPE marker genes [21][22][23][24] (Supplementary Table S3). For further analysis, we initially analyzed our samples of untreated control RPE cultures from the two runs (triplicates for 2W and 19W and duplicates for 9W cultures). We then separately analyzed the duplicate samples of our zinc-treated RPE cultures compared to the triplicate samples of untreated control RPE cultures of 2W and 19W.

scRNA-Seq
Approximately 7000 single cells per sample were processed with the Chromium system using the v3 single-cell reagent kit (10× Genomics, San Francisco, CA, USA). Barcoded libraries were pooled and sequenced on the NovaSeq platform (Illumina, San Diego, CA, USA), generating 150 bp paired-end reads as per 10× Genomics recommendations, with >30,000 reads per cell.

Bioinformatics
The raw scRNA-Seq data were processed using CellRanger version 3.0.0 (10× Genomics). The resulting filtered expression matrices were then imported into R for analyses using scRNA-Seq packages, Seurat (Version 3. Cells were filtered to exclude those with <1000 or >8000 genes, or with >20% of counts aligned to mitochondrial genes, or >40% counts aligned to ribosomal genes. Cells passing QC were downsampled randomly to 1000 cells per sample to prevent over-or under-representation of any sample. Each sample was log-normalized using default Seurat parameters, with the top 3000 highly variable genes used for Seurat iterative pairwise integration. The integrated dataset was scaled to regress variance arising from read depth and mitochondrial and ribosomal expression. Principal Component Analysis was then performed on the integrated dataset, and Seurat's JackStraw function was applied to determine the components used in UMAP and SNN clustering. Unsupervised clustering was run iteratively at resolutions ranging from 0.25 to 1, at increments of 0.25. At the highest resolution, a total of 13 clusters were detected. These clusters were observed in UMAP to form two overall, as-yet unannotated cell populations. Using untreated cells only, the average expression for the clusters was determined for a set of 213 canonical RPE marker genes [21][22][23][24] (Supplementary Table S3) to which hierarchical clustering was applied. The clusters were segregated into two distinct branches, exhibiting characteristics of more and less differentiated RPE, which matched the distinction observed in UMAP. As such, the 13 unsupervised clusters were annotated to reflect these two overall cellular populations for downstream differential expression analysis.
Seurat's Wilcoxon rank sum test was used for differential expression testing, using default FindMarkers parameters, with genes below 0.05 adjusted p-value considered significantly differentially expressed.
Monocle 3 [25] was used for pseudotime analysis, for which downsampled count data were imported from Seurat and independently processed and batch-corrected in Monocle using default parameters. For continuity, a pseudotime trajectory graph was calculated and projected on the UMAP coordinates preserved from Seurat analysis. The data were filtered to focus on the main less differentiated to more differentiated pseudotemporal trajectory, by excluding small branches not contributing to the main trajectory. This was followed by graph autocorrelation analysis to detect gene expression changes correlating with progress along the trajectory, filtered for significance at p-value and Q-value <0.05. Genes with expression significantly correlated with the trajectory were grouped into 'modules' of co-regulated genes and the average expression of each gene module calculated across pseudotime.

Functional Classification Pathway and Network Analysis
For pathway and network analysis, we used the GeneAnalytics (https://ga.genecards. org/#input; accessed on 14 March 2021) and STRING (Search Tool for the Retrieval of Interacting Genes/Proteins) 11.0 (https://string-db.org/; accessed on 14 March 2021) [26] in combination with Cytoscape [27]. GeneAnalytics uses binomial distribution to test the null hypothesis that the queried genes are not over-represented within any superpath, GO term, or compound in the GeneAnalytics data sources. The presented score in each section is a transformation of the resulting p-value, corrected for multiple comparisons using the false discovery rate (FDR) method, with higher scores indicating a better match. The bar color, indicating the matching quality-high (dark green), medium (light green), low (beige)-is common for all sections. STRING in combination with Cytoscape implements classification systems such as Gene Ontology, KEGG, and systems based on high-throughput text mining and the used reference dataset was the human genome. The identified functional protein association network was validated via text mining, database information, co-expression, and experimental evidence.

Immunofluorescence
For immunofluorescence analysis, the cells on the transwell membrane were permeabilized in 0.5% Triton-X (Merck, Darmstadt, Germany) in PBS for ten minutes at 4 • C and then washed in 0.1% Tween20 in PBS (PBST) (Merck, Darmstadt, Germany) and blocked with 5% goat sera (Merck, Darmstadt, Germany) in PBST for one hour at room temperature. Samples were then incubated with primary antibodies overnight: COL1A1 (Abcam plc, Cambridge, UK, dilution 1:200) and RPE65 (Merck Millipore, Darmstadt, Germany, 1:50), diluted in PBST containing 1% goat sera. Following washing with PBST, the samples were incubated with secondary antibodies in 1:200 in PBST with 1% goat sera for one hour in the dark at room temperature. Samples were washed with PBST for five minutes, then with PBS. Cell nuclei were then labelled with DAPI (Thermo Fisher Scientific, Waltham, MA, USA) diluted 1:1000 in PBS. Finally, samples were mounted onto Menzel-Glaser slides (Thermo Fisher Scientific, Waltham, MA, USA) in Vectashield (Vector Laboratories, Burlingame, CA, USA). For negative control, the primary antibody labelling was omitted. Cells were visualized using a Leica SP8 confocal microscope (Leica, Wetzlar, Germany).
Images were obtained and analyzed with Leica Application Suite X Image software (Leica, Wetzlar, Germany).

Maturation of RPE Cells in Culture
Primary human fetal RPE cells from a single donor were cultured for 2 weeks (2W = short term), 9 weeks (9W = medium term), and 19 weeks (19W = long term) (Supplementary Figure S1A). We used culture conditions that, in our hands, reproducibly recapitulated key aspects of RPE cells as described in previous studies [15][16][17]. As time in culture increased, RPE cells were observed to develop pigmentation, hexagonal morphology (Supplementary Figure  S1B), and a progressively increasing epithelial barrier function (112.9 ± 3.9 Ohm × cm 2 at 2W, 195.2 ± 16.6 Ohm × cm 2 at 9W, and 201.36 ±49 Ohm × cm 2 in 19W in culture). The cell cultures also began accumulating sub-RPE deposits (Supplementary Figure S1C) containing lipids and hydroxyapatite that we have shown earlier [15][16][17]. To identify the transcriptomic profiles of RPE cells at the three time points, we collected cells from three wells at 2W and 19W and two wells at 9W in culture (see Section 2 for detail). Approximately 3000-4000 cells were captured from each well and processed on the 10× Genomics Chromium v1.3 platform, with transcriptomes generated for a total of 30,000 cells.

Unsupervised Clustering Analysis
To ensure equal representation from all conditions, all samples were downsampled to include an equal number (1000) of randomly selected cells in Seurat 3.1 [28]. Based on 3417 differentially expressed transcripts (Supplementary Table S1), the cells were automatically allocated into thirteen clusters and visualized on a Uniform Manifold Approximation and Projection (UMAP) plot ( Figure 1A). The lists of cluster-specific 'marker' genes were input into GeneAnalytics. The gene set analysis tool identified significant cluster-specific canonical pathways [29], labelled as superpathways for the functional analysis of the cell populations. Supplementary Table  S2 contains information on the 'marker' genes, numbers of enriched pathways, and matched number of genes to the total number of genes in a pathway for each cluster. The software assigned the 312 genes in Cluster 0 to 18 superpathways, with respiratory electron transport and heat production of uncoupling proteins, metabolism, and visual cycle among the top five hits. The 205 genes in Cluster 1 were assigned to 44 superpathways, with degradation of extracellular matrix, ERK signalling, and phospholipase C pathway amongst the top five hits.
The 270 genes in Cluster 2 were associated with 19 superpathways with metabolism, The lists of cluster-specific 'marker' genes were input into GeneAnalytics. The gene set analysis tool identified significant cluster-specific canonical pathways [29], labelled as superpathways for the functional analysis of the cell populations. Supplementary Table S2 contains information on the 'marker' genes, numbers of enriched pathways, and matched number of genes to the total number of genes in a pathway for each cluster. The software assigned the 312 genes in Cluster 0 to 18 superpathways, with respiratory electron transport and heat production of uncoupling proteins, metabolism, and visual cycle among the top five hits. The 205 genes in Cluster 1 were assigned to 44 superpathways, with degradation of extracellular matrix, ERK signalling, and phospholipase C pathway amongst the top five hits.
The 270 genes in Cluster 2 were associated with 19 superpathways with metabolism, respiratory electron transport, heat production of uncoupling proteins, and visual cycle amongst the top five hits. In Cluster 3, the 210 differentially expressed genes were associated with 76 superpathways with cytoskeletal signalling, ERK signalling, and focal adhesion among the top five hits. The 87 differentially expressed genes in Cluster 4 were associated with 29 superpathways with cytoskeletal signalling, ERK signalling, and integrin signalling among the top five hits. The 30 differentially expressed genes in Cluster 5 were associated with two superpathways: melanin biosynthesis and tyrosine metabolism. The 270 differentially expressed genes in Cluster 6 were associated with 50 superpathways, degradation of extracellular matrix, metabolism of proteins, and cell adhesion and ECM remodelling amongst the top five hits. In Cluster 7, 313 differentially expressed genes were associated with 110 superpathways with cytoskeletal signalling, ERK signalling, and degradation of extracellular matrix among the top five hits. The 303 differentially expressed genes in Cluster 8 were associated with 61 superpathways, degradation of extracellular matrix, ERK signalling, and phospholipase C pathway among the top five hits. In Cluster 9, we identified 302 differentially expressed genes associated with 21 superpathways with metabolism, visual cycle, and copper homeostasis among the top five hits. In Cluster 10, 198 differentially expressed genes were associated with 24 superpathways with metabolism, visual cycle, and oxidative stress among the top five hits. The 807 differentially expressed genes in Cluster 11 were associated with 86 superpathways, degradation of extracellular matrix, protein processing in the endoplasmic reticulum, and cytoskeletal signalling amongst the top five hits. Finally, in Cluster 12, 164 differentially expressed genes were associated with 11 superpathways with organelle biogenesis and maintenance, intraflagellar transport, and mitotic cell cycle among the top five hits. The identification of thirteen clusters shows that cells in culture are not homogenous.

Hierarchical Clustering Analysis Using Markers of Mature RPE Cells
We aimed to identify which of the unsupervised clusters most resemble mature RPE. We separated cells that were deemed to be more differentiated based on the expression of 213 RPE-specific genes we identified from several publications [21][22][23][24] (Supplementary  Table S3). Hierarchical clustering based on the gene list divided the 13 clusters into two distinct groups ( Figure 1B). We annotated clusters 0, 2, 5, 9, 10, and 12 as 'more differentiated' RPE cells and the remaining clusters (1, 3, 4, 6, 7, 8, and 11) 'less differentiated' cells ( Figure 1B). We use the terms 'more differentiated' and 'less differentiated' from this point forward. The more and less differentiated cells are separated on the original UMAP ( Figure 1(C1)). We calculated the proportion of more and less differentiated cells at the 2W, 9W, and 19W. Interestingly, nearly half of the cells were more differentiated even as early as 2W or 9W in culture (2W = 41%, 9W = 41%). By 19W, the proportion of the more differentiated cells increased to 73% (Figure 1(C2)), with 27% remaining less differentiated. Supplementary Table S4 lists the genes that define the more and less differentiated groups. The expression levels of three highly expressed representative genes from each group are presented as violin plots and UMAP plots in Supplementary Figure S2, highlighting the enrichment but not the exclusive presence of these genes in one or the other group.
Next, we tested whether the protein products of the genes that distinguish more and less differentiated cells show differential expression. One of the highly expressed mRNAs in the less differentiated cells was Collagen Type I alpha 1 chain (COL1A1), fibril-forming collagen. The RPE secretes the protein encoded by this gene and it is found in the sub-RPE space [30]. We found that the expression of COL1A1 gradually increased in the less differentiated group and decreased in the more differentiated group (Figure 2A). In contrast, Retinoid Isomerohydrolase (RPE65), a visual cycle component marker for differentiated RPE, was mainly expressed in the more differentiated group (Figure 2A). Both genes were expressed in the other group but at a low level in the opposing groups ( Figure 2B). Next, we determined the immunolocalization of the COL1A1 and RPE65 proteins in the 19W RPE monolayer. In line with the gene expression results, the cells with a strong RPE65 immunolabelling also had weak intracellular immunoreactivity for COL1A1 proteins, and cells with strong immunolabelling for COL1A1 showed weak labelling for RPE65 ( Figure 2C; green: RPE65, red: COL1A1). Immunolabeling of COL1A1 is also present in the sub-RPE space. This extracellular immunoreactivity gradually increased with time in culture (Supplementary Figure S3), suggesting that the secreted COL1A1 accumulates as part of the developing extracellular sub-RPE material (Supplementary Figure S3).
As we identified more and less differentiated RPE cells in our hfRPE, we investigated whether more and less differentiated cells are also present in RPE cells directly isolated from human eyes. We used two independent previously published datasets: the scRNA-Seq data obtained from human embryos [31] or adult human eyes [32] (Supplementary Figure S4). We applied our cell grouping strategy based on 213 RPE-specific signature genes (Supplementary Table S2). Indeed, our analysis showed that both the embryonic RPE (Supplementary Figure S4(A1,A2)) and adult RPE (Supplementary Figure S4(B1,B2)) could be classified into more and less differentiated cell populations. Of note, the number of cells analyzed in the publication using adult RPE [32] was relatively low. Hence, clusters were less well separated. Next, we determined the immunolocalization of the COL1A1 and RPE65 proteins in the 19W RPE monolayer. In line with the gene expression results, the cells with a strong RPE65 immunolabelling also had weak intracellular immunoreactivity for COL1A1 proteins, and cells with strong immunolabelling for COL1A1 showed weak labelling for RPE65 ( Figure 2C; green: RPE65, red: COL1A1). Immunolabeling of COL1A1 is also present in the sub-RPE space. This extracellular immunoreactivity gradually increased with time in culture (Supplementary Figure S3), suggesting that the secreted COL1A1 accumulates as part of the developing extracellular sub-RPE material (Supplementary Figure S3).
As we identified more and less differentiated RPE cells in our hfRPE, we investigated whether more and less differentiated cells are also present in RPE cells directly isolated from human eyes. We used two independent previously published datasets: the scRNA-Seq data obtained from human embryos [31] or adult human eyes [32] (Supplementary Figure S4). We applied our cell grouping strategy based on 213 RPE-specific signature genes (Supplementary Table S2). Indeed, our analysis showed that both the embryonic RPE (Supplementary Figure S4(A1,A2)) and adult RPE (Supplementary Figure S4(B1,B2)) could be classified into more and less differentiated cell populations. Of note, the number of cells analyzed in the publication using adult RPE [32] was relatively low. Hence, clusters were less well separated.

Pseudotemporal Ordering of the Expressed RPE Genes
To identify the genes associated with transitioning from the less to the more differentiated cells, we performed a pseudotemporal ordering of our scRNA-Seq transcriptome profile using Monocle3 (Figure 3). This unsupervised analysis identified a main trajectory with 11 nodes (Figure 3(A1)). Based on the original cluster analysis depicted in Figure 1, node 1 corresponded to the less and node 10 to the more differentiated cells (Figure 3(A2)). The main trajectory was correlated with 537 variably expressed genes. Based on their pseudotemporal expression profile, these clustered into seven modules (Figure 3(B1); Supplementary Table S5). Modules 2 and 5 contained 175 genes with high expression at the early stages of the trajectory that gradually declined towards the end of the trajectory. Ge-neAnalytics identified 62 potential significant superpathways associated with these genes (Supplementary Table S6). Degradation of extracellular matrix, focal adhesion, and cell adhesion-endothelial cell contacts were amongst the top five ranked pathways. Modules 3, 4, and 6 contained 172 genes. These gradually increased towards the late stages of the trajectory. GeneAnalytics identified nine potential superpathways defined by these genes, including the transport of glucose, metabolism, and visual cycle among the top five hits (Supplementary Table S6). Modules 1 and 7 contained 190 genes. The expression of these genes transiently increased to a maximum at the middle of the trajectory followed by a decrease over pseudotime. These identified sixteen superpathways with degradation of extracellular matrix, ERK signalling, and cytoskeleton remodelling among the top five hits (Supplementary Table S6).

Pseudotemporal Ordering of the Expressed RPE Genes
To identify the genes associated with transitioning from the less to the more differentiated cells, we performed a pseudotemporal ordering of our scRNA-Seq transcriptome profile using Monocle3 (Figure 3). This unsupervised analysis identified a main trajectory with 11 nodes (Figure 3(A1)). Based on the original cluster analysis depicted in Figure 1, node 1 corresponded to the less and node 10 to the more differentiated cells (Figure 3(A2)). The main trajectory was correlated with 537 variably expressed genes. Based on their pseudotemporal expression profile, these clustered into seven modules (Figure 3(B1); Supplementary Table S5). Modules 2 and 5 contained 175 genes with high expression at the early stages of the trajectory that gradually declined towards the end of the trajectory. GeneAnalytics identified 62 potential significant superpathways associated with these genes (Supplementary Table S6). Degradation of extracellular matrix, focal adhesion, and cell adhesion-endothelial cell contacts were amongst the top five ranked pathways. Modules 3, 4, and 6 contained 172 genes. These gradually increased towards the late stages of the trajectory. GeneAnalytics identified nine potential superpathways defined by these genes, including the transport of glucose, metabolism, and visual cycle among the top five hits (Supplementary Table S6). Modules 1 and 7 contained 190 genes. The expression of these genes transiently increased to a maximum at the middle of the trajectory followed by a decrease over pseudotime. These identified sixteen superpathways with degradation of extracellular matrix, ERK signalling, and cytoskeleton remodelling among the top five hits (Supplementary Table S6). To identify potential transcriptional regulators of the pseudotemporal trajectory, GO term analysis was carried out in GeneAnalytics. This identified transcriptional regulator activity in two genes, ID1 and ID3, belonging to the combined Module 1 and 7, representing the transitional phase on the pseudotemporal trajectory. To identify potential transcriptional regulators of the pseudotemporal trajectory, GO term analysis was carried out in GeneAnalytics. This identified transcriptional regulator activity in two genes, ID1 and ID3, belonging to the combined Module 1 and 7, representing the transitional phase on the pseudotemporal trajectory.
Next, we examined the potential relationships between the most highly significantly correlated in the main trajectory using a cut-off value of Moran I = 0.5 (see Section 2). We identify 44 genes strongly influencing the main trajectory. Using GeneAnalytics in combination with STRING database and Cytoscape, we found that these genes do not appear to be randomly distributed. A total of 31 out of 44 genes showed a significant biological connection, validated via text mining, database information, co-expression, or experimental evidence (p < 1.0 × 10 −16 ; Figure 3(B2)). The 44 genes were associated with six potential superpathways, including visual cycle, extracellular matrix degradation, and cell adhesion-extracellular matrix remodelling as the top hits (Supplementary Table S6).

Transcriptional Changes in Response to Acute Zinc Supplementation
Previous studies have shown that chronic zinc supplementation has clinical benefit associated with molecular and cellular changes [13,15,17,33,34], but the effects of acute or short-term zinc supplementation had not been studied in detail. To identify the effects of short-term zinc supplementation, we treated our RPE cultures for one week with a zinc-supplemented medium, using the same approach as we described earlier [15]. This acute zinc supplementation was carried out on less differentiated cells starting at the end of the first week in culture. Then, cells were harvested at the end of 2W or more differentiated cells at the end of the 18th week, and cells were harvested at the end of 19W. Gene expression changes with zinc supplementation were compared to cells in culture without zinc supplementation for either 2W or 19W. Cells with and without zinc supplementation were clustered using the process used in Figure 1(C1). While acute zinc supplementation did not noticeably change the proportion of the more and the less differentiated cells ( Figure 4A), it significantly changed the expression of 472 genes in the more differentiated cells (Figure 4(B2)) and 149 genes in the less differentiated cells (Figure 4(B1)) at the twoweek time point (Supplementary Table S7).
At 19W, zinc altered the gene expression of 487 genes in the more differentiated cells (Figure 4(B4)) and 417 genes in the less differentiated cells (Figure 4(B3)) (Supplementary Table S7) (logFC > 0.25, adjusted p-value < 0.05). We displayed the four datasets in a four-way Venn diagram to further analyze specific temporal zinc-induced gene expression changes ( Figure 4C). We found 81 overlapping genes differentially expressed under all four conditions. Two-thirds of these 81 genes were identified as housekeeping genes by GeneAnalytics, confirming previous studies showing that zinc plays a role in regulating cellular homeostatic processes [35]. Relevant proteins include metallothioneins (MT1E, MT1F, and MT1X) that act as essential stress proteins to regulate immune homeostasis. In the more differentiated cells, 222 uniquely affected genes were at 2W and 163 at 19W ( Figure 4C). In the less differentiated cells, only four genes were specifically affected by zinc supplementation at 2W and 94 genes at 19W (Figure 4C).
At 2W, we identified superpathways only in the more differentiated cells; these were cytoskeleton remodelling, focal adhesion, and degradation of extracellular matrix among the top five superpathways (Supplementary Table S8). At 19W, in the less differentiated cells, we identified presenilin signalling, SMAD signalling, and antigen-presenting crosspresentation amongst the top five superpathways (Supplementary Table S8). In contrast, in the more differentiated cells, we identified metabolism, ferroptosis, and protein processing in the endoplasmic reticulum amongst the top five superpathways (Supplementary Table S8). Information on the magnitude and direction of zinc-associated change in transcript abundance of these gene lists is provided in Supplementary Table S7. The analysis of these five gene lists by GeneAnalytics to identify superpathways is listed in Supplementary Table S8.  Table S7) (logFC > 0.25, adjusted p-value < 0.05). We displayed the four datasets in a fourway Venn diagram to further analyze specific temporal zinc-induced gene expression changes ( Figure 4C). We found 81 overlapping genes differentially expressed under all four conditions. Two-thirds of these 81 genes were identified as housekeeping genes by GeneAnalytics, confirming previous studies showing that zinc plays a role in regulating cellular homeostatic processes [35]. Relevant proteins include metallothioneins (MT1E, MT1F, and MT1X) that act as essential stress proteins to regulate immune homeostasis. In the more differentiated cells, 222 uniquely affected genes were at 2W and 163 at 19W (Figure 4C). In the less differentiated cells, only four genes were specifically affected by zinc supplementation at 2W and 94 genes at 19W ( Figure 4C). ; volcano plot visualization of the differentially expressed genes in less (B1) and more (B2) differentiated cells in two-week culture and in less (B3) and more (B4) differentiated cells in 19-week culture; number of differentially expressed genes overlapped amongst the different cell types and culture times (C); overlap between pseudotimecorrelated genes and differentially expressed genes following acute zinc supplementation (D1) and a network representation of the 16 overlapping genes using Gene Analytics in combination with String and Cytoscape (D2), in which grey lines represent validated connections via text mining, database information, and co-expression, pink lines represent experimentally validated network connection, and the thickness of lines indicates the strength of supporting data.

Influence of Zinc on Transcription Dynamics
We next determined the overlap between the 537 genes identified in the main trajectory in the pseudotemporal analysis (Figure 3 Table S7). This comparison identified 16 common genes (Supplementary  Table S9). Using GeneAnalytics in combination with STRING database and Cytoscape, we found that these 16 genes show significantly (p-value < 1.0 × 10 −16 ) more interactions than expected, validated by text mining, database information, co-expression, and experimental evidence (Figure 4(D2)) that relates to the respiratory electron transport and response to metal ions as biological function (Supplementary Table S9).

Sub-RPE Deposition-Related Gene Expression Pattern Depends on Maturation State and Zinc Supplementation
Our hfRPE culture developed sub-RPE deposits even without photoreceptors and the supporting choriocapillaris (Supplementary Figure S1C). This allowed us to analyze the expression of genes potentially involved in the sub-RPE deposit formation process. We compiled lists of genes associated with various aspects of sub-RPE deposit formation and analyzed the changes in expression throughout cell maturation and zinc supplementation (Supplementary Table S11). Some genes belong to more than one gene list ( Figure 5).

Genelist 03
This contains cholesterol-metabolism-related genes ( Figure 5B, Supplementary Table S11 [40,41]). A total of 42 out of the 51 identified genes were expressed in hfRPE (Supplementary Table S11). The genes that were expressed higher at 2W were ABCG5, ANGPTL8, APOA1, CD36, LIPC, and STAR, while the genes expressed higher at 19W were APOC1, LDLRAP1, and NPC1. Among the genes expressed higher in less differentiated cells were CD36, LIPC, PCSK9, and STAR. In the more differentiated cells, the expressions of CYP27A1, LIPG, LPL, and PLTP were higher. PLTP, ANGPTL4, APOE, LRPAP1, VDAC2, and TSPO were significantly upregulated, and NPC2 was significantly downregulated in response to zinc supplementation. Interestingly, CYP27A1 showed significant upregulation at 2W and significant downregulation at 19W in response to zinc supplementation (Supplementary  Table S7).

Genelist 05
This contains genes that are related to pigmentation ( Figure 5E, Supplementary Table S11 [9,46,47]). Pigmentary abnormalities show strong correlation with sub-RPE deposit formation and the development of AMD, and we found that 19 out of the identified 21 genes were expressed in hpRPE (Supplementary Table S11). At 2W, we found no differentially expressed genes. At 19W, however, we found that AP3B1, AP3D1, BLOC1S6, HPS5, HPS6, and SLC24A5 were expressed higher. There were no highly expressed genes in less differentiated cells. In the more differentiated cells, the expression of OCA2 and SLC24A5 was higher. TYR, TYRP1, and DCT were significantly downregulated in response to zinc supplementation in our acute treatment (Supplementary Table S7).

Discussion
The RPE plays a pivotal role in maintaining the health of the retina, and changes in RPE function have been linked to the development and progression of AMD [48,49]. Optimal zinc balance is key for RPE function [50], and zinc deficiency contributes to AMD pathogenesis [51]. Based on these findings, it has been suggested that zinc supplementation can slow the progression of AMD [51,52], although the mechanism of this beneficial effect is not fully understood [53]. In this study, we used primary human fetal RPE cells and scRNA-Seq analysis to identify the transcriptomic changes and biologically plausible molecular pathways involved in the maturation of the RPE and the changes associated with zinc supplementation. The specific transcriptional changes and molecular pathways identified provide an improved understanding of RPE cell maturation and insight into how the function of RPE might be affected by acute zinc supplementation, which has relevance for the progression of AMD.

Study Rationale
Maturation of RPE cells is key to developing appropriate morphology, pigmentation [54], and production of key signature proteins that determine the function of these cells [55]. Different studies use a variety of sources to study RPE maturation and function, ranging from the immortalized ARPE-19 cells [56] to induced pluripotent stem-cell-derived RPE [57] and primary porcine [16] or human RPE [58]. As with all model systems, cellular models for RPE must replicate the in vivo situation as closely as possible. Recently we have shown that primary human fetal RPE cells develop the most critical features of native RPE, including the formation of pigmentation, tight junctions with high TEER values, and the expression of RPE signature genes and proteins [15,16]. Most importantly, the cells in culture can lay down sub-RPE deposits, a hallmark feature of AMD [15,16]. Despite demonstrating these in vivo-like features, the molecular signature for RPE maturation has not yet been fully explored. Previous studies have reported a variety of approaches to map molecular maturation. Earlier studies used microarrays [23,59] or bulk RNA sequencing. Most recently, a powerful tool capable of sequencing individual cells has been introduced. Single-cell RNA sequencing provides an unparalleled opportunity to identify cell heterogeneity [60]. Lidgerwood et al. [57] used pluripotent stem-cell-derived RPE to analyze transcriptomic changes after 1 month or 12 months in culture and analyzed these separately, then combined the data. In a subsequent study, the same group combined scRNA-Seq and proteomics in iPSC cells obtained from individuals with or without AMD to identify regulations in geographic atrophy [61]. Exciting opportunities are presented by scRNA-Seq studies using freshly isolated RPE from human eyes. RPE cells from both fetal and adult human eyes were analyzed in previous studies [31] and [32], respectively). Although both studies used a limited number of cells, they provide invaluable insight for cell-culture-based observations. In our study, we used primary fetal RPE cells that recapitulated features of RPE cells in vivo (Supplementary Figure S1). Despite their fetal origin, these cells developed sub-RPE deposits and varied pigmentation, suggesting that they recapitulate the hallmarks of AMD (Supplementary Figure S1) despite the relatively short time in culture (19W).

Heterogeneity of RPE Cells
The generation of scRNA-Seq data from a large number of cells allowed us to confidently determine that there is a significant degree of heterogeneity between the cells. A key observation was that some RPE cells could develop into more differentiated cells even after 2 weeks in culture, but even after 19 weeks, we still observed less differentiated cells (Supplementary Figure S1). Heterogeneity of RPE had been reported after multiple passages and over the years in culture [62] (Supplementary Figure S1), reflecting what had been reported for RPE in vivo [63][64][65] and in situ [62]. Despite the long-lasting heterogeneity, the melanosome precursor PMEL17 was expressed in both less and more differentiated cells. In fact, from the 19 pigmentation-related genes expressed in our cells, the only transcripts that showed elevated expression in the more differentiated RPE cells were OCA2 and SLC24A2 (Supplementary Figure S2(B1) and Supplementary Table S11, Figure 5E), suggesting that all cells could become pigmented [46,66].
COL1A1 was amongst the top transcripts in the less differentiated cells, and immunoreactivity of COL1A1 protein was able to distinguish the less differentiated cells from the more differentiated cells that express the RPE65 gene highly and are immunopositive for the RPE65 protein ( Figure 2C). Immunoreactivity to the COL1A1 protein gradually increased in the sub-RPE space with time in culture (Supplementary Figure S3), suggesting that the half-life of this extracellular matrix protein is long in our culture system. This increase in sub-RPE COL1A1 may correspond to the role this protein plays in forming the extracellular matrix of Bruch's membrane [67]. Other collagens were also expressed highly in the less differentiated cell population (Supplementary Table S4), reflecting their reported involvement in increased attachment and spread of RPE cells [68]. The only highly expressed transcript for collagen in the more differentiated cells was COL8A1 (Supplementary Table S4). The COL8A1 protein is a component of basement membranes in the eye and contributes to the formation of the basement membrane of RPE [21,69] and a genetic risk variant of AMD [70]. The findings on COL1A1 and RPE65 might be mechanistically important: the mature RPE cells (RPE65 expressing) could enable the performance of the visual cycle, while the less differentiated cells (COL1A1 expressing) can support the formation of ECM throughout life.

Transition from Less to More Differentiated RPE
As more and less differentiated cells are present at all three time points, we combined the scRNA-Seq data from the three time points and analyzed these datasets together, an approach different from a previous study [57]. This integrated approach helped us to identify a pseudotemporal trajectory of gene expression from less to more differentiated cells (Figure 3(A1)). This approach identified a well-defined main trajectory (Figure 3(A2)). The top genes with the highest score in the main trajectory were associated with regulating the visual cycle (RPE65, LRAT, TTR, RDH5) (Supplementary Table S6). Transcriptomic analysis of the bulk RNA isolated from RPE cells from aging human donor eyes recently reported a positive feedback mechanism between the upregulation of visual cycle genes and the accumulation of retinoid by-products [71]. As visual cycle-related bisretinoids are constituents of the accumulating lipofuscin in RPE [72], this upregulation could eventually lead to AMD-like pathogenesis [73] in this cell culture model. Indeed, there are ongoing clinical trials for visual cycle modulators as therapeutic options for AMD [74], and our cell culture model has the potential to serve as a preclinical tool for testing novel compounds.

Genes Involved in Transitioning RPE from Less to More Differentiated Cells
The genes associated with the main trajectory could be clustered into seven modules based on their transcriptional change along the pseudotemporal trajectory (Supplementary Table S5). The transcripts whose expression is transiently upregulated on the pseudotemporal trajectory likely represent the genes mediating the transition from the less to the more differentiated cells (Supplementary Table S5). These genes were associated with cellular and extracellular remodelling and metabolic pathways (Supplementary Table S6). Therefore, our data support the hypothesis that extracellular matrix remodelling of the Bruch's membrane could become a therapeutic target to combat RPE loss [75] due to topographic changes in the RPE-Bruch's membrane interface [68]. Alterations of the extracellular matrix may impact immune response as well as the secretion of pro-inflammatory cytokines, such as MCP-1 and IL-8 [68], and promote sub-RPE deposit formation [76][77][78][79][80]. Our data highlights potential molecular targets to achieve a regulation of this process.
Among the transiently expressed genes, we identified ID1 and ID3 (Supplementary Table S6). The corresponding helix-loop-helix (HLH) proteins form heterodimers with members of the basic HLH family of transcription factors, inhibiting DNA binding and preventing the formation of active transcriptional complexes [81]. ID proteins promote cell cycle progression and cell migration and restrict cellular senescence and the differentiation of a number of progenitor cell types [82,83]. Recent results indicate that the expression of ID family proteins may play an important role in regulating retinal progenitor cell proliferation and differentiation [84]. ID genes and proteins showed increased expression levels in the retina at embryonic and early postnatal stages and declined in the adult [84]. ID protein expression is silenced in many adult tissues but is re-activated in diverse disease processes [83,85,86]. ID proteins appear to play a crucial role in the angiogenic processes. It was proposed that inhibition of expression and/or function of ID1 and ID3 may be of therapeutic value for conditions associated with pathological angiogenesis [87]. In fact, the deletion of Id1/Id3 reduced ocular neovascularization in a mouse model of neovascular AMD [81]. In conclusion, drugs targeting ID1/ID3 could modulate RPE maturation and pathological changes in AMD.

Response to Acute Zinc Supplementation
Treatment with zinc has been reported to prevent progression to advanced AMD (for review, see [88]), at least partly due to a direct effect of zinc on the RPE [15,89,90]. In previous in vitro studies, we investigated long-term supplementation with zinc and found altered selective gene expression, protein secretion, and increased pigmentation and barrier function [15,17]. We identified several molecular pathways, such as cell adhesion/polarity, extracellular matrix organization, protein processing/transport, and oxidative stress response, involved in the beneficial effects of chronic zinc supplementation on the RPE. However, these studies could not address the complexity associated with cell heterogeneity and detailed temporal changes. We were particularly interested in exploring how zinc supplementation could affect the less and more differentiated cells in the short term to understand the potential to develop a more targeted intervention through supplementation.
To decipher the effects of acute zinc supplementation, RPE cells were treated with elevated zinc for 1 week following the protocols we used previously [17]. We found that acute zinc supplementation induced significant changes in gene expression in both shortand long-term cultures (Figure 4(B1-B4)) regardless of the temporal stage of the cells. We also identified 81 zinc-responsive transcripts ( Figure 4C) that were common amongst all groups. These transcripts were enriched in housekeeping genes and contained transcripts for metallothioneins, ribosomal protein, and ATP synthases (Supplementary Table S8), indicating that zinc affects the cellular homeostasis of the RPE, similar to that of other systems [91].
Apart from the shared genes, specific changes were associated with the more or the less differentiated cell groups at both 2W and 19W in culture ( Figure 4C). The four specific genes affected by short-term zinc supplementation in the less differentiated cell group (Supplementary Table S8) are genes linked to the integrity of Bruch's membrane (COL8A1) [70], epithelial-mesenchymal transition (KRT17) [92], phagocytic activity and the rescue of the RPE (MFGE8) [93,94], and activity of heparan sulfate (SULF1) [22], suggesting that zinc might influence interaction with the local extracellular environment. In the more differentiated cell group in the 2W cultures, zinc affected biological processes including extracellular matrix organization, cellular polarity, and visual processes (Supplementary Table S8) that are critical for supporting the photoreceptors [95].
At 19W in culture, zinc affected the less differentiated cells via modulating proteolysis, DNA replication and RNA transcription, and amino acid metabolisms (Supplementary  Table S8), probably to mitigate oxidative stress, one of the AMD-associated biological functions [96]. In the more differentiated cells at 19W in culture, zinc supplementation affected several metabolic pathways (Supplementary Table S8). Dysregulation of metabolic pathways is an important contributor to AMD pathophysiology [97]. This may directly explain the benefit of zinc supplementation in patients in the AREDS study [13,51,98]. Therefore, zinc supplementation has a multitude of effects on RPE, with some specific effects depending on cell differentiation and maturity. Identifying the specific molecular changes may help redefine treatment strategies based on zinc supplementation or nutritional interventions.

The Effects of Zinc on the Genes in the Pseudotemporal Trajectory
Earlier we identified 537 genes (Supplementary Table S5) in the main pseudotemporal trajectory (Figure 3(A2)). Zinc supplementation did not affect 240 genes (Figure 4(D1)). Of the remaining 297 genes, 16 were housekeeping genes (Supplementary Table S8; Figure 4C) associated with the mitochondrion, the activation of cytochrome-c oxidase and ubiquinone, and response to metal ions (Supplementary Table S9). This is in line with a previous observation that zinc supplementation can protect the RPE from oxidative-stress-induced cell death by improving mitochondrial function [89], and this could be behind the positive effect of zinc supplementation in the AREDS studies [13,51,98] or increased zinc intake through diet [34,99]. Metallothioneins (MT1F and MT1E) that belong to this group (Figure 4(D2), Supplementary Tables S7 and S9) are well-recognized mediators of zinc supplementation in the RPE [100] via mediating oxidative-stress-induced RPE damage [90] and differentiation of RPE [57].
The remaining 281 genes in the main trajectory (Figure 4(D1)) were associated with various biological processes including extracellular matrix organization, angiogenesis, collagen fibril organization, and visual perception (Supplementary Table S10). The composition of extracellular matrix has a profound effect on how the RPE attaches to the Bruch's membrane [101]. Thus, gene expression modification by zinc could directly affect sub-RPE deposit formation [76].
We also found that acute zinc supplementation upregulated the expression of transcriptional regulators ID1 and ID3, a finding that had not been reported before. In addition, in a previous study, we identified TGFB1 as a potential upstream regulator effect of chronic zinc supplementation [17]. In our current study, we found that TGFB1 expression was also upregulated by acute zinc supplementation. Therefore, we carried out an Upstream Analysis in Ingenuity Pathway Analysis (QIAGEN, Redwood City) for the 190 transiently expressed genes in the combined pseudotime-correlated groups 1 and 7 (Supplementary  Table S5). We identified a strong relationship for TGFB1 (p < 6.98 × 10 −19 ) and also for ID1 (p < 3.59 × 10 −5 ) and ID3 (p < 1.23 × 10 −3 ) as potential upstream regulators for a group of genes among the transiently expressed group. In fact, TGFB1 was an upstream regulatory element for ID1 and ID3 (Supplementary Table S12). A direct molecular link between ID1 and TGFB1 had already been suggested [102]. Therefore, the positive effects of zinc supplementation could be directly through TGFB1 signalling, which involves ID1 and ID3. The receptor of TGFB1, TGFBR1, is an AMD genetic risk variant [36], suggesting that these findings are directly relevant to further studies on AMD.

AMD-Specific Gene Expression Changes
Based on literature searches, we generated gene lists that have been shown to contribute to the pathological changes associated with AMD and we examined the effects of cell maturation and zinc supplementation on these genes ( Figure 5, Supplementary Table S11). Specific attention was paid to the activation complement system and lipidmetabolism-related genes, as these were the genetically most significantly associated pathways with AMD [36,103]. We also scrutinized genes associated with pigmentary changes and mineralization-associated genes due to their potential link with RPE function and/or sub-RPE deposit formation in AMD [5,45].
Not all genes involved in complement regulation were expressed in RPE cells ( Figure 5B, Supplementary Table S11). This is perhaps not surprising, as the local activity of the complement cascade is influenced by a complicated mix of local and systemic regulatory factors, which is altered in AMD retina [4,104,105]. However, some complement genes that were expressed in the RPE were affected by acute zinc supplementation, including CFH, C1R, CFHR1, and CLU ( Figure 5B, Supplementary Table S7). These transcriptomic changes are in line with our previous reports that zinc supplementation has a functional effect on CFH secretion [17] as well as oligomerization and activity [106], and zinc levels can regulate interferon gamma systematically, which, in turn, regulates expression of complement genes [107,108]. Apart from CFH, several complement proteins can bind zinc, and this binding alters their activity [109,110]. In addition, network analysis has highlighted elements of the complement regulation as potential targets for nutrient-affected pathways [111]. Finally, there is also clinical evidence that zinc supplementation can directly inhibit complement activation in AMD patients [104], suggesting that modulation of the complement system could be one of the ways that zinc supplementation affects the progression to AMD.
Of the 42 genes expressed in our RPE culture associated with cholesterol metabolism ( Figure 5C, Supplementary Table S11), ANGPTL4, LRPAP1, VDAC2, APOE, PLTP, NPC2, TSPO, and CYP27A1 were altered in response to acute zinc supplementation ( Figure 5C, Supplementary Table S7) [112][113][114]. These findings corroborate our previously reported effect of long-term zinc supplementation on lipid metabolism [17]. ANGPTL4 is a lipidinducible feedback regulator of LPL-mediated lipid uptake. However it is also a multifunctional cytokine, regulating vascular permeability, angiogenesis, and inflammation [115]. The systemic level of ANGPTL4 is associated with NV AMD [116]. Reportedly, this protein indirectly induces RPE barrier breakdown [117]. LRPAP1 is a chaperon protein, generally controlling the folding and ligand-receptor interaction expression of the LRP receptors [118]. Its role in RPE and AMD remains elusive. VDAC2 is a ceramide sensor integrated into the mitochondrial membrane and its function relates to regulation of mitochondrial apoptosis [119,120]. Increased ceramide levels affect non-polarized RPE cells found in late stages of AMD [121]. APOE, a lipophilic glycoprotein with a major role in lipid transport, is one of the many constituents of the sub-RPE deposits and has been associated with increased AMD risk [122,123]. PLTP is a phospholipid transfer protein and is one of the main players in lipid homeostasis in ApoB-containing particles and high-density lipoprotein metabolism. PLTP plasma levels are associated with AMD [124], but their potential role in drusen formation remains elusive. NPC2 is a cholesterol transporter, effluxing cholesterol out of late endosomes in RPE. The lack of this protein is associated with age-related maculopathies [125]. TSPO is a translocator protein that transfers cholesterol from the mitochondrial outer membrane to the mitochondrial inner membrane and also plays role in oxidative stress and inflammation. It was recently implicated as a highly relevant drug target for immunomodulatory and antioxidant therapies of AMD [126,127]. CYP27A1 is involved in the elimination of 7-ketocholesterol from RPE, a toxic product of cholesterol auto-oxidation, which accumulates in drusen [128,129]. In summary, the aforementioned affected gene expressions in response to zinc suggest that zinc has an impact on sub-RPE cholesterol accumulation, oxidative stress, inflammation, and angiogenesis via the regulation of lipid-membrane interaction, lipid transport, and the elimination of toxic lipid byproducts.
In our cultures, we found 80 RPE-expressed genes associated with mineralization ( Figure 5D, Supplementary Table S11). Out of these, we found that COL1A1, POSTN, CD63, LAMP1, BMP4, SLC20A1, SOX9, and BMP7 were altered in response to acute zinc supplementation ( Figure 5D, Supplementary Table S7). The POSTN gene encodes a secreted extracellular matrix protein that functions in tissue development and regeneration and a potential anti-fibrotic therapeutic target for NV AMD [130]. CD63 is involved in the regulation of cell development, activation, growth, and motility [131], and together with LAMP1, it plays a role in autophagy, exosome secretion, and drusen formation [132,133]. BMP4 has been implicated in the disruption of RPE cell migration and barrier disruption in NV AMD [134]. The protein encoded by SLC20A1 is a sodium-phosphate symporter involved in vascular calcification but not reported in association with RPE function or AMD [135]. SOX9 plays a key role in regulating visual cycle gene expression in RPE [136] but also plays a role in the prevention of calcification [137]. BMP7 is hypothesized to be critical for the differentiation of the retinal pigmented epithelium during development [138]. It also has been implicated in prevention of vascular calcification [139]. Zinc supplementation is reported to inhibit phosphate-induced vascular calcification [140], but, as our results indicate, it may also have a (indirect) role in the prevention of drusen calcification.
In our cultures, most pigmentation-related genes were detected and their expression level either remained constant or increased throughout the culture time ( Figure 5E, Supplementary Table S11). Only TYR, TYRP1, and DCT were altered in response to acute zinc supplementation ( Figure 5D, Supplementary Table S7). TYR, TYRP1, and DCT are key to the production of melanin [46] and pigmentary abnormalities show a strong correlation with sub-RPE deposit formation and development of AMD [9]. TYR catalyzes the production of melanin from tyrosine, in which L-DOPA is produced as an intermediate [141,142]. The function of TYRP1 is in the biosynthesis of melanin from tyrosine, whilst TYRP1 catalyzes the oxidation of 5-6-dihydroxyindole-2-carboxylic acid to an indole, whilst DCT catalyzes the conversion of L-dopachrome into 5-6-dihydroxyindole-2-carboxylic acid [142]. These events lead to the activation of GPR143 signaling and may initiate several downstream effects, such PEDF, VEGF secretion, and/or exosome release [46]. Since we found an influence of zinc on the expression of these pigmentation-related genes, and given the data from literature above, zinc might also have an influence on GPR143 signaling. Surprisingly, acute zinc treatment resulted in downregulation of the aforementioned genes, despite long-term zinc supplementation enhancing RPE pigmentation [15]. At the transcriptional level, long-term zinc supplementation significantly altered the expression of 18 out of the 21 pigmentationrelated genes (Supplementary Table S11, [17]), of which the majority were also downregulated, except for HPS5, HPS6, and LYST. These three upregulated genes are all related to intracellular trafficking, such as lysosomes and melanosomes [143][144][145].The negative effect of acute zinc supplementation on the gene expression of other pigmentation-related genes needs to be further investigated.

Conclusions
Primary hfRPE cultures that recapitulate the main phenotypes of aged RPE in vivo can help to dissect the molecular changes associated with RPE maturation and experimental manipulation, such as zinc supplementation. This cellular model provides an excellent platform for further preclinical studies to identify new treatment strategies for AMD. As reported in vivo, these cells retain a high degree of heterogeneity even after extended time in culture, which may help to understand the role of this heterogeneity in the human eyes. Identifying the transcriptional machinery, including transcriptional regulators ID1 and ID3, may help us to target pathways previously not considered for AMD. The data also show that the differentiation of RPE into cells that resemble those in vivo requires an extended time in culture, and experimental manipulation will need to consider this. The wide-ranging effects of zinc supplementation, from the regulation of housekeeping genes to very specific AMD-associated transcripts, build confidence that this intervention could indeed be a suitable intervention strategy to slow the progression to advanced-stage AMD, as suggested by the AREDS studies.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/cells12050773/s1, Figure S1: In vitro primary fetal RPE culture; Figure S2: Comparing gene expression patterns of more and less differentiated groups; Figure S3: Basal immunolabeling of COL1A1 protein increased in the RPE culture experiments; Figure S4: Gene expressional patterns of single-cell populations in the developmental and adult RPE ex vivo; Table S1: Unsupervised clusters marker genes; Table S2: Unsupervised clusters Gene Analytics; Table S3: RPE  specific genes; Table S4: Supplementary Table S4 More and Less differentiated cells marker genes;  Table S5: Pseudotime main trajectory genes; Table S6: Pseudotime main trajectory Gene Analytics;  Table S7: Acute zinc differentially expressed genes list; Table S8: Acute zinc differentially expressed Gene Analytics; Table S9: Acute zinc pseudotime overlap Gene Analytics 01; Table S10: Acute zinc pseudotime overlap Gene Analytics 02; Table S11: AMD Genelists; Table S12: Overlapping genes for  TGFB1 Upstream regulator table run

Data Availability Statement:
The data that support the findings of this study will be openly available in the GEO database public repository upon publication, which does not issue DOIs.